
==== Front
R Soc Open Sci
R Soc Open Sci
RSOS
royopensci
Royal Society Open Science
2054-5703
The Royal Society

39086829
rsos240557
10.1098/rsos.240557
1001100160197Ecology, Conservation, and Global Change Biology
Research Articles
Immigration allows population persistence and maintains genetic diversity despite an attempted experimental extinction
Immigration allows population persistence and maintains genetic diversity despite an attempted experimental extinction
https://orcid.org/0000-0002-5461-1339
Park Keon Young 1 Conceptualization Data curation Formal analysis Funding acquisition Investigation Visualization Writing – original draft kpark258@uwo.ca

Lucas Mel 1 Data curation Methodology Validation Writing – original draft mlucas26@uwo.ca

Chaulk Andrew 1 2 Data curation Methodology Validation achaulk@uwo.ca

Matter Stephen F. 3 Conceptualization Methodology Resources Writing – review and editing mattersf@ucmail.uc.edu

Roland Jens 4 Conceptualization Methodology Writing – review and editing jroland@ualberta.ca

Keyghobadi Nusha 1 Conceptualization Funding acquisition Methodology Project administration Supervision Writing – review and editing nkeyghob@uwo.ca

1 Department of Biology, Western University , London, Ontario N6A 5B7, Canada
2 Department of Biology, Memorial University of Newfoundland , St John’s, Newfoundland A1C 5S7, Canada
3 Department of Biological Sciences, University of Cincinnati , Cincinnati, OH 45221, USA
4 Department of Biological Sciences, University of Alberta , Edmonton, Alberta T6G 2E9, Canada
Electronic supplementary material is available online https://doi.org/10.6084/m9.figshare.c.7370663.

7 2024
31 7 2024 July 31, 2024
31 7 2024 July 31, 2024
11 7 24055708 4 2024 April 8, 2024
07 7 2024 July 7, 2024
08 7 2024 July 8, 2024
© 2024 The Author(s).
2024
https://creativecommons.org/licenses/by/4.0/ Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.

Widespread fragmentation and degradation of habitats make organisms increasingly vulnerable to declines in population size. Immigration is a key process potentially affecting the rescue and persistence of populations in the face of such pressures. Field research addressing severe demographic declines in the context of immigration among interconnected local populations is limited owing to difficulties in detecting such demographic events and the need for long-term monitoring of populations. In a 17-subpopulation metapopulation of the butterfly, Parnassius smintheus, all adults observed in two adjacent patches were removed over eight consecutive generations. Despite this severe and long-term reduction in survival and reproduction, the targeted populations did not go extinct. Here, we use genetic data to assess the role of immigration versus in situ reproduction in allowing the persistence of these populations. We genotyped 471 samples collected from the targeted populations throughout the removal experiment at 152 single nucleotide polymorphisms. We found no reduction in the genetic diversity of the targeted populations over time, but a decrease in the number of loci in Hardy–Weinberg equilibrium, consistent with a high level of immigration from multiple surrounding populations. Our results highlight the role of connectivity and movement in making metapopulations resilient to even severe and protracted localized population reductions.

immigration
; population rescue
; genetic maintenance
; genetic diversity
; artificial removal
; SNP
Alberta Conservation Association http://dx.doi.org/10.13039/100007583 Natural Sciences and Engineering Research Council of Canada http://dx.doi.org/10.13039/501100000038 Directorate for Biological Sciences http://dx.doi.org/10.13039/100000076
==== Body
pmc1. Introduction

In the current age of widespread habitat loss and degradation resulting from anthropogenic activities, populations of many species across the globe are experiencing, or are vulnerable to, severe declines in size [1,2]. A severe decline in the number of reproducing individuals increases demographic stochasticity, reducing the likelihood of a population’s persistence, especially for organisms that depend on a threshold population number for successful breeding or survival (Allee effect) [3,4]. In addition, the loss of genetic diversity caused either by an initial or temporary reduction in population size (i.e. a bottleneck) or continued exposure to high levels of genetic drift can lead to inbreeding depression and reduced evolutionary potential of populations, exacerbating the risk of population extinction [5,6]. However, the effects of demographic and genetic stochasticity in small populations can potentially be countered through various forms of rescue, frequently involving the immigration of individuals into the population either naturally or through translocation by humans [7].

Three primary forms of population rescue are generally recognized: evolutionary, demographic and genetic. Evolutionary rescue is the recovery of a population by adaptation to the conditions that caused the initial decline [7,8]. In the strictest sense, evolutionary rescue occurs by natural selection acting on standing genetic variation and does not involve immigration from other populations [7,9]. However, a broader definition of evolutionary rescue can include selection acting on novel genetic variants introduced by immigrants that are particularly well suited to the altered or degraded environmental conditions [10,11]. Although more often reported under controlled laboratory conditions [12], examples of strict-sense evolutionary rescue do exist in nature, especially in organisms with rapid generation times, such as insects [13] and rodents [14]. Demographic rescue involves an influx of individuals through immigration that acts to directly augment population size and protect the population from stochastic demographic events and Allee effects [15,16]. Finally, genetic rescue refers to an increase in population fitness caused by the genetic contributions of immigrants that reduce the extent of inbreeding and inbreeding depression [17–19]. A broader definition of genetic rescue can include a role for genetic variation introduced by migrants in contributing to evolutionary potential, and therefore, there is some overlap in evolutionary and genetic rescue in their broadest senses [7]. Regardless, immigration probably plays a central role in many cases of population rescue [7,18], whether occurring through artificial translocation or natural dispersal.

Several studies demonstrate that connectivity among local populations, which increases the likelihood of between-population dispersal and immigration, is a key factor in the recovery or maintenance of local population numbers in the face of demographic crashes or bottlenecks [20–22]. Furthermore, populations rescued by immigration have shown increases in both genetic diversity and fitness [8,22,23]. An extreme example of genetic rescue is the rapid, albeit temporary, recovery of the isolated wolf population on Isle Royale, triggered by the dispersal of a single individual from the mainland population [24]. Overall, the positive effects of both demographic and genetic rescue have been established in various natural and experimental systems [18,25]. Furthermore, the importance of immigration to rescue processes and to the persistence of populations and metapopulations more generally, is widely recognized [26,27]. However, the ability of natural populations to persist under sustained conditions that severely elevate mortality or reduce reproduction and the potential for immigration to contribute to persistence under such conditions are less explored. It has been suggested that detailed, empirical assessments of the prevalence and power of immigration in maintaining populations lag behind our theoretical understanding [28].

Given the multitude of processes simultaneously contributing to rapid and extreme environmental changes across the globe, including climate change, land use change and invasive species, many populations may be confronted by sustained and severe pressures that reduce their abundance. It is, therefore, important to understand the key factors contributing to persistence in the face of such pressures to be able to inform management and conservation actions [29]. Here, we assess the role of immigration in allowing populations to persist despite a continuous and extreme elevation of mortality, along with reduction of reproductive output sustained over multiple generations. We genotype samples collected across an eight-generation-long experiment involving the continual removal of individuals from target populations within a natural metapopulation system of an alpine butterfly, the Rocky Mountain Apollo (Parnassius smintheus) [30].

The study metapopulation is located on Jumpingpound Ridge in the Rocky Mountains of Kananaskis Country, Alberta, Canada, where 17 subpopulations in alpine meadow habitat are separated by varying amounts of intervening forest [31]. Starting in 2001 and ending in 2008, the south-eastern portion of the network (figure 1) was subject to a long-term experiment in which all adults observed in two adjacent patches (P and Q, figure 1) were captured and removed each year, to examine effects of a simulated local extinction on the population dynamics of the surrounding patches [30,32,33]. Specifically, the aim of the removals was to mimic or induce an extinction event in the target patches, to test the hypothesis that their extinction would reduce immigration into neighbouring patches and thereby lead to reduced population sizes, increased extinction risk, and decreased growth in the surrounding populations proportional to their connectivity to the target patches [33]. While the removals aimed to induce local extinctions in patches P and Q, these target populations did not go extinct, nor did the likelihood of extinction in surrounding patches increase [30,33]. However, the reproductive output of the focal patches P and Q was reduced by more than 75% by the removals [30].

Figure 1. Map of the Jumpingpound Ridge study system located in Kananaskis Country, Alberta, Canada showing outlines of meadow habitat patches (above approximately 2000 m). Each habitat patch is labelled with a unique letter label. The dotted circle contains the section of the population network studied here, with patches P and Q (meadows filled in dark) being the target of artificial removals and other patches in the dotted circle (meadows filled in light grey) being within potential dispersal distance of the butterfly, Parnassius smintheus (i.e. within 2 km of patches P and Q).

Map of the Jumpingpound Ridge study system located in Kananaskis Country, Alberta, Canada

The populations in the target experimental patches persisted, therefore, despite the yearly removal of all observed adults over eight generations [30]. Two not mutually exclusive mechanisms may explain the survival of these populations. First, it is possible that not all adults in the two populations were removed each year and the offspring of those that were left behind experienced a higher per capita survival and reproductive success thanks to release from density-dependent factors [34]. Second, the removals may have created sinks that were recolonized by immigration from surrounding populations in the network. Here, we test these competing hypotheses by tracking genetic diversity of the target populations through the course of the experiment across the 8 years. If persistence occurred primarily through in situ reproduction from within the targeted populations, this would have created successive bottleneck events and we would therefore expect the local populations to decline significantly in genetic diversity over time [35,36]. Furthermore, we would expect a substantially greater reduction in allelic diversity relative to heterozygosity, as the former is more sensitive to bottleneck events [36]. However, to the extent that populations in the targeted patches P and Q persisted primarily owing to rescue through immigration from multiple different source populations, we would expect them to better maintain levels of genetic diversity despite the continuous removal of individuals. In this case, we would also predict decreasing Hardy–Weinberg equilibrium (HWE) over time owing to the continuous pooling of immigrants from genetically differentiated source populations (i.e. Wahlund effect) [37]. We, therefore, use genetic data collected over multiple generations to assess the relative roles of in situ reproduction versus immigration in allowing population persistence in the face of the attempted experimental extinctions.

2. Methods

2.1. Study organism, system and sample collection

The Rocky Mountain Apollo butterfly inhabits primarily high-altitude alpine meadow habitats throughout the Rocky Mountains in North America. The larvae’s main host plant is the perennial succulent, Sedum lanceolatum [38], though Rhodiola integrifolia is also occasionally used. The species is univoltine, completing a single generation each year. The larvae pupate after completing five instars, from which the adults emerge around July to August to mate and oviposit [39]. This annual adult flight season is the only time individuals can disperse from their natal habitat patch to another, as the larvae are not known to be able to move between patches. Dispersal is generally limited to neighbouring habitat patches; average recapture distances of marked adults within a flight season are only about 130 m for both males and females, and less than 10% of individuals are typically recaptured outside of the habitat patch in which they were originally captured and marked [31].

The Jumpingpound Ridge metapopulation occurs within a network of 17 patches of alpine meadow habitat ranging in size from 0.2 to 22.7 ha [31]. This system has been the subject of long-term population monitoring and genetic sample collection since 1995, with mark-recapture performed during the annual flight season to monitor population size and movement [31]. Populations in this system are genetically structured and normally display isolation by distance [40]. However, occasional network-wide bottlenecks in population size, thought to be driven by unfavourable over-wintering conditions, temporarily lead to increased genetic differentiation and loss of isolation by distance, as different alleles persist in different local populations; this differentiation declines rapidly following the bottleneck events and isolation by distance re-establishes as gene flow among the populations redistributes the surviving genetic variation [41,42]. The removal experiment on patches P and Q was conducted from 2001 to 2008 to investigate the effects of severe and continuous population reduction in patches P and Q on the population dynamics in surrounding, neighbouring patches [30,33]. Through the course of the experiment, patches P and Q were surveyed every 1–3 days during the flight season each year, and all observed butterflies were captured by hand netting and removed from the site. In total, 4830 butterflies were removed from P and Q over the 8 years. The combined number of individuals removed from both P and Q was highest in the first year of the experiment (approximately 1200) and lowest in 2003 (<100), a year in which the entire population network experienced a severe demographic bottleneck (fig. 2 in [30]). The removed individuals were labelled with the date and location of capture, and then stored, dried and pinned in a collection housed at the University of Cincinnati.

Figure 2. Changes in population genetic metrics over time within local populations in patches P (dark solid lines) and Q (light solid lines), through the course of an 8 year experiment in which all adults observed in each patch were removed each year from 2001 (year 1) until 2008 (year 8): (a) mean allelic richness, (b) expected heterozygosity, (c) percentage of single nucleotide polymorphism (SNP) loci not in HWE, and (d) percentage of SNP locus pairs in linkage disequilibrium. The dotted lines (red) represent the best-fit linear mixed model for each genetic response variable as a function of time (year of the experiment) as the fixed effect, with patch (P or Q) as a random effect.

Changes in population genetic metrics over time within local populations.

We sampled from this collection of pinned specimens from patches P and Q, removing one leg per individual across each year of the removal experiment. For years in which the number of captured individuals in a patch was lower than 30, we sampled all available individuals. Otherwise, we sampled between 30 and 60 individuals per patch per year. Samples from 2001 represent the basal or initial state of the populations, as they would not yet have responded to the experimental manipulation. In addition, leg tissues from 42 frozen whole butterflies, removed during the 2008 flight season were also used; in combination with pinned samples from 2008, these represent the final year of the experimental removals. In total, we sampled 471 individuals from patches P and Q across 8 years (table 1). The forest separating patches P and Q forms only an incomplete barrier (figure 1) and individuals can move readily between these patches. To be consistent with previous studies in this system, including the removal experiments, we analysed samples from patches P and Q separately.

Table 1. Sample size and estimates of basic population genetic metrics (based on 152 polymorphic SNP loci) for samples from patches P and Q for each year of the removal experiment (2001–2008). (MAR, mean allelic richness; HE, expected heterozygosity; ‘loci out of HWE’, the proportion of loci significantly out of HWE.)

patch	year of collection	sample size	MAR	HE	loci out of HWE	
P	2001	24	1.76	0.24	0.027	
P	2002	20	1.75	0.21	0.011	
P	2003	28	1.78	0.21	0.037	
P	2004	29	1.77	0.22	0.059	
P	2005	27	1.77	0.24	0.011	
P	2006	30	1.77	0.24	0.080	
P	2007	29	1.79	0.24	0.054	
P	2008	43	1.78	0.24	0.064	
Q	2001	21	1.71	0.22	0.021	
Q	2002	30	1.77	0.20	0.037	
Q	2003	14	1.78	0.22	0.005	
Q	2004	31	1.79	0.22	0.043	
Q	2005	28	1.78	0.24	0.043	
Q	2006	26	1.76	0.23	0.043	
Q	2007	30	1.79	0.24	0.064	
Q	2008	61	1.78	0.24	0.123	

To assess relationships of individuals removed from the target patches to populations that could have been potential sources of immigrants, we also genotyped samples from surrounding patches most likely to contribute migrants to P and Q. Based on a priori knowledge that individual P. smintheus travel approximately 130 m on average, to an observed maximum of 2 km within one flight season [31], we determined that the five neighbouring patches (M, N, O, R and S) within 2 km of P and Q contained potential source populations (figure 1). This is further substantiated by previous estimates of patch connectivity indicating that only patches M, O, R and S have sufficient connectivity to patches P and Q to be influenced by any changes in the population sizes of the latter [41]. Patch N is a small patch that harbours a very small local population (with population size typically in the single digits and subject to frequent extinction) [43]. Patch N is therefore an unlikely source of immigrants; furthermore, because of its low population size, no or very few samples are available from that patch in most years. While connected by movement of adults [31], these surrounding populations are semi-independent of the experimental populations, and of each other, as reflected in patterns of genetic differentiation [40–42] and a degree of asynchrony in their population dynamics [33].

Non-lethal samples of wing tissue were removed from individuals in surrounding populations during the course of mark-recapture studies, as described previously [31]. Briefly, a small (approximately 2 mm2) piece of wing tissue was removed from newly marked individuals and stored immediately in 100% ethanol until DNA extraction. With the exception of samples collected in 1995 [40], we did not have tissue samples available from the surrounding populations prior to 2005 as we were not regularly engaging in wing clip sampling before that time [42]. Here, we include genetic data from surrounding populations M, O, R and S, which are the only likely sources of direct immigrants to P and Q, in 2005 and 2008 (the final year of the experiment); these correspond to years in which genetic diversity and differentiation across the ridge were previously analysed using microsatellite data [42].

2.2. DNA extraction and single nucleotide polymorphism genotyping

We extracted DNA from leg and wing tissues using the DNeasy Blood and Tissue Kit (Qiagen, Germantown, MD). Each tissue sample was first homogenized in a microcentrifuge tube containing lysis buffer using a disposable pestle, before being incubated with proteinase K at 56°C for approximately 18 hours. The purified DNA was eluted from the Qiagen spin columns with 200 μl of purified water, preheated to 37°C. We genotyped each individual at a panel of single nucleotide polymorphism (SNP) loci using the Agena iPlex Gold MASSarray platform (Sequenom, San Diego, CA), which can perform robust genotyping assays using samples of low DNA concentrations (<0.01 ng μl−1).

Our SNP panel contained 171 unlinked SNP loci and was developed based on a double-digest restriction-site associated DNA sequencing (ddRADseq) dataset for the species [44]. Briefly, multiple ddRADseq libraries were constructed from 80 to 100 individuals each (for a total of 501 individuals, sampled from Jumpingpound Ridge and the surrounding region [45] and not overlapping with the individuals removed from patches P and Q) by digesting genomic DNA using the restriction enzymes NlaIII and EcoRI, then selecting and amplifying fragments between 200 and 500 bp and sequencing the libraries on an Illumina HiSeq 2500 sequencer. We identified SNPs from the resulting sequences using the program STACKS [46] with the following parameters: minimum stack depth (m) = 3, mismatches allowed between putative catalogue loci (n) = 3, mismatches allowed between putative alleles (M) = 2, mismatches allowed to align secondary reads (N) = 4, and a maximum allowed missing data of 50%. This resulted in a ddRADSeq dataset of 8814 SNP loci. From the ddRADSeq sequences, we filtered SNPs with at least 40 bp upstream and downstream flanking sequences to allow space for designing primers for the iPlex Gold MASSarray assay, and used the software PLINK to identify and remove statistically linked loci. To identify SNPs that might putatively be functional, we aligned the filtered ddRADSeq sequences to a transcriptome of thorax tissue of adult butterflies captured during flight [47] using Magic-BLAST [48]. We defined fragments with at least a 90% sequence match as putatively expressed loci. Potential functions of these putatively expressed loci were identified using Trinotate [49], and 35 loci were chosen to reflect a variety of cellular functions such as metabolism, transcription regulation, transmembrane proteins, protein modification and insect development. An additional 130 loci having both less than 10% sequence match to the transcriptome and a greater than 90% sequence match to a shotgun-sequenced genome [50] were chosen at random to include in the final SNP panel. Finally, an additional six non-synonymous SNPs were included from the coding region of phosphoglucose isomerase (Pgi) [47], a gene associated with dispersal and flight metabolism in other insects, including the Glanville fritillary butterfly [51,52].

2.3. Genetic variation in the target populations over time

We estimated allele and genotype frequencies, and metrics of genetic diversity and differentiation, in each year of the removal experiment, separately for patches P and Q. All analyses were performed in the R statistical platform [53].

2.3.1. Genetic diversity

We estimated allelic richness, rarefied to 14 individuals (the lowest sample size, in patch Q in year 3), using the package ‘PopGenReport’ [54] and conducted a Wilcoxon test for the significant differences in mean allelic richness (averaged across loci) between year 1 of the experiment (i.e. original population unaffected by the yearly removals) and year 8 (final year of removals). We estimated expected heterozygosity and tested for a significant difference in mean expected heterozygosity across loci between the first and last year of the experiment using the package ‘adegenet’ [55]. We also fitted each of the mean allelic richness and mean expected heterozygosity in a separate linear mixed model in the ‘nlme’ package [56] with year as a numeric predictor (i.e. reflecting time since the start of the removal experiment) and patch identity (P or Q) as a random factor.

2.3.2. Hardy–Weinberg and linkage

We tested for HWE at each locus in each year using the package ‘pegas’ [57], with the number of Monte Carlo replicates set to 1000 and α = 0.01 to determine significance. Linkage disequilibrium tests were conducted for each pair of loci in each year with Genepop 4.7.5 [58], with the number of dememorizations set to 10 000 and α = 0.01 to determine significance. We fitted the proportion of SNPs out of HWE and the proportion of locus pairs in linkage disequilibrium in separate linear mixed models in the ‘nlme’ package [56] with year as a numeric predictor (i.e. reflecting time since the start of the removal experiment) and with patch identity (P or Q) as a random factor.

2.3.3. Genotype and allele frequency changes across the removal experiment

To capture changes in allele frequency distributions across all loci over time, we used the R package ‘hierfstat’ [59] to estimate pairwise F ST between samples from consecutive years of the experiment (i.e. year 1–2, 2–3, etc.), as well as between year 1 and each subsequent year (i.e. year 1–2, 1–3, etc.), within each of patches P and Q.

To test for changes in genotype frequencies of individual loci over time within each of patches P and Q we used a series of multinomial logit models, implemented in R using the package ‘mclogit’ [60]. We coded genotypes at each locus in a trinomial form, with the minor (less common) homozygote, heterozygote and the major (more common) homozygote assigned a value of ‘0’, ‘1’ and ‘2’ respectively. Owing to the trinomial form of the multinomial model, all assayed SNPs that were monomorphic or lacked the minor (less common) allele homozygote in the dataset were also removed, such that we assessed changes in genotype frequencies at 145 SNPs in patch P and 135 SNPs in patch Q. The observations of each of the three possible genotypes were then modelled as a function of time (years), separately for each SNP and in each patch. Significance values were adjusted for multiple tests using the Benjamini–Hochberg correction [61].

2.4. Relationships to putative sources of immigrants

2.4.1. Pairwise F ST relative to surrounding patches

We estimated pairwise F ST between populations in each of the experimental patches (P, Q) and each of the surrounding patches that are potential sources of immigrants (M, O, R and S) in 2005 (the first year of the experiment for which genotype data were available for surrounding patches) and 2008 (the final year for the experiment) using the ‘hierfstat’ package [59].

2.4.2. Population assignment of potential dispersers

We used the approach implemented in Geneclass 2.0 [62], with the Rannala & Mountain [63] Bayesian method of assignment, to assign individuals removed from patches P and Q to putative patches of origin from the set of neighbouring, potential source populations. Importantly, these assignment analyses were not intended to test or validate immigration as the main factor allowing the persistence of populations P and Q. Instead, once we had inferred a primary role of immigration through temporal changes in genetic diversity and HWE, the assignment tests were conducted to examine potential sources of those immigrants. In particular, our goal was to assess which neighbouring patches provided the majority of immigrants, and if the source of immigrants was stable or variable over time. These analyses make the assumption that all individuals captured and removed from patches P and Q are first-generation migrants from surrounding patches; while some removed individuals may have been born in P and Q, our analysis of genetic diversity over time suggested that a large proportion removed each year were likely to be first-generation immigrants. We only assigned individuals removed from P and Q in 2005 and 2008, years in which genotype data from neighbouring populations were also available, and we only considered assignments to surrounding populations using data from the same year (i.e. samples removed in 2005 from P and Q were assigned to surrounding populations using 2005 data from those populations). To assess the quality of assignments provided by our SNP dataset, we also assigned all individuals sampled in the potential surrounding source populations to the same set of populations (M, O, R and S), separately for 2005 and 2008; the extent to which individuals are assigned back to their ‘own’ population provides an index of the quality of assignments in a given dataset [62].

3. Results

3.1. Genetic variation in the target populations over time

Of the 171 SNPs targeted for genotyping, assays for three SNPs failed altogether, seven SNPs failed across more than 70% of samples, and nine SNPs were monomorphic across all samples in all years from both patches P and Q. These 19 SNPs were removed, and we used the genotype data from the remaining 152 polymorphic SNPs for further analyses. Of the 477 individuals genotyped from P and Q, six individuals were successfully genotyped at less than 80% of all SNPs and were also removed from the dataset. The total genotype failure rate across all remaining 471 individuals and 152 SNP loci used for further analyses was 4.1%.

3.1.1. Genetic diversity

Populations in both patches P and Q showed a slightly increasing trend in mean allelic richness over the course of the removal experiment, with only minor oscillations (figure 2a ). The mean allelic richness in a given year ranged from 1.764 to 1.827 (s.e. = 0.005) for patch P and 1.764 to 1.816 (s.e. = 0.005) for patch Q. A Wilcoxon test comparing the first and final years of the experiment (i.e. years 1 and 8) indicated no significant change in mean allelic richness through the course of the entire experiment for either P (p = 0.359) or Q (p = 0.543). However, the year was a significant predictor of mean allelic richness across the entire dataset with a positive slope coefficient indicating a trend of increasing allelic richness over time (p = 0.023, β = 0.004 ± 0.001; figure 2a ).

The mean expected heterozygosity was also mostly stable throughout the experiment with a slightly increasing trend (figure 2b ). The mean expected heterozygosity in a given year during the experiment ranged from 0.219 to 0.230 (s.e. = 0.001) for patch P, and 0.207 to 0.229 (s.e. = 0.002) for patch Q. Testing for a significant difference between the expected heterozygosity of the first and final years (i.e. years 1 and 8) of the experiment, we found no significant difference in patch P (p = 0.798) but did find one in Q (p = 0.014), where the mean expected heterozygosity increased from 0.207 in year 1 to 0.224 in year 8. The year was not a significant predictor of mean expected heterozygosity across the entire dataset (p = 0.209, β = 0.76 × 10–3 ± 0.57 × 10–3).

3.1.2. Hardy–Weinberg and linkage

The proportion of SNPs out of HWE in each patch fluctuated throughout the experiment but showed a strong increasing trend towards the end of the experiment, particularly in patch Q (table 1; figure 2 c ). In the initial year of the experiment (i.e. representing the ‘pre-removal’ conditions), 3.11% and 1.86% of loci were out of HWE in patches P and Q, respectively; by the final year, these values increased to 7.45% and 13.04%, respectively. The year was a significant positive predictor of the number of loci out of HWE overall (p = 0.001, β = 1.040 ± 0.249).

The proportion of SNP pairs in linkage disequilibrium also fluctuated over time but displayed a strong increasing trend towards the end of the experiment in both patches P and Q (figure 2d ). In the initial year of the experiment (i.e. representing the ‘pre-removal’ conditions), approximately 0.028% and 0.009% of pairs of loci were in linkage disequilibrium in patches P and Q, respectively; by the final year, these values increased by approximately three and ninefold to 0.095% and 0.077%, respectively. Year was a significant positive predictor of the number of pairs of loci in linkage disequilibrium (p < 0.001, β = 0.756 × 10−3 ± 0.0001).

3.1.3. Genotype and allele frequency changes across the removal experiment

Year-to-year (i.e. year 1–2, 2–3, etc.) pairwise F ST estimates were low to moderate in both patches, ranging between 0.003 and 0.012 (mean = 0.007 ± 0.001) in patch P, and between 0.003 and 0.01 (mean = 0.006 ± 0.0008) in patch Q. Pairwise F ST estimates between year 1 and each subsequent year of the experiment were similarly low, ranging from −0.004 to 0.008 (mean = 0.002 ± 0.0014) in patch P, and from −0.002 to 0.005 (mean = 0.002 ± 0.0008) in patch Q (electronic supplementary material, figure S1). Pairwise F ST estimates between years 1 and 8, the first and final years of the study, were also low (P: F ST = −0.007; Q: F ST = −0.005). None of the pairwise F ST estimates, including comparisons between the first and final years of the experiment, were significantly different from zero (i.e. the 95% confidence intervals (CIs) bracketed zero), with the exception of the comparison between years 1 and 6 in patch P (F ST = 0.008; 95% CI = 0.0008–0.017).

After controlling for a false discovery rate, we did not detect any SNPs that showed significant change in genotype frequency across years of the experimental removals.

3.2. Relationships to putative sources of immigrants

3.2.1. Pairwise F ST relative to surrounding patches

Pairwise F ST between each of the experimental patches (P, Q) and immediately surrounding patches (M, O, R and S) were low to moderate within each year tested and declined in value overall from 2005 to 2008. In 2005, F ST estimates ranged between 0.038 and 0.079 for all pairwise comparisons involving patch P (mean = 0.050 ± 0.010) and between 0.037 and 0.074 for all pair-wise comparisons involving patch Q (mean = 0.050 ± 0.008). In 2008, the F ST values ranged between 0.026 to 0.032 for all pairwise comparisons involving patch P (mean = 0.030 ± 0.001) and between 0.026 to 0.031 for all pairwise comparisons involving patch Q (mean = 0.030 ± 0.001). Differentiation among populations in the non-experimental patches themselves (M, O, R and S) declined from 2005 to 2008, with mean pairwise F ST among only those populations at 0.046 ± 0.016 in 2005 and 0.006 ± 0.002 in 2008.

3.2.2. Population assignment of potential dispersers

Individuals removed from the experimental patches could be assigned to a given neighbouring patch with generally high confidence. The assignment score (likelihood of a given population as the source of the individual relative to the sum of likelihoods for all populations, expressed as a percentage; [62]) to the most likely potential source population, averaged among all individuals collected from P and Q (in both 2005 and 2008), was 84.29 ± 2.26% and 82.47 ± 1.61%, respectively, compared with 14.51 ± 1.96% and 13.33 ± 1.30% for the second most likely source population. For the samples collected from the surrounding patches M, O, R and S, all individuals in 2005 and all but two individuals in 2008 were assigned to their 'own' patch with the highest likelihood (in 2008, one individual sampled in O was assigned to M with an assignment score of 95.5% and one individual sampled in S assigned to R with an assignment score of 75.7%). The mean assignment score for individuals from the surrounding patches to their own patch (i.e. where they were sampled), which provides an index of the quality of the assignments [62]), was 97.5% in 2005 and 95.1% in 2008.

In 2005, the majority of individuals removed from patch P were assigned with the highest likelihood to neighbouring patches M (37.0% of individuals) and S (51.9% of individuals), with 11.1% assigned to R and no individuals assigned to patch O (out of 27 individuals sampled from P in 2005). However, the relative contribution of inferred immigrants to patch P from patches M and S declined in 2008 (proportion of sampled individuals assigned with the highest likelihood to M: 11.6%, and to S: 32.6%), while patch O and R’s contributions increased substantially with 32.6% and 23.3% of individuals assigned, respectively (figure 3).

Figure 3. Proportion of samples collected in (a) patch P, and (b) patch Q, in the 2005 (dark bars) and 2008 (light bars) flight seasons, assigned to each of four surrounding populations that are potential sources of immigrants (in neighbouring patches M, O, R and S, located within 2 km of the target patches). Each individual removed from patch P or Q was assigned to the most likely of the potential source populations based on genotypes at 152 polymorphic SNP loci, using the Rannala & Mountain [63] Bayesian method of assignment implemented in Geneclass 2.0.

Proportion of samples collected in (a) patch P, and (b) patch Q, in the 2005

About half of the individuals removed from patch Q in 2005 were assigned with the highest likelihood to patch M (46.4% of individuals). In 2005, only one individual captured in patch Q was assigned with the highest likelihood to O, out of 28 total individuals. However, patch Q also displayed a large apparent decrease in inferred immigration from patch M, and an increase in immigration from patch O, from 2005 to 2008. The proportion of individuals captured in Q that assigned with highest likelihood to patch M more than halved from 2005 to 2008, while those assigned to patch O increased by almost 20% points, making the inferred contribution of immigrants to experimental patch Q fairly even across the source populations in 2008, compared with 2005. (figure 3).

4. Discussion

The continual removal of all observed adult individuals from two populations over eight consecutive generations did not result in the extinction of those populations [30,33]; we show here that those removals also did not lead to a reduction of genetic diversity in the populations. Our results support immigration from neighbouring populations as the main process allowing the persistence despite the experimental removals and highlight the resilience of this metapopulation system to severe, localized reductions in survival and reproduction.

Any isolated population experiencing the high level of mortality and reduced reproductive output induced by the experimental removals consistently over several generations [30], would be expected to show a significant decline in genetic diversity [34,64]. If such an isolated population was not driven to extinction and persisted primarily through in situ reproduction, the resulting bottleneck events over generations would cause a significant decline in allelic diversity, in particular [36,65]. However, mean allelic richness and expected heterozygosity for both experimental patches P and Q did not decline across the years, with allelic richness even showing a moderate but significantly increasing trend over time (figure 2a,b ). The maintenance of genetic diversity, especially allelic richness, supports the alternative mechanism for the persistence of these populations: consistent demographic recovery by incoming dispersal from other populations in the network.

Further supporting the immigration hypothesis is the significant increase in the proportion of loci out of HWE within the experimental populations over time (figure 2c ); a continuous, yearly influx and pooling of dispersers from multiple partially differentiated source populations would be expected to lead to increasing deviations from HWE through the Wahlund effect [37]. By contrast, successive bottlenecks would not be expected to lead to increasing deviation from HWE (i.e. a mismatch between observed allele and genotype frequencies) with each generation, as long as the probability of an individual being removed from the population was random with respect to genotype and matings among surviving adults were also largely random. Consistent with an increasing number of loci out of HWE, we also observed a substantial increase over time in the numbers of pairs of loci in linkage disequilibrium (figure 2d ), although this pattern would be expected under both hypotheses of mixing of immigrants from multiple gene pools and successive bottlenecks [66].

Additional evidence in support of the immigration hypothesis is the lack of differentiation of successive generations in each of the experimental patches from their starting state (i.e. F ST not significantly >0 between even the first and final year of the experimental removals). In the absence of significant immigration, we would expect successive local bottlenecks to lead to rapidly increasing genetic distance between the initial population state and each subsequent generation. By contrast, we would expect a continual source of immigrants from a given set of surrounding populations to result in much more stable allele frequencies over time, as was observed here, with perhaps an initial increase in genetic distance only after the first year of removals.

We also observed a decline in genetic differentiation over time (i.e. in 2005 versus 2008) between each of the populations in P and Q and their surrounding populations. This pattern is concordant with a previously described decline in genetic differentiation among all populations on the ridge over the same time interval, associated with synchronous fluctuations in population size [41,42]. Specifically, the entire ridge experienced a severe collapse in adult population numbers in 2003, followed by a recovery of population sizes over the following years. Previous work has shown that genetic differentiation among populations on the ridge increased after the initial collapse (as measured in 2005) as different alleles survived in different local populations, but rapidly declined over subsequent years (as measured in 2008) through gene flow [41,42]. The increased genetic similarity between each of P and Q and their surrounding populations between 2005 and 2008 may partly reflect the changes in genetic differentiation across all the non-experimental populations on the ridge as result of population cycles. Importantly, however, the differentiation between each of P and Q and their surrounding populations declined only slightly from 2005 to 2008 (by approximately 40% on average), while the differentiation among all the surrounding populations to each other declined considerably more (almost eightfold; electronic supplementary material, table S1). Therefore, although allele frequencies in the surrounding populations (and relationships of those populations to each other) changed considerably from 2005 to 2008, the relationships between each of P and Q to those surrounding populations were much more stable over the same time. This observation is also consistent with a hypothesis of high levels of sustained immigration over that period from all the surrounding populations into patches P and Q.

Importantly, we do not suggest that no individuals were left behind in patches P and Q after the yearly removals or that there was zero in situ reproduction. Indeed, it is very likely that at least some adults present in the populations were left behind each year of the experiment and that they did reproduce locally. However, our results do indicate that immigration from surrounding populations was the dominant process, relative to in situ reproduction and population growth, that allowed populations to persist and maintain the sizes of those populations through the experiment.

There are several examples of populations showing no or minimal loss of genetic diversity following observed demographic bottlenecks as a result of immigration from external sources [20,21,67]. Furthermore, previous work on P. smintheus in the Jumpingpound Ridge system found that patch connectivity, which facilitates immigration, was strongly associated with faster recovery of allelic diversity within local patches after the ridge-wide demographic collapse in 2010 [68]. Our current study highlights the role of immigration in maintaining genetic diversity and allowing population persistence in the face of even extreme and lengthy demographic pressures, further supporting the importance of patch connectivity for both demographic and genetic resilience. Although other studies have demonstrated the ability of immigration to counter high local mortality in populations, including from lethal culling of managed populations [69,70], our study represents an extreme case of severe population reduction (all observed adults removed from the population) occurring long term (over eight consecutive generations), and represents a particularly strong test of the ability of immigration to maintain local populations.

With immigration implicated as the main driver of population persistence in patches P and Q across the experiment, we investigated the relative contributions of potential source populations to the immigrant pool using genetic assignment tests. The assignment tests suggested some degree of immigration from all potential sources, consistent with the maintenance of genetic diversity in patches P and Q and with the increase in numbers of loci out of HWE (which may be ascribed to a Wahlund effect), over time. Therefore, immigration into the target patches appeared to follow largely a migrant pool model with the mixing of individuals from multiple potential sources, as opposed to a propagule pool model with migration from a single source [71].

On the one hand, minimal genetic differentiation between samples from consecutive years in patches P and Q (i.e. non-significant temporal F ST estimates) suggests a relatively stable pool of potential immigrants over time (electronic supplementary material, figure S1). Nonetheless, the genetic assignment of individuals captured in patches P and Q in 2005 and 2008 indicated potentially subtle variation in the source of migrants. Specifically, we observed a decrease in the relative contribution of inferred dispersers from patch M to both patches P and Q between those years, and an increase in the contribution from patch O. Patches R and S showed a relatively stable contribution of immigrants between 2005 and 2008, with their contribution into patch P and Q changing only moderately between these years. Patch M is the largest in the spatial network and, among the potential source patches for P and Q, contains the largest population of P. smintheus in most years [43]. Interestingly, despite the decline in patch M’s inferred contribution of immigrants to patches P and Q from 2005 to 2008, population size in M itself did not decline in that time period based on mark-recapture estimates [43,68]. In fact, in 2008 the largest population size estimate for M over the last two decades was observed, with the estimated local population size more than triple that estimated in 2005 [43]. Negative density-dependent dispersal has been observed in this species, with individuals more likely to emigrate from less dense populations and less likely to leave more dense populations [31]; this effect may explain decreased immigration from patch M in the year in which population size there was larger. Regardless, the assignment test results support the hypothesis that the persistence of populations in patches P and Q was not dependent on a single primary source population, but rather on the broader network with multiple sources contributing to the yearly pool of immigrants.

We recognize important caveats and assumptions in our assignment of individuals to potential source populations. First, we are assuming that all sampled individuals (i.e. all individuals removed from the experimental populations) are immigrants. Given the important role of immigration suggested by our analyses of genetic diversity and HWE over time, most individuals are indeed likely to be immigrants or direct descendants of immigrants; nonetheless, at least some adults removed from the experimental population may have been born locally. Second, the metapopulation is characterized by isolation-by-distance and adjacent populations are genetically similar [40,42]. This could limit our power to accurately assign individuals between adjacent potential source populations; we may have more confidence in assigning individuals generally to one side or the other of the targeted, experimental populations (i.e. M and O versus R and S) as those sets of populations are more strongly differentiated. Given these caveats, we do not interpret any individual assignment as unequivocally representing a specific immigration event. Instead, our aim was to assess general trends in immigration to the experimental patches over space and time, on aggregate and overall inferred assignments.

The increased representation of immigrants within P and Q in response to the experimental removals could potentially have been associated with selection for traits that facilitate immigration or colonization. For example, in a metapopulation of the Glanville fritillary butterfly (Melitaea cinxia), more isolated local populations show significantly higher frequencies of alleles at the Pgi locus associated with flight metabolism and ability compared to less isolated populations, suggesting selection for these alleles during dispersal or colonization [72,73]. We did not detect significant directional changes in the frequencies of genotypes at any individual loci in response to the experimental removals in patches P and Q, and therefore no evidence of selection in response to either the removal itself or subsequent immigration; however, this result is not surprising considering the moderate size of our SNP panel.

We present evidence for demographic rescue facilitated by immigration in response to a long-term, multi-generation artificial removal experiment. Interestingly, although loss of immigrants from the target patches P and Q to surrounding populations because of the experimental removals had little impact on the persistence of those surrounding populations [30,33], our results here indicate that immigration into patches P and Q, by contrast, was probably a key factor buffering those experimental populations against extinction. The lack of an effect of removals from P and Q on the persistence of surrounding populations may partially reflect ongoing movement and connectivity among the surrounding populations. Our results provide further evidence for the resilience of this metapopulation to both local and network-wide demographic decline [41,68] and highlight the importance of movement and connectivity for the persistence of organisms existing in fragmented landscapes [74].

Acknowledgements

We respectfully acknowledge that our field sampling occurred in the Treaty seven region of southern Alberta, which encompasses traditional territories of the Niitsitapi (Blackfoot) Confederacy (Siksika, Kainai and Piikani First Nations), the Tsuut’ina First Nation, the Îyâxe Nakoda Nations (Chiniki, Wesley and Bearspaw), and the Métis Nation (Region III). We thank Royal Society Open Science for the review of the manuscript and the staff at the Genome Québec sequencing centre for technical support. We also thank the many field assistants who conducted the experimental removals and mark-recapture studies in the Jumpingpound metapopulation, as well as the University of Calgary Biogeoscience Institute (Barrier Lake Field Station) for housing the crews during the field season.

Ethics

This work did not require ethical approval from a human subject or animal welfare committee.

Data accessibility

A copy of the SNP dataset is available on Dryad [75].

Supplementary material is available online [76].

Declaration of AI use

We have not used AI-assisted technologies in creating this article.

Authors’ contributions

K.Y.P.: conceptualization, data curation, formal analysis, funding acquisition, investigation, visualization, writing—original draft; M.L.: data curation, methodology, validation, writing—original draft; A.C.: data curation, methodology, validation; S.F.M.: conceptualization, methodology, resources, writing—review and editing; J.R.: conceptualization, methodology, writing—review and editing; N.K.: conceptualization, funding acquisition, methodology, project administration, supervision, writing—review and editing.

All authors gave final approval for publication and agreed to be held accountable for the work performed therein.

Conflict of interest declaration

We declare we have no competing interests.

Funding

Funding was provided by the Alberta Conservation Association grants in Biodiversity (K.Y.P.), NSF grants DEB-0326957 and 0918929 (S.F.M) and NSERC Discovery grants (N.K. and J.R.).
==== Refs
References

1. Thomas CD , Singer MC , Boughton DA . 1996 Catastrophic extinction of population sources in a butterfly metapopulation. Am. Nat. 148 , 957–975. (10.1086/285966)
2. Dirzo R , Young HS , Galetti M , Ceballos G , Isaac NJB , Collen B . 2014 Defaunation in the Anthropocene. Science 345 , 401–406. (10.1126/science.1251817)25061202
3. Allee WC . 1931 Animal aggregations: a study in general sociology. Chicago, IL: University of Chicago Press. (10.5962/bhl.title.7313)
4. Angulo E , Luque GM , Gregory SD , Wenzel JW , Bessa-Gomes C , Berec L , Courchamp F . 2018 Review: Allee effects in social species. J. Anim. Ecol. 87 , 47–58. (10.1111/1365-2656.12759)28940239
5. Frankham R . 2008 Genetic adaptation to captivity in species conservation programs. Mol. Ecol. 17 , 325–333. (10.1111/j.1365-294X.2007.03399.x)18173504
6. Charlesworth D , Willis JH . 2009 The genetics of inbreeding depression. Nat. Rev. Genet. 10 , 783–796. (10.1038/nrg2664)19834483
7. Hufbauer RA , Szűcs M , Kasyon E , Youngberg C , Koontz MJ , Richards C , Tuff T , Melbourne BA . 2015 Three types of rescue can avert extinction in a changing environment. Proc. Natl Acad. Sci. USA 112 , 10557–10562. (10.1073/pnas.1504732112)26240320
8. Whiteley AR , Fitzpatrick SW , Funk WC , Tallmon DA . 2015 Genetic rescue to the rescue. Trends Ecol. Evol. 30 , 42–49. (10.1016/j.tree.2014.10.009)25435267
9. Gomulkiewicz R , Holt RD . 1995 When does evolution by natural selection prevent extinction? Evolution 49 , 201–207. (10.1111/j.1558-5646.1995.tb05971.x)28593677
10. Gonzalez A , Ronce O , Ferriere R , Hochberg ME . 2013 Evolutionary rescue: an emerging focus at the intersection between ecology and evolution. Phil. Trans. R. Soc. B 368 , 20120404. (10.1098/rstb.2012.0404)23209175
11. Carlson SM , Cunningham CJ , Westley PAH . 2014 Evolutionary rescue in a changing world. Trends Ecol. Evol. 29 , 521–530. (10.1016/j.tree.2014.06.005)25038023
12. Samani P , Bell G . 2010 Adaptation of experimental yeast populations to stressful conditions in relation to population size. J. Evol. Biol. 23 , 791–796. (10.1111/j.1420-9101.2010.01945.x)20149025
13. Cohen ZP , François O , Schoville SD . 2022 Museum genomics of an agricultural super-pest, the Colorado potato beetle, Leptinotarsa decemlineata (Chrysomelidae), provides evidence of adaptation from standing variation. Integr. Comp. Biol. 62 , 1827–1837. (10.1093/icb/icac137)36036479
14. Vander Wal E , Garant D , Festa-Bianchet M , Pelletier F . 2013 Evolutionary rescue in vertebrates: evidence, applications and uncertainty. Phil. Trans. R. Soc. B 368 , 20120090. (10.1098/rstb.2012.0090)23209171
15. Brown JH , Kodric-Brown A . 1977 Turnover rates in insular biogeography: effect of immigration on extinction. Ecology 58 , 445–449. (10.2307/1935620)
16. Kanarek AR , Webb CT , Barfield M , Holt RD . 2015 Overcoming Allee effects through evolutionary, genetic, and demographic rescue. J. Biol. Dyn. 9 , 15–33. (10.1080/17513758.2014.978399)25421449
17. Wright S . 1931 Evolution in Mendelian populations. Genetics 16 , 97–159. (10.1093/genetics/16.2.97)17246615
18. Frankham R . 2015 Genetic rescue of small inbred populations: meta-analysis reveals large and consistent benefits of gene flow. Mol. Ecol. 24 , 2610–2618. (10.1111/mec.13139)25740414
19. Hedrick PW , Garcia-Dorado A . 2016 Understanding inbreeding depression, purging, and genetic rescue. Trends Ecol. Evol. 31 , 940–952. (10.1016/j.tree.2016.09.005)27743611
20. Kuo CH , Janzen FJ . 2004 Genetic effects of a persistent bottleneck on a natural population of ornate box turtles (Terrapene ornata). Conserv. Genet. 5 , 425–437. (10.1023/B:COGE.0000041020.54140.45)
21. Ortego J , Aparicio JM , Calabuig G , Cordero PJ . 2007 Increase of heterozygosity in a growing population of lesser kestrels. Biol. Lett. 3 , 585–588. (10.1098/rsbl.2007.0268)17609170
22. Vergeer P , Sonderen E , Ouborg NJ . 2004 Introduction strategies put to the test: local adaptation versus heterosis. Conserv. Biol. 18 , 812–821. (10.1111/j.1523-1739.2004.00562.x)
23. Bijlsma R , Westerhof MDD , Roekx LP , Pen I . 2010 Dynamics of genetic rescue in inbred Drosophila melanogaster populations. Conserv. Genet. 11 , 449–462. (10.1007/s10592-010-0058-z)
24. Vilà C et al . 2003 Rescue of a severely bottlenecked wolf (Canis lupus) population by a single immigrant. Proc. R. Soc. B 270 , 91–97. (10.1098/rspb.2002.2184)
25. Fitzpatrick SW , Bradburd GS , Kremer CT , Salerno PE , Angeloni LM , Funk WC . 2020 Genomic and fitness consequences of genetic rescue in wild populations. Curr. Biol. 30 , 517–522.e5.(10.1016/j.cub.2019.11.062)31902732
26. Hanski I . 1998 Metapopulation dynamics. Nature 396 , 41–49. (10.1038/23876)
27. Tallmon D , Luikart G , Waples R . 2004 The alluring simplicity and complex reality of genetic rescue. Trends Ecol. Evol. 19 , 489–496. (10.1016/j.tree.2004.07.003)16701312
28. Millon A , Lambin X , Devillard S , Schaub M . 2019 Quantifying the contribution of immigration to population dynamics: a review of methods, evidence and perspectives in birds and mammals. Biol. Rev. 94 , 2049–2067. (10.1111/brv.12549)31385391
29. Virtanen EA , Söderholm M , Moilanen A . 2022 How threats inform conservation planning—a systematic review protocol. PLoS One 17 , e0269107. (10.1371/journal.pone.0269107)35639722
30. Matter SF , Roland J . 2010 Effects of experimental population extinction for the spatial population dynamics of the butterfly Parnassius smintheus . Oikos 119 , 1961–1969. (10.1111/j.1600-0706.2010.18666.x)
31. Roland J , Keyghobadi N , Fownes S . 2000 Alpine Parnassius butterfly dispersal: effects of landscape and population size. Ecology 81 , 1642–1653. (10.1890/0012-9658(2000)081[1642:APBDEO]2.0.CO;2)
32. Matter SF , Roland J . 2002 An experimental examination of the effect of habitat quality on the dispersal and local abundance of Parnassius smintheus. Ecol. Entomol. 27 , 308–316. (10.1046/j.1365-2311.2002.00407.x)
33. Matter SF , Roland J . 2010 Local extinction synchronizes population dynamics in spatial networks. Proc. R. Soc. B 277 , 729–737. (10.1098/rspb.2009.1520)
34. Matter SF , Roland J . 2013 Mating failure of female Parnassius smintheus butterflies: a component but not a demographic Allee effect. Entomol. Exp. et Appl. 146 , 93–102. (10.1111/j.1570-7458.2012.01279.x)
35. Lande R . 1988 Genetics and demography in biological conservation. Science 241 , 1455–1460. (10.1126/science.3420403)3420403
36. Nei M , Maruyama T , Chakraborty R . 1975 The bottleneck effect and genetic variability in populations. Evolution 29 , 1–10. (10.1111/j.1558-5646.1975.tb00807.x)28563291
37. Hedrick PW . 2010 Genetics of populations, p. 376, 4th edn. Sudbury, MA: Jones and Bartlett Publishers.
38. Roslin T , Syrjälä H , Roland J , Harrison PJ , Fownes S , Matter SF . 2008 Caterpillars on the run - induced defences create spatial patterns in host plant damage. Ecography 31 , 335–347. (10.1111/j.0906-7590.2008.05365.x)
39. Sperling FAH , Kondla NG . 1991 Alberta swallowtails and Parnassians: natural history, keys, and distribution. Blue Jay 49 , 183–191. (10.29173/bluejay5078)
40. Keyghobadi N , Roland J , Strobeck C . 1999 Influence of landscape on the population genetic structure of the alpine butterfly Parnassius smintheus (Papilionidae). Mol. Ecol. 8 , 1481–1495. (10.1046/j.1365-294x.1999.00726.x)10564454
41. Caplins SA , Gilbert KJ , Ciotir C , Roland J , Matter SF , Keyghobadi N . 2014 Landscape structure and the genetic effects of a population collapse. Proc. R. Soc. B 281 , 20141798. (10.1098/rspb.2014.1798)
42. Jangjoo M , Matter SF , Roland J , Keyghobadi N . 2020 Demographic fluctuations lead to rapid and cyclic shifts in genetic structure among populations of an Alpine butterfly, Parnassius smintheus. J. Evol. Biol. 33 , 668–681. (10.1111/jeb.13603)32052525
43. Matter SF , Keyghobadi N , Roland J . 2014 Ten years of abundance data within a spatial population network of the Alpine butterfly, Parnassius smintheus. Ecology 95 , 2985–2985. (10.1890/14-1054.1)
44. Lucas M . 2022 Effects of spatial and temporal heterogeneity on the genetic diversity of the alpine butterfly Parnassius smintheus. PhD thesis, Western University, London, Ontario, Canada.
45. Keyghobadi N , Roland J , Strobeck C . 2005 Genetic differentiation and gene flow among populations of the alpine butterfly Parnassius smintheus, vary with landscape connectivity. Mol. Ecol. 14 , 1897–1909. (10.1111/j.1365-294X.2005.02563.x)15910314
46. Catchen J , Hohenlohe PA , Bassham S , Amores A , Cresko WA . 2013 Stacks: an analysis tool set for population genomics. Mol. Ecol. 22 , 3124–3140. (10.1111/mec.12354)23701397
47. Jangjoo M . 2018 Spatial and temporal patterns of neutral and adaptive genetic variation in the alpine butterfly, Parnassius smintheus. PhD thesis, Western University, London, Ontario, Canada.
48. Boratyn GM , Thierry-Mieg J , Thierry-Mieg D , Busby B , Madden TL . 2019 Magic-BLAST, an accurate RNA-Seq aligner for long and short reads. BMC Bioinformatics 20 , 405. (10.1186/s12859-019-2996-x)31345161
49. Bryant DM et al . 2017 A tissue-mapped axolotl de novo transcriptome enables identification of limb regeneration factors. Cell Rep. 18 , 762–776. (10.1016/j.celrep.2016.12.063)28099853
50. Allio R , Scornavacca C , Nabholz B , Clamens AL , Sperling FA , Condamine FL . 2020 Whole genome shotgun phylogenomics resolves the pattern and timing of swallowtail butterfly evolution. Syst. Biol. 69 , 38–60. (10.1093/sysbio/syz030)31062850
51. Mitikka V , Hanski I . 2010 Pgi genotype influences flight metabolism at the expanding range margin of the European map butterfly. Ann. Zool. Fenn. 47 , 1–14. (10.5735/086.047.0101)
52. Niitepõld K , Smith AD , Osborne JL , Reynolds DR , Carreck NL , Martin AP , Marden JH , Ovaskainen O , Hanski I . 2009 Flight metabolic rate and PGI genotype influence butterfly dispersal rate in the field. Ecology 90 , 2223–2232. (10.1890/08-1498.1)19739384
53. R Core Team . 2021 The R project for statistical computing . Vienna, Austria: R Foundation for Statistical Computing. See https://www.r-project.org.
54. Adamack AT , Gruber B . 2014 Popgenreport: simplifying basic population genetic analyses in R. Methods Ecol. Evol. 5 , 384–387. (10.1111/2041-210X.12158)
55. Jombart T . 2008 Adegenet: a R package for the multivariate analysis of genetic markers. Bioinformatics 24 , 1403–1405. (10.1093/bioinformatics/btn129)18397895
56. Pinheiro JC . 2000 Linear mixed-effects models: basic concepts and examples. In Mixed-effects models in S and S-PLUS (ed. DM Bates ), pp. 3–56. New York: Springer. (10.1007/b98882)
57. Paradis E . 2010 Pegas: an R package for population genetics with an integrated–modular approach. Bioinformatics 26 , 419–420. (10.1093/bioinformatics/btp696)20080509
58. Rousset F . 2008 Genepop’007: a complete re-implementation of the Genepop software for windows and Linux. Mol. Ecol. Resour. 8 , 103–106. (10.1111/j.1471-8286.2007.01931.x)21585727
59. Goudet J. 2005 Hierfstat, a package for R to compute and test hierarchical F-statistics. Mol. Ecol. Notes 5 , 184–186. (10.1111/j.1471-8286.2004.00828.x)
60. Elff M . 2022 Mclogit: multinomial logit models, with or without random effects or overdispersion. R package version 0.9.6. See https://CRAN.R-project.org/package=mclogit.
61. Benjamini Y , Hochberg Y . 1995 Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Stat. Soc. B Met. 57 , 289–300. (10.1111/j.2517-6161.1995.tb02031.x)
62. Piry S , Alapetite A , Cornuet JM , Paetkau D , Baudouin L , Estoup A . 2004 GENECLASS2: a software for genetic assignment and first-generation migrant detection. J. Hered. 95 , 536–539. (10.1093/jhered/esh074)15475402
63. Rannala B , Mountain JL . 1997 Detecting immigration by using multilocus genotypes. Proc. Natl Acad. Sci. USA 94 , 9197–9201. (10.1073/pnas.94.17.9197)9256459
64. Garza JC , Williamson EG . 2001 Detection of reduction in population size using data from microsatellite loci. Mol. Ecol. 10 , 305–318. (10.1046/j.1365-294x.2001.01190.x)11298947
65. Kirkpatrick M , Jarne P . 2000 The effects of a bottleneck on inbreeding depression and the genetic load. Am. Nat. 155 , 154–167. (10.1086/303312)10686158
66. Schaper E , Eriksson A , Rafajlovic M , Sagitov S , Mehlig B . 2012 Linkage disequilibrium under recurrent bottlenecks. Genetics 190 , 217–229. (10.1534/genetics.111.134437)22048021
67. McEachern MB , Van Vuren DH , Floyd CH , May B , Eadie JM . 2011 Bottlenecks and rescue effects in a fluctuating population of golden-mantled ground squirrels (Spermophilus lateralis) . Conserv. Genet. 12 , 285–296. (10.1007/s10592-010-0139-z)
68. Jangjoo M , Matter SF , Roland J , Keyghobadi N . 2016 Connectivity rescues genetic diversity after a demographic bottleneck in a butterfly population network. Proc. Natl Acad. Sci. USA 113 , 10914–10919. (10.1073/pnas.1600865113)27621433
69. Coulson JC , Duncan N , Thomas C . 1982 Changes in the breeding biology of the herring gull (Larus argentatus) induced by reduction in the size and density of the colony. J. Anim. Ecol. 51 , 739–756. (10.2307/4002)
70. Bosch M , Pocino N , Carrera E . 2019 Effects of age and culling on movements and dispersal rates of yellow-legged gulls (Larus michahellis) from a western mediterranean colony. Waterbirds 42 , 179. (10.1675/063.042.0204)
71. Wade MJ . 1978 A critical review of the models of group selection. Q. Rev. Biol. 53 , 101–114. (10.1086/410450)
72. Haag CR , Saastamoinen M , Marden JH , Hanski I . 2005 A candidate locus for variation in dispersal rate in a butterfly metapopulation. Proc. R. Soc. B 272 , 2449–2456. (10.1098/rspb.2005.3235)
73. Hanski I , Mononen T . 2011 Eco-evolutionary dynamics of dispersal in spatially heterogeneous environments. Ecol. Lett. 14 , 1025–1034. (10.1111/j.1461-0248.2011.01671.x)21794053
74. Wang S , Haegeman B , Loreau M . 2015 Dispersal and metapopulation stability. PeerJ 3 , e1295. (10.7717/peerj.1295)26557427
75. Park KY , Lucas M , Chaulk A , Matter SF , Roland J , Keyghobadi N . 2024 Data from: Immigration allows population persistence and maintains genetic diversity despite an attempted experimental extinction. Dryad. (10.5061/dryad.8cz8w9gzz)
76. Park KY , Lucas M , Chaulk A , Matter SF , Roland J , Keyghobadi N . 2024 Data from: Immigration allows population persistence and maintains genetic diversity despite an attempted experimental extinction. Figshare. (10.6084/m9.figshare.c.7370663)
