An evolutionary perspective on genetic load in small, isolated populations as informed by whole genome resequencing and forward-time simulations

Samples, sequencing, sequence processing, and genotype calling

Arizona samples were taken from donated hunter-harvested wings (Randel et al., 2019). Lysed blood tissues were collected from West Texas in 2018–2020 and from Central Texas in 2020 under a Texas Parks and Wildlife Department Scientific Permit SPR-0410-139 for research. Lysed blood from Mexico samples was collected under scientific collection permit number SGPA/ DGVS/ 06160/19 and were imported to the United States under USDA APHIS permit 140134 issued on December 08, 2019. Note that many of the details associated with our methods are in the extensive supplement and here we only provide a brief overview of our approach.

Genomic DNA was extracted from 98 Montezuma Quail samples from Arizona (AZ; N = 60), West Texas (WTX; N = 31), Central Texas (CTX; N = 3), and Mexico (MX; N = 4) using a modified potassium acetate protocol (Nicholls, Double, Rowell, & Magrath, 2000). Paired end libraries (2 × 150 bp) were created using a hyper Library construction kit from Kapa biosystems (Roche) and sequenced on S4 lanes using Illumina NovaSeq 6000. Raw sequences were trimmed to remove adapter sequences and low-quality bases, then mapped to the chromosome-level chicken genome (Gallus gallus assembly GRCg6a, genome size = 1.065 Gb). We used GATK v4.1.6 (McKenna et al., 2010) to discover Single Nucleotide Polymorphisms (SNPs) and call genotypes. After quality filtering, we analyzed genotypes at 9,324,827 SNPs in 66 individuals (AZ = 28, WTX = 31, CTX = 3, MX = 4) that were mapped on chicken autosomes 1–33 with 11.2× ± 2.0× (mean ± SD) depth of coverage (Data S1). Variant statistics, nucleotide diversity (π), pairwise relatedness, and heterozygosity (measured as the proportion of heterozygous nucleotides per individual genome) were calculated using VCFtools (Danecek et al., 2011).

Identification of runs of homozygosity

We used PLINK v.1.9 (Chang et al., 2015) to identify runs of homozygosities (ROHs) using multiple sets of parameter values (“Supplementary Methods” section). To identify the best set of parameters for our dataset, we performed a sensitivity analysis using package sensitivity (Pujol et al., 2015) in R ((R Core Team 2013); see “Supplementary Methods” section). We chose the final parameters that minimized the individual variation in F_ROH estimates between runs (Supplementary Methods). For the final ROH identification, we used the parameters --homozyg-window-het 2, --homozyg-snp 50, --homozyg-kb 100, --homozyg-window-snp 20 --homozyg-het 2 --homozyg-density 50 --homozyg-gap 1000 --homozyg-window-missing 2 --homozyg-window-threshold 0.05. We used detectRUNS package (Biscarini et al., 2018) in R to visualize PLINK output files and identify ROHs within each individual. We categorized ROHs into different classes based on length: short ROHs (200–300 kb), medium ROHs (300–600 kb), and long ROHs (>600 kb) to delineate short-term and long-term demographic effects (Ceballos et al., 2018). In general, long ROHs represent recent inbreeding whereas shorter ROHs are indicative of ancient inbreeding that has been broken down by generations of mutation and recombination.

Functional variant annotation, fitness effect classification, and genetic load estimation

We used Ensembl variant effect predictor (VEP) v101.0 (McLaren et al., 2016) to identify synonymous and non-synonymous (missense + nonsense) mutations within the exome. We considered only non-synonymous SNPs as potentially functional because they are easily identifiable using bioinformatic approaches and tractable to characterize across the annotated genome. We acknowledge that coding sequences represents only <2% of the whole genome and our approach ignores the putative functional effects of regulatory or non-coding intergenic regions. However, unlike synonymous and non-coding substitutions, non-synonymous mutations change the associated amino acid sequences (and thus their associated physical properties, which influence protein structure and function) whose long-term functional impact could be quantified via evolutionary constraint across the phylogeny. We predicted deleteriousness of potentially functional mutations using SIFT (Sorting Intolerant from Tolerant) scores (Ng, 2003). SIFT scores were generated for all possible amino acid changes relative to the chicken reference genome (assembly GRCg6a) and then categorized based on their deleteriousness (i.e., “deleterious” or “tolerated”). Clearly, SIFT scores do not directly estimate the fitness effects across many thousands of mutations but they do provide a quantitative proxy for the genotype-fitness relationship. Deleterious mutations were non-synonymous mutations with SIFT score < 0.05, whereas tolerated mutations were non-synonymous mutations with SIFT score ∈ [0.05, 0.1]. Tolerated mutations are expected to impart nearly neutral effects relative to deleterious mutations and thus subject to more relaxed purifying selection. These thresholds are conservative estimates and prescribed as best-practices by the tool developers (Ng, 2003); furthermore, these same demarcations have been used in similar studies in other species (Robinson et al., 2019^; Tennessen et al., 2012) and thus allow comparisons among studies.

Populations or genomes with more deleterious alleles across loci have a higher genetic load than populations or genomes with fewer deleterious loci. One major issue with comparing genetic load based on the number of deleterious mutations is that it is biased against larger populations with more diverse genomes (Ralls, Sunnucks, Lacy, & Frankham, 2020). Individuals from more diverse populations are expected to carry not only a greater number of neutral mutations but also a greater number of deleterious mutations across their genomes (Kim et al., 2017). To estimate individual genetic load, we not only compared the number of homozygous and heterozygous deleterious SNPs within individual genomes from both populations but also compared individual potential load (Load_P) and realized load (Load_R). These two metrics measure the proportion of deleterious sites present among all non-synonymous SNPs within an individual genome (Load_P) and the proportion of deleterious alleles that are in homozygous state within individual genomes (Load_R), respectively (see Supplementary Note 1; Mathur & DeWoody, 2021).

Besides the number of deleterious mutations, populations with higher frequency of deleterious alleles would have higher load and lower fitness than populations with the same deleterious mutations but at lower allele frequencies (Lynch & Gabriel, 2017). R_XY is a population-specific statistical estimate of genetic load that measures the relative frequency of a mutation between a pair of populations (Do et al., 2015). We estimated R_XY from allele frequency data of deleterious and tolerated functional SNPs using custom scripts (see “Data Accessibility” section).

Haplotype phasing, shared ancestry, and genealogical history of Montezuma Quail populations

Avian species lack PRDM9-directed recombination (Baker et al., 2017^; Singhal et al., 2015) and thus avian genomes are relatively stable over evolutionary time (Singhal et al., 2015). Fixed recombination hotspots and the conservation of synteny between quail and chicken (Morris et al., 2020) allowed us to use the chicken recombination map (Groenen et al., 2009). We also assumed that sex-linked regions are syntenic between chicken and Montezuma quail genomes. We phased Montezuma Quail genotypes into haplotypes and imputed missing genotypes using recombination maps for chicken autosomes 1–28 (except for chromosome 16 due to lack of polymorphic markers; Groenen et al., 2009^; Morris et al., 2020) using Beagle v.5.1 (Browning et al., 2018). We ran Beagle for each chromosome separately with burnin = 200, iterations = 500, phase-states = 200, and ancestral Ne = 100,000 as our best demographic model indicated that Montezuma Quail had a large ancestral population (see “Results” section). To assess the accuracy of phasing, we estimated switch error rate (i.e., the fraction of genotypes with switched haplotypes at a SNP for each individual between a pair of runs) from 10 independent runs of the phasing algorithm using VCFtools. The mean switch error was estimated to be 7.4% (SD = 1.5%). We used fineSTRUCTURE v4.1.1 (Lawson et al., 2012) to estimate the co-ancestry among individuals using phased haplotypes (N = 132 haplotypes; “Supplementary Methods” section).

Reconstructing recent demographic history

We used Genetic Algorithm for Demographic Model Analysis (GADMA) (Noskova et al., 2020) to find the best demographic model for two populations: Arizona and Texas (“Supplementary Methods” section). To infer a demographic model that best fits the joint distributions of the observed site frequency spectra (SFS), we simulated SFS under multiple possible demographic models to identify one that could best explain the observed SFS. We only considered synonymous mutations (N = 210,505) for calculating observed SFS as they are expected to segregate purely due to demographic changes and not due to selective changes. Here, we combined CTX and WTX samples into a single population (TX) due to high shared co-ancestry between the two populations (see “Results” section). We used AZ (N = 28) and TX (N = 34) samples and projected the observed SFS in [25 × 25] grid size. To calculate confidence intervals (CI) for estimated parameters, we bootstrapped our SFS 100 times by random sampling of 25 individuals from each population. For each model, we started with a series of simple initial structure that includes the number of time intervals that occur before and after each single splitting event and ran forward-time simulations using moments within GADMA (“Supplementary Methods” section). We used a generation time of one year as Montezuma Quail are effectively an annual species and a point mutation rate of µ = 3.14 × 10^–9 bp^–1 year^–1 as previously estimated for Montezuma Quail (Mathur & DeWoody, 2021).

Estimating genealogical age of deleterious mutations

We leveraged the phased data from 132 haplotypes to estimate the age of deleterious alleles by inferring the time to the most recent common ancestor (TMRCA) of deleterious mutation. We used Genealogical Estimation of Variant Age (GEVA) (Albers & McVean, 2020) to reconstruct genealogical trees of variable sites along each chromosome and estimate the age of deleterious mutations. GEVA estimates age of a mutation by identifying identical by descent genomic regions that are shared between a pair of haplotypes. The allele is assumed to be derived from a mutation event in the genome of the common ancestor and shared by all descendent haplotypes. A probabilistic estimate of the time of origin of a mutation is obtained by combining the cumulative distributions for pairs of haplotypes that share the mutation (“concordant pairs”) and the pairs that do not (“discordant pairs”) (Albers & McVean, 2020). We used recombination rate estimates from the chicken (Groenen et al., 2009) and an estimated point mutation rate of 3.14 × 10^–9 (Mathur & DeWoody, 2021). We used the effective population sizes inferred from EM optimization by ChromoPainter as a scaling parameter and joint molecular clock to estimate the age of deleterious mutations in different Montezuma Quail populations (“Supplementary Methods” section). We ran GEVA analysis 10 independent times to estimate the mean age and 95% CI.

Simulating temporal dynamics of non-neutral mutations

To better understand how drift and selection shape the distribution of non-neutral mutations during recent changes in demographic history, we performed forward-time simulations using SLiM v.3.6 (Haller & Messer, 2019). We used SLiM to model the distribution of functional mutations that arise (via mutation) and segregate conditioned on Montezuma Quail demographic model of the AZ and TX populations (see “Supplementary Methods and Results” section). The genetic and demographic parameters of the evolutionary genomics model simulated here are illustrated in Figure S1. We simulated the exome of a Galliform-like bird that carries a total of 2,256 genes of equal length (1,500 bp) distributed proportionally to the number of genes found in each chicken chromosome. Recombination between genes occurred at a uniform rate of r = 1 × 10^–3 crossovers per bp per generation (Kyriazis et al., 2020). Mutations arose in our simulation at a rate of µ = 3.14 × 10^–8 corresponding to Montezuma Quail estimates (Mathur & DeWoody, 2021); both r and µ were scaled up by a factor of 10 to account for the 1/10th population size of the simulated model (see below). We simulated neutral (s = 0), deleterious (s < 0), and beneficial mutations (s > 0) with the ratio of neutral:deleterious:beneficial mutations fixed at 1000:2310:2.31 (Figure S1A; Kim et al., 2017). The associated selection coefficients for the deleterious mutations were drawn from a gamma distribution with mean = –0.03 and shape parameter (α) = 0.186 as estimated for humans (Kim et al., 2017). For beneficial mutations, we kept the same distribution and α parameter but changed the mean to +0.03. We did not remove beneficial or deleterious mutations once fixed in the population as they still contribute to individual fitness (i.e., fixed non-neutral mutations can still increase/decrease individual fitness as homozygotes). We modeled a complex but biologically relevant dominance scheme that determines the dominance coefficient (h) based on the selection coefficients (s). The most deleterious mutations (s <= –0.05) were considered completely recessive (h = 0) whereas the most beneficial mutations (s >= 0.05) were considered completely dominant (h = 1). Dominance coefficients for the remaining mutations (–0.05 < s < 0.05) was determined as a step function between h = 0 and h = 1 with h increasing by 0.1 for every 0.01 increase in s; truly neutral mutations were completely additive (h = 0.5; Figure S1B).

Since SLiM explicitly models each individual within the population, it is memory intensive and thus, we rescaled population sizes to improve simulation performance (Haller & Messer, 2019). We rescaled our effective populations size (N_e) to 1/10th of that inferred from GADMA, and increased population genetic estimates of mutation rate (µ), recombination rate (r), and selection coefficient (s) 10× so that the products N_eµ, N_er, and N_es remain unchanged (Figure S1C). Rescaling genetic parameters allowed us to simulate a feasible miniaturized demographic model, but of course this skews absolute values (e.g., artificially reducing the time required to fix beneficial alleles) as compared to true population sizes (see SLiM documentation for details; https://github.com/MesserLab/SLiM/releases/download/v4.0.1/SLiM_Manual.pdf). However, since rescaling should affect both populations equally, our relative comparisons of genetic metrics between the two populations should still be meaningful. At generation 1, a single ancestral population (p₀) was created with size = 1,1750 that ran for 1,17,500 generations at a constant population size to reach equilibrium (“burn-in”). At the end of the 117,500th generation, p₀ expanded instantaneously to 1,20,000 individuals and evolves at constant size for another 7,500 generations. At the end of 1,25,000^th generation, p₀ splits in an 80:20 ratio to create two new populations p₁ and p₂. In our scenario, p₁ represents the larger AZ population (“AZ-like”) and p₂ represents the smaller TX population (“TX-like”). Immediately after the split, TX-like population size changes to = 2,750 and remains constant until the last 200 years when it begins to exponential decline (λ_decay = 1.10 × 10^–3) until present day ( = 2211). The AZ-like population starts at = 96,000 post-split and exponentially declines till present day to = 3,600 (λ_decay = 2.19 × 10^–3). The migration rates were set to m₁₂ = 4 × 10^–3 and m₂₁ = 4 × 10^–4 where m₁₂ represents migration from p₁ to p₂ and vice-versa. To simulate recent complete lack of gene flow between the two populations, we set both migration rates to 0 in the last 200 years (Figure S1C). During the simulations, after the founding split, we kept track of genetic diversity indices (heterozygosity, F_ROH), relative fitness, and average number of different impact type mutations from random sampling of 100 individuals from each population (p₁, p₂) every 50 generations following Kyriazis et al., (2020). Relative fitness is measured as , i.e., mean absolute fitness of 100 sampled individuals within a population ( ) at a given time (T = t) as compared to mean absolute fitness immediately before the split (T = t₀). Thus, relative fitness > 0 represents that individuals in the population are relatively more fit than they were at the time of the founding bottleneck. The entire simulation post-split was run 100 times to capture the variation between individual runs.

Statistical analysis

All pairwise comparisons of means were performed using Wilcoxon rank sum test using false discovery rate method of p-value correction in R unless mentioned otherwise. Statistical details are provided in the Supplementary Tables.

Contemporary genomic variation in Montezuma Quail populations

Individuals from the larger AZ population had significantly higher genome-wide estimates of heterozygosity and lower inbreeding as compared to quail from the smaller and more isolated WTX population (Tables S1–S3 and Figure S2). As expected, samples from the smallest and most isolated Montezuma Quail population in Central Texas (CTX) were the least diverse and most inbred samples among all Montezuma Quail populations (Figure S2). Heterozygosity in the MX samples was intermediate between AZ and WTX, and genome-wide inbreeding was intermediate between CTX and WTX (Tables S1–S3 and Figure S2). We found a significant pattern of reduced genome-wide heterozygosity in all TX samples (Figure S3) which is also indicated by more of the TX quail genome covered by short ROHs (200–300 kb) relative to AZ quail genomes (Figure 1B). All Montezuma Quail individuals carry long ROHs (>600 kb) but medium ROHs (300–600 kb) are more abundant in Texas quail (Figure 1B).

Shared ancestry and demographic history of Montezuma Quail populations

Genealogical trees based on shared ancestry indicate that AZ and TX populations are each monophyletic clades and that the MX population is a sister clade to the AZ population (Figure 1C). The CTX population has an ancestry tied closely to that of a single WTX clade. All individuals within a population share more ancestry than among populations, meaning there is a lack of contemporary gene flow among populations (Figure 1D). We also identified a few familial groups that shared highest co-ancestry within WTX samples, but excluding those samples had no effect on our diversity estimates (see “Supplementary Methods” section).

We inferred the demographic trajectory of the two divergent Montezuma Quail lineages by simulating the site frequency spectrum (SFS) of a 2-population demographic model with migration and dynamic population sizes. The most probable demographic model over the last ~100,000 years (Figure 2A and Figure S4) indicates that the ancestral Montezuma Quail population expanded ~90.5 kya (95% CI = 81,835–1,03,372) and was relatively large and stable for ~75,000 years until it split around 16.6 kya (14,515–28,951; Figure 2A and Table S4), ultimately producing the modern AZ and TX populations. The model estimates that Arizona and Texas populations split in a ~80:20 ratio, meaning that modern Texas populations were formed from fewer founders than the Arizona population (i.e., Texas birds went through a stronger bottleneck, ~4-fold more intense than Arizona; Figure 2A and Table S4). Since the split, Texas populations have consistently maintained smaller sizes (N_TX = 27,711 [24,391–33,532]), whereas the modern Arizona population has resulted from an exponentially decline to present size (N_AZ = 35,656 [33,510–44,141]). These results show the different population dynamics of AZ and TX populations where Texas populations have remained small for ~16,000 years after a strong initial bottleneck, whereas the modern AZ population (which is currently much larger than TX) has been progressively declining in size over the same timeframe (possibly via successive founder events associated with range expansion). Our model also recovered an unequal rate of historic gene flow between the two populations, indicating a ~10× higher historical migration rate from AZ to TX relative to gene flow from TX to AZ (Figure 2A and Table S4).

Open multimedia modal

Contemporary genetic load in Montezuma Quail populations

Genetic load of a population is largely dependent on the number of deleterious mutations across the genome and the frequencies of deleterious alleles at each locus within the population (Bertorelle et al., 2022). The contemporary AZ population not only contains more functional non-synonymous SNPs (deleterious, neutral, and putatively beneficial mutations) across loci but also more deleterious mutations across loci than in the smaller contemporary WTX population (Figure 3A and Table S5). This is largely consistent with recent evidence that larger populations not only have higher genetic diversity but also have more loci harboring deleterious alleles than smaller populations. However, the relative frequency of deleterious mutations (R_XY) is higher in the smaller WTX population than larger AZ population (R_XY = 1.10 ± 0.002; Figure 3B) which means deleterious mutations are comparatively rarer in the larger AZ population. However, all functional mutations mostly exist at similar frequencies in both populations (R_XY-Tolerated = 1.08 ± 0.002; R_{XY-Nonsynonymous} = 1.02 ± 0.001). Deleterious mutations become common in small populations as weak purifying selection plus pronounced drift increase their frequencies over time, compared to larger populations with more efficient selection and less drift where deleterious mutations are maintained at lower frequencies. Thus, even though the smaller WTX population has fewer deleterious loci than larger AZ population, alleles at those loci are segregating at higher frequencies in the smaller WTX population.

Open multimedia modal

From an individual perspective, WTX quail carry significantly more deleterious alleles masked in heterozygotes than do AZ quail, but there is no difference between tolerated alleles in heterozygote (Figure 3C). Since AZ individuals are more genetically diverse, they have significantly more overall functional heterozygosity. WTX individuals have more deleterious and tolerated alleles in homozygous states than do more outbred AZ individuals (Figure 3D). To compare if individuals in the larger AZ population carry a disproportionally higher number of functional mutations that are deleterious, we compared their individual potential load (Load_P). We found that individuals in both AZ and WTX populations have comparable Load_P but individuals in the larger AZ population carry a higher proportion of tolerated functional SNPs than individuals from the smaller WTX population (Figure S5 and Tables S6 and S7). Since WTX individuals are more inbred, homozygosity is higher and this is reflected in Load_R for both deleterious and tolerated SNPs (Figure S5). When we compared individual Load_R of homozygous deleterious mutations present in ROHs, we found that long ROHs carry significantly more deleterious mutations than tolerated mutations (Figure S6). Moreover, medium and short ROHs are more abundant in tolerated mutations than deleterious mutations (Figure S6). These results indicate that recent inbreeding is more likely to increase the genetic load of deleterious mutations as they accumulate disproportionately in longer ROHs. Details about genetic load in CTX and MX samples are provided in Figure S5 and Table S6.

Demographic impact on distribution of deleterious mutations in Montezuma Quail populations

We estimated the age of deleterious mutations that are currently present in AZ and WTX populations to understand how drift and selection shaped the number of deleterious alleles and their frequencies considering demographic histories of different populations. Our results indicate that almost all the deleterious mutations in Montezuma Quail populations originated within the past 150 kya (Figure S7). Within the last 100 kya, more deleterious mutations arose in AZ (N = 9400 ± 35) than in the WTX population (N = 8590 ± 23; Figure 4A). In both populations, more than half of the deleterious mutations currently segregating are ancestral; meaning they arose ~25 kya well before the population split. When categorizing into age bins, we found that both populations carry fewer deleterious mutations that arose around the time of bottleneck (15–25 kya) than at other time epochs, and the smaller WTX population has fewer mutations that arose post-bottleneck than the larger AZ population (Figure 4A). Next, we looked at how many of the deleterious mutations within each population are shared with the other population versus the number of private mutations exclusively within each population. Most of the deleterious mutations that are currently segregating in each population are young and private in origin, whereas most of the deleterious mutations that are shared between the two populations arose in the ancestral population (Figure S7). Most of the deleterious SNPs are shared but AZ has ~25% more private SNPs than WTX (Figure 4B). However, with more inbreeding in Texas quail (Figure 1B), WTX individuals have significantly more shared and private deleterious alleles as homozygotes (Figure 4C). The distribution of allele frequencies is similar in both populations with respect to deleterious alleles that arose in the ancestral population and are shared between the two populations (Figure 4D and Figure S8). However, younger private mutations are more abundant in the smaller WTX population (Figure S8) and private deleterious mutations reached higher frequency in the smaller WTX population over time (Figure 4E).

Open multimedia modal

Forward-time, evolutionary genomic simulations of populations with different demographies

We simulated the evolution of the exome of a Galliformes-like species with neutral, deleterious, and beneficial mutations that evolved under a demographic model inferred from real Montezuma Quail populations (Figure S1). After initial expansion, the ancestral population split into two unequal-sized populations: the smaller population (“TX-like”) remained small and stable over evolutionary time, whereas the larger population (“AZ-like”) went through a slow exponential decline in population sizes. Both populations had very low and unequal rates of gene flow. Within the last 200 years, the TX-like population again started to slowly decline, and the two populations became completely isolated.

Our simulations show that in the first 500 generations after the founding bottleneck, the small TX-like population that underwent the stronger bottleneck lost genetic diversity more rapidly with a concurrent increase in inbreeding than the AZ-like population (Figure 5A and B). Genetic diversity in the larger AZ-like population responded much more slowly immediately post-bottleneck. However, despite a stronger bottleneck, the TX-like population still maintained relative fitness (i.e., fitness of a population at a given time relative to fitness at the time of bottleneck) in the first 500 years (Figure 5C). Our simulations indicate that inbreeding effectively exposed deleterious mutations in homozygotes and, consequently, that purifying selection purged deleterious mutations more rapidly from the TX-like population (Figure 5D and E) while simultaneously maintaining strongly beneficial mutations (Figure 5F). Thus, our simulations suggest that small populations that survive severe bottlenecks may prevent immediate loss of fitness by purging deleterious alleles present in the ancestral population without necessarily losing ancestral adaptive alleles.

Open multimedia modal

Over the next 800 years (i.e., from T = 500 –1,300 years post-bottleneck), as the small TX-like population remained small, genetic drift was the dominant evolutionary force acting on the TX-like population that led to a slow loss of genetic diversity and increased inbreeding, which stabilized over time (Figure 5A and B). After the initial purging of strongly deleterious mutations, genetic drift continued to stochastically reduce both the number of weakly deleterious (Figure 5E) and weakly beneficial mutations over time in the small TX-like population (Figure 5G). Relative fitness of both populations also fluctuated more stochastically in this time (Figure 5C). For the larger but slowly shrinking AZ-like populations, diversity declines as inbreeding gradually increases and deleterious alleles are purged (Figure 5D and E), which leads to increases in relative fitness (Figure 5C). The AZ-like populations were less likely to lose weakly beneficial mutations (Figure 5G) because of the increased efficiency of selection when population sizes are larger. This means that as populations continue their independent demographic trajectory (with limited gene flow), larger populations are less impacted by drift and more impacted by selection. In the last 200 years (i.e., T = 1,300–1,500 years post-bottleneck), the TX-like populations were isolated and again started losing genetic diversity and became more inbred as the population declined (Figure 5A and B).

When we examined the allele frequencies of non-neutral mutations in the two populations over time, relative frequencies of deleterious mutations remain higher (R_XY > 1) but stable over time in the small, isolated TX-like population as compared to the larger AZ-like population (Figure S9A). However, frequencies of beneficial mutations increased more rapidly in the smaller TX-like population before stabilizing when beneficial alleles were fixed (Figure S10). Our simulation results suggest that purifying selection against rare deleterious mutations may remain effective in small populations (i.e., N_e may be small, but s is large so N_e(s) > 1) but higher drift and/or positive selection acting on beneficial alleles maintained in small populations post-bottleneck fixes these common beneficial mutations more rapidly in the smaller populations as compared to the larger population where time of fixation is relatively slower (Figure S10). This is also explained by the levels of differentiation between the two simulated populations captured at functional and neutral mutations over time (Figure S9B). F_ST increases as populations diverge over time, but the level of differentiation is much higher for adaptive alleles (i.e., beneficial mutations). In contrast, the differentiation at deleterious mutations was comparable to neutral mutations (Figure S9B). When we investigated the distribution and frequency of non-neutral mutations at the end of the simulation, we observed that there were more considerably more deleterious mutations (N = 12571 ± 138) than beneficial mutations (211 ± 6) but on average beneficial mutations are considerably older and mostly fixed between the two populations (Figure S10). Just as in the empirical dataset, most deleterious mutations that are segregating in the contemporary populations at the end of the simulation are young and private (Figure S11).

Montezuma Quail populations as incipient ring species

Our comparisons of genome-wide variation patterns, shared ancestry, and demographic models all suggest that Montezuma Quail likely evolved similar to a ring species (Figure 2B; (Irwin et al., 2001a). Ring species evolve as populations diverge from their common ancestor under a scenario where an unsuitable landscape creates divergent populations that could subsequently speciate over time despite historic gene flow (Irwin et al., 2001b). In the case of Montezuma Quail, our demographic model indicates that a relatively small eastern population in the Sierra Madre Oriental mountains was separated around 15–25 thousand years ago from the large ancestral population by the Central Mexico Plateau (Figures 1A and 2). Following this separation, eastern populations remained relatively small, stable, and isolated to ultimately form the contemporary Texas populations. This scenario is evidenced by the modern-day distribution of the species (Figure 1A) and its suitable habitat, the demographic model (Figure 2A), the shorter ROHs in modern Texas genomes (Figure 1B), and the reduction of genome-wide diversity in Texas populations (Figure S3). In species with limited dispersal capacity, like ground-dwelling montane birds, range expansions are characterized by loss of genetic diversity along the expansion axis due to recurrent bottleneck effects (Excoffier et al., 2009) as seen in Arizona. Arizona populations of Montezuma quail were likely founded by serialized bottlenecks as the larger ancestral population along the Sierra Madre Occidental mountains in the West of Mexico continued to colonize new habitat (Leopold & McCabe, 1957). Geographic conditions suitable for ring species formation with reproductive isolation are rare (Alcaide et al., 2014), but stable terminal populations with historic gene flow can exist if the populations are locally adapted (Martins et al., 2013). Natural history studies of Montezuma quail have identified distinct morphological features between individuals from West or East of the Central Plateau in Mexico (Leopold & McCabe, 1957) that may result from local adaptation associated with their specialist habitat of montane piñon-oak grasslands (Heffelfinger & Olding, 2000). We do not mean to suggest that there is reproductive isolation between contemporary Montezuma quail populations in Arizona and Texas, that they now represent distinct species, or that they will speciate any time soon. Instead, we suggest that these populations are on incipient evolutionary trajectories similar to those that ultimately give rise to ring species (e.g., Greenish Warblers; Irwin et al. 2001a). Future phylogeographic and landscape ecology studies in Montezuma Quail could help identify admixture zones (like populations in New Mexico) and help illustrate the population dynamics of a species that evolved peripherally around a habitat barrier (Irwin et al., 2001b^; Petren et al., 2005).

Genetic purging and persistence of small populations following strong bottlenecks

Montezuma quail populations in Texas were established ~16.5 kya as part of a substantial bottleneck (i.e., an 80% reduction in size from the ancestral population; Figure 2A). This pattern is similar to other long-term small populations (or species) like polar bears (Benazzo et al., 2017), Channel Island foxes (Robinson et al., 2018), Tibetean wolves (Robinson et al., 2019), and Western mountain gorillas (van der Valk et al., 2019). Western mountain gorillas and island foxes have maintained small population sizes and reduced genetic diversity since they split from Eastern gorillas (T = 20 kya) and mainland gray foxes (T = 9 kya), respectively (Hofman et al., 2015^; Xue et al., 2015). Our simulations show strong bottlenecks like those associated with founder effects or speciation events can purge ancestral deleterious mutations rapidly and that surviving small populations suffer from elevated inbreeding and lower genetic diversity relative to populations that experienced less pronounced bottlenecks (Figure 5A–E). These empirical examples (quail, fox, gorilla) and our simulations are from natural populations that underwent historic bottlenecks, but reintroduced populations provide another set of examples of small populations that underwent a strong, recent bottleneck (e.g., Alpine Ibex reintroductions (Biebach & Keller, 2009) or the introduced steelhead trout (Willoughby et al., 2018)). These small, reintroduced populations also have lower overall genomic diversity (Willoughby et al., 2018) and have purged deleterious mutations that still occur in larger populations (Grossen et al., 2020). Thus, strong bottlenecks reduce overall genomic diversity, including the reduction of deleterious mutations that comprise the ancestral genetic load (Figure 5D and E).

Small populations can be stable populations despite reduced overall diversity and increased load (Ralls et al., 2020^; Reed, 2010^; Robinson et al., 2018). By including beneficial mutations in our simulations, we show that the small populations exibit lower overall diversity (e.g., heterozygosity), but can maintain relative fitness by maintaining and fixing strongly beneficial alleles via positive selection (i.e., “genetic thrust”). In our simulations, small populations purged strongly deleterious, but maintained strongly beneficial mutations immediately following the bottleneck (i.e. within 500 generations; Figure 5D–G). Rapid adaptation despite loss of diversity due to strong bottlenecks has been documented in introduced populations of steelhead trout (Willoughby et al., 2018). Theoretical and empirical evidence suggest that populations can maintain adaptive genetic variance immediately post-bottleneck (Fitzpatrick et al., 2020^; Mularo et al., 2022), but adaptive variation is more stochastic long-term depending on the evolutionary history and on the quantitative traits underlying fitness (Kardos et al., 2021). Taken together, our simulations suggest that strong bottlenecks like founder events or (re)introductions increase individual load due to inbreeding but ultimately lead to increases in population fitness by purging the ancestral deleterious alleles from the population. However, as witnessed from both empirical data (Figures 3B and 4E) and simulations (Figure S9A), deleterious mutations that are not purged from the small population tend to rise more rapidly over time and segregate at higher frequencies in the smaller population as compared to the larger population. Since overall genetic load is proportional to the product of the number of deleterious mutations and their allele frequencies (Bertorelle et al., 2022), small populations may reach a fitness equilibrium (Figure 5C), and thus, bottlenecked populations can remain stable despite higher population load.

Historically large populations risk inbreeding depression due to recent declines, but historically small populations also risk loss of adaptive potential

Large populations harbor more deleterious mutations by masking them in heterozygotes and escape genetic purging by avoiding inbreeding (Bertorelle et al., 2022). For many species, recent population declines and small dispersal distances associated with fragmented habitat raise concerns about recent inbreeding. Consanguineous matings create longer ROHs (Ceballos et al., 2018) and since deleterious mutations disproportionally appear in longer ROHs (Figure S6; (Szpiech et al., 2013), individuals of species with large historic population sizes that experienced recent inbreeding have higher individual load and carry signatures of inbreeding depression (Abascal et al., 2016^; Benazzo et al., 2017^; Stoffel et al., 2021).

Newly established populations require time to generate new de novo deleterious mutations, but these segregating deleterious mutations can rise in frequency over evolutionary timescales in small populations as shown in Montezuma quail from WTX and in other species ((e.g., Eastern massasauga rattlesnakes (Ochoa & Gibbs, 2021) or Mountain gorrilas (van der Valk et al., 2019)). Even in reintroduced populations of Alpine Ibex or endangered populations like Indian tigers or Chinese crocodile lizards, the relative frequency of functional non-synonymous variation was higher in smaller populations (Grossen et al., 2020^; Khan et al., 2021^; Xie et al., 2022). Weaker selection (i.e., N_e(s) ~ 1), either due to moderate selective coefficients and/or small population sizes, leads to less effective purging, which results in more private deleterious alleles that become more frequent as small populations continue to grow more isolated. Besides gaining population load over time, our simulations show that, relative to large populations, small isolated populations lose weakly beneficial mutations (Figure 5G), as expected by high drift and small population dynamics (Bijlsma & Loeschcke, 2012^; Frankham, 2005). Long-term small populations may be demographically stable, but their reduced levels of standing variation puts them at risk of losing adaptive potential because adaptation is strongly correlated with levels of genetic variation and is better predicted by genetic variation than levels of inbreeding (Ørsted et al., 2019).

Genealogical dating, simulations, and their limitations

Our goal was to use distributions of mutational ages to assess allele frequency trends of pre- and post-bottleneck deleterious variants in the two populations; we were mostly interested in the comparisons between populations with differing demographies rather than the absolute age estimates of specific variants. Estimating geneaological ages of deleterious mutations is dependent on the accuracy of haplotype phasing and improves with the size of the sequencing datasets (Albers & McVean, 2020). A wider distribution of independently sampled genomes improves genealogical resolution, but the evaluations of multiple runs capture the variation in the mutation age estimates and suggest our overall age distributions of deleterious mutations are robust (Figure S12). We acknowledge that our evolutionary genomics simulations of a miniaturized demographic model (Figure S1) could exaggerate the effect of drift despite our scaling the genetic parameters to account for small sizes (Haller et al., 2019). Thus, the absolute number of SNPs and their allele frequency estimates might be inaccurate as opposed to estimates from the simulations of true population sizes, but our population genetic inferences should be well-supported as they are either proportional to the genome sizes (heterozygosity, F_ROH) or relative to another population (relative fitness, R_XY, F_ST). This is evidenced by the fact that the temporal distribution of deleterious mutations currently segregating in the final populations at the end of our simulation was similar to our observations from the empirical genomic data (Figure 4B, D, and E and Figure S11), so our evolutionary genomic model recapitulates the evolutionary history of Montezuma Quail populations. A general limitation of simulation studies of natural populations is the lack of true estimates for genetic parameters like the distribution of fitness effects (DFE) of non-neutral mutations or the dominance-selection relationships (Ralls et al., 2020), which can lead to underestimates of load in empirical studies (Kardos et al., 2021). Herein, we used the same DFE estimates as other simulation studies (Kyriazis et al., 2020^; Robinson et al., 2019^; Xie et al., 2022), but we incorporated a more realistic dominance-selection relationship and a slow bottleneck based on the inferred demographic history. Finally, we only simulated the exome (i.e., protein coding regions) that represents only 1–2% of the whole genome and ignores the putative functional effects of regulatory or non-coding intergenic regions. Simulations like ours have limitations, but they can also provide powerful insights as to the impact of population trajectories on genetic effects. Future simulation studies should consider epistatic interactions between different genes (Whitlock et al., 1995) and the role of dominance in maintaining additive genetic diversity (Mularo et al., 2022).

Evolutionary history, demography, genetic load, and conservation

Our study highlights that genetic load is complex and context-dependent. The evolutionary history (e.g., phylogeography and demography; Figures 1 and 2) of a population or species has a profound impact on its overall genomic variation, including the genetic load of deleterious recessive alleles. Recent inbreeding and short-term demographics shape individual realized load (Figure S5), but long-term population size is deterministic of population load and adaptive potential (Ørsted et al., 2019). In the future, assisted conservation strategies such as supplementation or translocation of individuals from donor populations should continue to consider ecological limitations such as carrying capacity but may also consider evolutionary attributes (e.g., demographic history and genetic load) of the donor population(s). Our empirical data from natural populations indicates that natural selection against ancestral deleterious mutations purges them in small populations, much like artificial selection can (López-Cortegano et al., 2021^; Templeton & Reed, 1984).

Small isolated populations have been traditionally conserved via assisted gene flow from genetically diverse populations (Frankham, 2015^; Ralls et al., 2018). Such a strategy may risk the inadvertent introduction of new or previously purged deleterious mutations into the small recipient population, which could further increase realized genetic load (Kyriazis et al., 2020^; Robinson et al., 2019). We think that future assessments of genetic load could reduce the likelihood of conservationists inadvertently introducing a Trojan horse full of deleterious alleles into a genetically fragile recipient population. These are important considerations for genetic rescue, but major knowledge gaps remain in our understanding of fitness in natural populations. Conservationists contemplating genetic rescue must try to maximize adaptive potential while minimizing reductions to fitness by introducing novel deleterious alleles. Accordingly, we think Load_P measures the trade-off between the number of deleterious alleles/loci and other functional alleles within an individual and, at least in principle, translocations using individuals with lower Load_P values could increase the number of putatively functional variants in the recipient population. The influx of functional but non-deleterious alleles could help replenish the lost adaptive variation in small populations and improve their adaptive potential. From a practical perspective, we think that our results suggest that in cases where genetic rescue is a plausible conservation scenario, biologists should consider the demographic history, inbreeding, and genetic load in donor populations as part of the conservation planning.

Acknowledgements

The authors thank J. Heffelfinger at the Arizona Department of Game and Fish for the hunter-harvested wings. The authors also thank Genaro Olmos, Erasmo Olmos, and Juan Felipe Montoya for their Mexican field work support; Alberto Macías Duarte, Angel B. Montoya, James Weaver, Ryan Schmidt, Joyce Moore, and Mike Miller for Texas field work support. The authors also thank Drs Ryan Luna and Louis Harveson from Sul Ross State University; as well as Kristyn Stewart and Dr Fidel Hernandez from Texas A&M University-Kingsville.