Purging of deleterious burden in the endangered Iberian lynx.

SignificanceThe dynamics of deleterious variation under contrasting demographic scenarios remain poorly understood in spite of their relevance in evolutionary and conservation terms. Here we apply a genomic approach to study differences in the burden of deleterious alleles between the endangered Iberian lynx (Lynx pardinus) and the widespread Eurasian lynx (Lynx lynx). Our analysis unveils a significantly lower deleterious burden in the former species that should be ascribed to genetic purging, that is, to the increased opportunities of selection against recessive homozygotes due to the inbreeding caused by its smaller population size, as illustrated by our analytical predictions. This research provides theoretical and empirical evidence on the evolutionary relevance of genetic purging under certain demographic conditions.

In the last few centuries, the ever-growing anthropogenic pressures (1) have significantly accelerated the extinction of populations and species across a wide spectrum of taxa (2–4). Small populations are particularly susceptible to environmental and demographic stochastic events (5), as well as to negative genetic processes such as genetic drift, inbreeding depression, the accumulation of deleterious variation, and the loss of adaptive potential, which jeopardize their persistence (6). Most worryingly, these genetic hazards may persist even when the original cause of population decline is removed (7), creating a positive feedback loop as the population spirals toward extinction (8).

Albeit highly relevant for conservation and evolutionary biology, the dynamics of genetic load remain poorly understood. Classically, this load has been defined in terms of consequences on population fitness, although it can also be defined at the individual level. It includes the reduction of expected fitness due to segregating deleterious alleles (i.e., the expressed load), the load from recessive deleterious effects hidden in heterozygosis (i.e., the inbreeding load) responsible for inbreeding depression, and the slow fitness decline due to the continuous fixation of deleterious mutations (9). Characterizing these components of the genetic load is not an easy task, since measuring fitness in wild populations is extremely difficult (10). However, nowadays, whole-genome sequencing offers the opportunity to gain insights into the load of wild populations using a number of genomic summaries, such as the average number of (segregating plus fixed) derived alleles per individual (10), hereafter the derived count. This is possible provided there exists a good reference annotation that allows arranging the identified mutations into rough deleteriousness categories according to the genomic feature they lay in, their predicted effect on protein function, and the degree of conservation of the site across taxa (11). For deleterious categories, the derived count in individuals provides an appropriate and robust measure of the deleterious burden that can be interpreted at the individual or population level. Since it takes into account the burden from alleles occurring in heterozygosis, this measure is related to the whole fitness load, including the hidden inbreeding load. Therefore, the more recessive the deleterious effects, the looser the relationship between this count and the fitness load expressed in the individuals becomes (see ref. 10). When applied to populations with a common origin but contrasting demographies, such a genomic empirical approach can improve our understanding of how genetic and demographic factors shape the deleterious burden of small, endangered populations (12), which, in turn, may lead to more effective management. In these populations, genetic purging, that is, the increased selective pressure against (partially) recessive deleterious alleles due to increased homozygosis, is expected to translate into a reduction in both the inbreeding load (13) and the derived count in deleterious categories with prevailing recessive fitness effects. At the same time, the reduced efficacy of purifying selection caused by drift is expected to increase the derived count of alleles with small deleterious effects (14, 15). Therefore, genomic analysis can help to understand the relative roles of the purging of deleterious alleles and the drift-driven relaxation of natural selection in the evolution of the deleterious burden after a reduction in population size, and provide insights into their fitness consequences. In recent years, several studies focused mainly on the genetic consequences of the Out-Of-Africa bottleneck have contributed to the ongoing debate on whether human populations differ in their genetic load (10, 16–20). Only a few studies, however, have seized upon the power of genomics to assess the deleterious burden in small, endangered wild populations of nonmodel organisms (21–23).

The Iberian lynx (Lynx pardinus) and the Eurasian lynx (Lynx lynx) are the two extant species of the Lynx genus in Eurasia. After their divergence, dated around 1.1 million years ago (24, 25), they followed quite parallel demographic trajectories; however, Iberian lynx population sizes remained consistently smaller than those of the Eurasian lynx (26). Nowadays, they show highly contrasting distribution areas and conservation statuses. The Iberian lynx, a habitat and prey specialist endemic to the Iberian Peninsula, went through several bottlenecks and has maintained small population sizes throughout its recent history. The extreme and well-documented decline during the second half of the 20th century (27–29) resulted in the extirpation of all populations except for two in southern Spain: 1) a central population in Sierra de Andújar in the Sierra Morena range (AND), which maintained a comparatively larger population size and remained well connected to other populations until the mid-20th century, and 2) a smaller, peripheral population in the protected area of Doñana (DON), which has been effectively isolated for around two centuries (30, 31). The species was classified as “critically endangered” in the 2002 and 2008 International Union for Conservation of Nature (IUCN) red lists (32) after its global population dwindled to less than 100 individuals. Nonetheless, active conservation efforts, which include in situ, ex situ, and reintroduction programs (33, 34), have since kick-started a remarkable recovery which led to the downlisting of the species to “endangered” in the 2015 IUCN red list (32), and have propelled its population to 1,100 free-living lynxes in 2020 (35). Genetic erosion in the species is well documented (26, 31, 36, 37), with studies on mitochondrial, microsatellite, and genome-wide data providing evidence of a high inbreeding rate, high differentiation between the two populations, and a whole-genome species-wide nucleotide diversity (π = 0.026%) that is among the lowest reported thus far, and comparable to the single most genetically eroded population of the island fox (Urocyon littoralis) (38). Furthermore, there are documented signs of high inbreeding load in the species (39) and inbreeding depression in the DON population (36).

In contrast, the Eurasian lynx is one of the most widespread felids in the world, with a distribution area that extends from eastern Asia to central Europe and covers a wide range of habitats. It is currently classified at the global level as “least concern” by the IUCN. However, recently inferred demographic trajectories have revealed a generalized decline across its range during the last millennia (40). This decline has been particularly intense in the European part as a consequence of higher anthropogenic pressure (40), to the point where many of the western-most populations became extirpated during the 19th and 20th centuries, and some of the remaining ones, such as those in Norway (NOR) and the Białowieża and Knyszyn Primeval Forests in northeast Poland (POL), underwent severe bottlenecks and complete genetic isolation for decades (41, 42). Meanwhile, the Kirov (European Russia) population (KIR), which is part of a large continuous patch of the species in eastern Europe, has maintained a relatively large size and good connectivity with neighboring populations (40, 43). Even so, genetic diversity for KIR (π = 0.046%) is barely twice as high as that of the Iberian lynx and still below that of severely bottlenecked populations of other species, such as the Apennine brown bear, the Scandinavian wolverine, and the Amur tiger (40).

The contrasting demographic histories of the Iberian and the Eurasian lynx, coupled with the availability of a reference genome for the former species, showcases this dyad as a suitable model for carrying out genome-wide studies of the dynamics of the deleterious burden in bottlenecked populations. We specifically intend to address four questions: 1) whether there are empirically measurable differences in the amount or the patterns of the derived count between the Iberian lynx and the Eurasian lynx; 2) whether such differences extend to populations within each species; 3) whether differences in the derived count are dominated by the relaxation of purifying selection or by genetic purging; and 4) how these differences vary across genomic features, mutation effect categories, and genomic regions with contrasting efficiency of natural selection. To better interpret the observed patterns, we also extended the previous theoretical model on the purging of the inbreeding genetic load (13, 44) to develop analytical expressions predicting the evolution of the derived count, and we obtained theoretical projections under demographic scenarios similar to those inferred for the two lynx species. Beyond broadening our understanding of genetic processes in populations at the verge of extinction, such insights could be critical to predict the fate of small lynx populations, and to improve the efficacy of ongoing conservation programs.

Results

Genomic Diversity.

We analyzed 48 genomes sampled from the two remnant populations of the highly endangered Iberian lynx and from three demographically distinct populations of the widespread Eurasian lynx in order to unravel the dynamics of their deleterious burden. The mean genome-wide heterozygosity was 2.27 × 10−4 (±4.61 × 10−6) per base pair in the Iberian lynx (n = 20), which is roughly half that found in the Eurasian lynx (4.10 × 10−4 ± 7.10 × 10−6; n = 28). The genetic diversity of the populations (Fig. 1) ranked parallel to the available estimates of recent effective population sizes (31, 37, 40): KIR, 4.74 × 10−4 (±5.95 × 10−6); POL, 3.96 × 10−4 (±1.07 × 10−5); NOR, 3.28 × 10−4 (±1.76 × 10−5); AND, 2.75 × 10−4 (±6.46 × 10−6); and DON, 1.54 × 10−4 (±4.48 × 10−6).

Fig. 1.

Average nucleotide heterozygosity per annotation category and population, with error bars. Note the different ordinate scales between Upper and Lower. Blue, Eurasian lynx populations; red, Iberian lynx populations; dark tones, largest population of each species; light tones, smaller populations. syn., synonymous; mod., moderate.

Variant Annotation.

Bottlenecks and persistently small population sizes are predicted to alter the efficacy of selection, leaving an uneven footprint across sites under different pressures from natural selection. To analyze the effects of diverse demographic scenarios across populations, we considered several mutation categories that are expected to differ in their effects on fitness. Out of the total 4,245,741 high-confidence (i.e., postfiltering) biallelic single-nucleotide polymorphisms (SNPs) that we identified across all samples, we annotated 2,743,496 (64.62%) in intergenic regions (which are estimated to account for 58.08% of the callable genome after filtering), 1,404,222 (33.07%) in introns (34.31%), 53,182 (1.25%) in coding sequences (CDS; 2.22%), and 889 (0.02%) in ultraconserved noncoding elements (UCNEs; 0.09%). Within coding regions, SnpEff (45) identified 30,084 (56.6%) synonymous, 22,844 (43.0%) missense, and 560 (0.01%) loss-of-function (LoF) mutations. Among missense variants, 13,657 mutations with low predicted impact (Provean score > −2.5) were subclassified as tolerated, and 9,168 with high predicted impact (Provean score ≤ −2.5) were subclassified as deleterious. Similarly, UCNE variants were grouped according to sequence conservation-based Genomic Evolutionary Rate Profiling (GERP) scores (46) into two subcategories: moderate effect (moderately constrained; 2 < GERP ≤ 5), with 339 mutations, and high effect (highly constrained; GERP > 5), with 395 mutations.

Genetic Diversity within Annotation Categories.

Genetic diversity matched expectations of reduced heterozygosity in functional regions relative to neutral regions, particularly for the putatively more constrained categories (missense deleterious, LoF, and UCNE high effect). Although between-population differences were smaller in less diverse categories, the population ranking for genome-wide nucleotide diversity held across all annotation categories (Fig. 1 and ). Patterns of derived homozygosity were the inverse of diversity patterns in the case of neutral categories, but, for the putatively deleterious ones, differences between populations disappeared or even changed sign ().

Population Unfolded Site Frequency Spectra.

Focusing only on sites that are segregating within each species, we found a general trend of more right-skewed frequency distributions for mutations with a higher predicted impact on function (missense deleterious, LoF, and UCNE categories) compared to intergenic and intronic mutations (), as expected for more efficient selection against more deleterious alleles. In contrast, the unfolded site frequency spectra are roughly uniform across all categories in the highly bottlenecked DON population. Additionally, in the smaller populations of each species (NOR, POL, and DON), we found a higher proportion of derived alleles that are locally fixed compared to the larger populations (KIR and AND), particularly in the more neutral categories, as expected from classical theory (14).

Derived Count Estimates.

The derived count (i.e., the overall number of deleterious derived alleles per individual, both segregating and fixed) is a natural measurement of the deleterious burden. Using this statistic for each annotation category, we found similar average counts in the two species for both intergenic and introns, as expected for putatively neutral categories. In contrast, we found a deficit of LoF (∼29%) and missense deleterious (∼11%)—but not missense tolerated—mutations in Iberian lynx with respect to Eurasian lynx (Fig. 2 and Table 1), both of which remained significant even after applying Bonferroni and false discovery rate (FDR) adjustments (47). Results were nearly identical after applying a correction based on the intronic category (), and similar, although with slightly lower statistical significance, when correcting by the synonymous values (). Results for the high-confidence LoF subsets (which excluded mutations with annotation warnings as well as those in genes at the upper end of the mutation rate distribution) were nearly identical to those for the main LoF dataset (), suggesting this result is not driven by annotation errors or the inclusion of pseudogenes. Similar differences between species were also observed for one-to-one orthologs and genes intolerant to variation (i.e., two subsets of purportedly functionally relevant genes), although not for LoF mutations in one-to-one orthologs (). We also observed fewer UCNE high-effect mutations and an unexpected slight deficit of synonymous derived alleles in Iberian lynx with respect to Eurasian individuals, although neither reached statistical significance. Notably, populations within each species showed no derived count differences across any category (Fig. 2), despite their distinct recent demographic histories and genetic diversities.

Fig. 2.

Average number of derived alleles per annotation category. Population averages (represented by the dots) are all relative to the Kirov population average and shown with their error bars, with the darker tones representing the largest population of each species and the lighter tones representing the smaller populations. Species averages (represented by the horizontal solid lines) are relative to the Eurasian lynx average, with the colored ribbon representing the error. Asterisks mark the statistically significant comparisons between the two species. Blue, Eurasian lynx; red, Iberian lynx; m., missense; tol., tolerated; del., deleterious.

Table 1.

Summary of the statistical analysis of derived count differences between the Eurasian lynx and the Iberian lynx

Category	Mean derived count of Iberian relative to Eurasian lynx	Z value	P value (two-tailed)
Synonymous	0.967	1.097	0.273
Missense tolerated	0.999	0.026	0.979
Missense deleterious*	0.893	2.612	0.009
LoF*	0.710	2.659	0.007
UCNE	1.029	−0.221	0.824

The table shows, for the five tested genomic categories, the mean derived count of Iberian lynx relative to Eurasian lynx, the value of the Z test performed using the overall error of the species mean σ(M), and the resulting two-tailed P value.

The hypothesis of equal derived count for both species would be rejected in favor of the hypothesis that loads are different even after a 5% significance sequential Bonferroni multitest adjustment or using a 5% FDR.

Focusing on autosomal regions with contrasting recombination rates (), we found a deficit of missense deleterious mutations in Iberian lynx relative to Eurasian lynx in the highly recombining subtelomeric regions (although it did not remain significant after multiple test correction; ), but not in the low-recombination pericentromeric regions (). Moreover, we also found significant differences between the two species for the variable D, which results from subtracting the pericentromeric from the subtelomeric species average of the derived count. In other words, the reduction in the telomeric derived count compared to the centromeric derived count is larger for the Iberian lynx than for the Eurasian lynx (). For the X chromosome, the excess of missense deleterious mutations that we found in Iberian lynx () contrasts with the deficit observed in the genome-wide results, although the lack of contiguous blocks for the bootstrap prevented significance testing.

We also computed, for comparative purposes, other population statistics that have typically been used to estimate genetic load, such as the π/π and the P/P ratios (48), and the R statistic (22, 49). However, as we briefly discuss in , the first two () are highly dependent on demographic and methodological factors and do not necessarily reflect differences in selection regimes, and, in the case of the R ratio (), while it gives values consistent with our derived count results, its interpretation in terms of genetic theory is less straightforward.

Fixation Rate.

The Iberian lynx showed a significantly higher proportion of fixation of derived alleles across all neutral and putatively selected categories, with the exception of the LoF class, where the difference was nonsignificant (Fig. 3). As discussed below, this result is consistent with theoretical predictions, according to which purging is expected to induce quicker selection against deleterious alleles in the species with the smaller population size, but is not able to reduce their final fixation load driven by drift.

Fig. 3.

Average number of fixed derived alleles (with error bars) per annotation category and species relative to the Eurasian lynx values. Blue, Eurasian lynx; red, Iberian lynx.

Theoretical Predictions of the Derived Count.

The aim of these predictions is to provide a general framework to interpret the evolution of the derived count. They correspond to a simplified demographic model () inspired by the demography previously inferred for our two species (; see also figure 1 in ref. 26) where, after an initial common decay in effective size about 105 generations ago, one of the species undergoes a series of bottlenecks modeled after those of the Iberian lynx, while the other one is assumed to persist at a large constant size (a qualitatively informative simplification of the Eurasian lynx history). These analytical results should therefore be interpreted as a qualitative illustration of the dynamics of the derived count after a bottleneck, rather than as precise quantitative predictions for these species. To illustrate the characteristics of the mutations that could be purged in the putatively deleterious categories under these demographic scenarios, predictions were computed for different combinations of representative values of the homozygous deleterious effect (s) and the dominance coefficient (h). The per gamete mutation rate (λ) was always set to one so that, to adjust predictions considering a particular joint distribution for s and h, each predicted burden should be multiplied by the gametic deleterious mutation rate considered more plausible for each s,h class. While it is unlikely that any single s,h pair could adequately represent the distribution of deleterious effects within any of the deleterious categories considered in our genomic analysis, these predictions can help us to understand some characteristics of such effect distributions. Notwithstanding, at least in the case of the putatively more deleterious categories (missense deleterious and LoF), it could be hypothesized that they are enriched in (partially) recessive and moderately to strongly deleterious mutations.

Modeling the expected dynamics of the derived count (or deleterious burden) of the two species showed that, after a bottleneck, purging only produces a large reduction of the derived count of the smaller population (relative to that of the larger one) for deleterious mutations with h ≤ 0.25 (Fig. 4). Conversely, for mutations that are roughly additive (h = 0.45), reductions in the effective size (N) such that N < 5 result in increased burden accumulation in the long term. Within the represented time scale, the increase in burden for roughly additive mutations is smaller than the reduction from purging observed for recessive deleterious mutations (h = 0.05). Importantly, the latter contribute a much greater derived count relative to their mutation rate than those that are nearly additive (note the different scaling of the ordinates axis in Fig. 4).

Fig. 4.

Predicted evolution of the derived count against the number of generations before present (in decimal logarithm) for two populations derived from a common ancestor that was close to the mutation–selection balance (effective size n = 106). Thin line, Iberian lynx; thick line, Eurasian lynx. The demography of both populations, represented in , is modeled after that of the two lynx species analyzed here. The effective sizes after each bottleneck are shown in the first panel. Panels reflect different combinations of h,s values, always assuming a per gamete mutation rate of λ = 1.

Discussion

In this study, we analyzed whole-genome sequences of the endangered Iberian lynx and the widespread Eurasian lynx, sister species with historical differences in effective population size and highly contrasting recent demographies and distribution ranges, to compare patterns of deleterious burden at the species and population levels. We observed a significant depletion in the number of derived deleterious alleles per individual (adding together both segregating and fixed alleles), along with lower genetic diversity and a higher fixation rate, in the Iberian lynx relative to its sister species. Populations within each species showed no differences in overall derived count levels, but we observed higher homozygosity and fixation rate in the smallest populations.

Species-Level Analysis.

Our main finding was a significantly smaller derived count in the endangered Iberian lynx than in the widespread Eurasian lynx for the two main putatively deleterious categories—missense deleterious and LoF—which contrasts with the absence of differences for the putatively neutral categories.

As Fig. 4 illustrates, a smaller derived count in the historically smaller population should be ascribed to genetic purging of deleterious alleles with high recessivity. These purging-driven differences in the derived count between the two species are a consequence of the higher rate of selection against homozygotes following a more pronounced reduction in effective size and, as shown in the figure, are not restricted to very small populations. However, since the ultimate fixation rate of both additive and recessive deleterious alleles drops with N (14, 15), for small enough populations (say N < 5/s), the continuous fixation of deleterious mutations will, at some point, overcome the reduction of burden from purging. This point is not reached in Fig. 4, due to the relatively short periods during which the declining population had very small N values, but it is illustrated by the example presented in Fig. 5. Therefore, a larger derived count in a smaller population does not exclude the possibility that purging may have occurred in the past, mitigating, to some extent, its fitness inbreeding depression. Fig. 5 also shows that the transiently smaller derived count of the smaller population does not translate into a higher fitness because of the concomitant increase in homozygosis. Nevertheless, increased purging in the smaller population will result in some reduction of its fitness inbreeding depression and its inbreeding load (which is responsible for the inbreeding depression that could occur in the future) and, therefore, in a somewhat increased resilience to subsequent bottlenecks. For roughly additive mutations, when N < 5, increased drift in the smaller population leads to a larger derived count, which accumulates very slowly but indefinitely (note that, in Fig. 4, the apparent deceleration of this increase is due to the logarithmic scale used for generations). In contrast, for recessive mutations, purging effects are limited to the transition period from the ancestral to the new mutation–selection–drift (MSD) balance, which is approached faster for larger s(1/2 – h) and for smaller N values (see ref. 13). Thus, the detection of differences in derived count due to purging depends on the specific demographies and on the length of the periods involved.

Fig. 5.

Example illustrating the role of purging and drift in determining the evolution of the derived count and the fitness average through generations for two populations with effective sizes N = 104 and 2,500 (thick and thin lines, respectively) derived from a common population at the MSD balance with N = 106. (A) Evolution of the haploid derived count. (B) Evolution of the average fitness: Solid lines give predictions matching those of the derived count in A (i.e., computed considering purging, mutation, and standard selection by using equation 14 from the Full Model approach in ref. 13); dashed lines give analogous fitness predictions assuming no purging. Per gamete mutation rate of λ = 0.1, s = 0.001, h = 0.05.

As most deleterious mutations are thought to have small effects, and since it is usually accepted that mutations with smaller s values tend to be less recessive, it may seem unexpected that purging causes a substantial reduction of the derived count. However, the estimates of the dependence of expected h on s are mainly based on the analysis of mutations whose effect can be experimentally detected (50–52). Information about h for deleterious mutations with s ≤ 10−2 is scarce, with a recent study suggesting h close to zero for s > 4 × 10−4 (53) and thus effectively trivializing additive gene action. Furthermore, even if most mildly deleterious new mutations were additive, the smaller share of recessive ones could still make a large relative contribution to burden, as Fig. 4 illustrates, and lead to a net reduction in the overall derived count following a population bottleneck.

Our main result should thus be interpreted as evidence of genetic purging in the Iberian lynx. Beyond the prevailing recessivity within the more deleterious categories (which is not at odds with many new mutations being roughly additive), these results also imply some historical reduction of 1) the inbreeding depression in the Iberian lynx, which may have enabled its long-term persistence despite extensive genetic erosion, and 2) its inbreeding load responsible for future inbreeding depression.

Although well supported by evolutionary experiments (54–59), and simulations (e.g., refs. 60–63), including some conducted for wild bottlenecked populations under plausible models of demographic history (21, 23), the purging of deleterious variation in small populations has only been explored from a genomics standpoint in the last few years. In humans, Narasimhan et al. (49) found significantly fewer LoF alleles per individual in Finns and other newly founded populations when compared to nonbottlenecked populations. Similarly, Xue et al. (22) observed a significantly lower burden of LoF mutations per individual in the bottlenecked eastern species of gorilla (Gorilla beringei) than in its western counterpart (Gorilla gorilla). Grossen et al. (21) recently found signs of increased purging in the most bottlenecked Alpine ibex (Capra ibex) population, although significance testing was performed only for the accompanying simulations. In contrast, in our study, the differences in deleterious burden between the Iberian and the Eurasian lynx are backed by a statistical approach that incorporates both the error due to individual sampling and evolutionary stochasticity and is further underpinned by multiple testing correction procedures. The conspicuous purging in Iberian lynx, possibly among the strongest reported so far in the literature, coupled with its low levels of genetic diversity, is fully consistent with the long-term small population sizes and serial bottlenecks reported for the species (26), as shown by our predictions. Notwithstanding, it should be noted that LoF mutations are often enriched in mapping, genotyping, and annotation errors compared to other variants (64); thus, our derived count results for this category must be viewed with some caution despite their statistical significance and the reassuring outcomes of the additional validation tests. Nevertheless, the significant purging effects detected for the missense deleterious category constitute a particularly reliable and relevant finding.

To further test the robustness of our observations of augmented purging in Iberian lynx, we analyzed genome compartments with disparate predicted efficiency of selection due to differences in either recombination rate or sex-dependent ploidy. Briefly, the lower recombination rates in pericentromeric regions (65) can slow down the action of purging (66), compared to high recombination subtelomeric regions (65). Accordingly, we did not find any difference between the two species for pericentromeric derived count, but the burden of missense deleterious mutations in subtelomeric regions was significantly smaller in the Iberian lynx than in the Eurasian lynx. On the other hand, we did not observe any deficit of missense deleterious mutations in the Iberian lynx relative to the Eurasian lynx in the X chromosome. In mammals, the (partially) recessive deleterious alleles in this chromosome are continuously exposed to selection in males regardless of population size, so that there should be very little burden ascribed to substantially recessive deleterious alleles to be exposed by inbreeding (67, 68). Thus, purging is not expected to produce appreciable differences between species in the deleterious burden of this chromosome. Taken together, these results indicate that the reduced burden in the Iberian lynx is associated with regions where a more effective purging is predicted and, at the same time, is not observed where the derived count ascribed to recessive alleles is expected to be negligible in both species. These observations reinforce the purging hypothesis.

Deleterious effects in the missense tolerated category are expected to be less severe than in the missense deleterious one, and, in this study, their patterns are virtually identical to those of the neutral categories. A parsimonious explanation for this result is that the evolution of missense tolerated mutations is driven mainly by genetic drift in both species, which suggests, on the basis of theory (14), that their s values are usually, at most, on the order of the inverse of the historical effective size (s ≤ 1/N) for the Eurasian lynx (i.e., the most abundant species), estimated to be below 20,000 in the last ∼50,000 y (26). There are other alternative explanations for this result. Missense tolerated mutations could have roughly additive effects large enough to allow efficient selection that leads to similar counts in both species, although this seems unexpected for such a putatively less severe category. Alternatively, the accumulation of additive deleterious mutations behaving as neutral in the smaller population could be compensated by a serendipitously similar fraction of high-effect recessive variants being purged. The starkly different results that we obtained for the missense subcategories illustrate the wide variability in dominance and selection coefficients among these mutations, and highlight the pertinence of the more fine-grained approaches that we have applied to this (as well as the UCNE) category, as previous studies have suggested (10, 69) or have shown (e.g., refs. 20 and 21).

Regarding the UCNE category, where most mutations are expected to be strongly deleterious, the number of variants is so small that SEs are too large to allow any relevant conclusions. That said, in the case of the UCNE high-effect subcategory, which is enriched in highly conserved sites, the relative deficit of derived copies in Iberian lynx (although nonsignificant) is also suggestive of purging, whereas no such reduction is observed in UCNE moderate-effect sites. Surprisingly, we also found a slight relative deficit of derived copies in Iberian lynx for the synonymous category. Although there is increasing evidence that some synonymous mutations have deleterious effects (70–72), these are expected to be small, and, to our knowledge, ubiquitous purging within this category has not been described in the literature thus far. Rather, neutral processes such as GC-biased gene conversion have been identified as the main driver of synonymous variation patterns in mammals with small populations (73, 74). Since differences between the two lynx species for the synonymous category are nonsignificant, they can nevertheless be attributed to stochasticity.

Within-Species Population-Level Analysis.

Contrary to what we observed in the species-level comparison, we did not find significant differences in the derived count between populations of the same species, despite their distinct recent demographic histories. There are some plausible explanations for this. First, the differences in recent demographic history between the sampled populations remain poorly known beyond their general trend of decline (particularly in the case of the Eurasian lynx populations), including the number, intensity, and duration of the bottlenecks, the patterns of migration and isolation, and the possibility of sporadic interspersed growth periods—and some of these events have opposing expected effects on the deleterious burden. In the case of the Iberian lynx, purging might have been boosted by increased fragmentation during the period of accelerated demographic decline, although this may have occurred mostly before the isolation of its two remnant populations. Second, the joint effect of genetic drift and selection is appreciably driven by the latter only for alleles whose homozygous deleterious effect is at least on the order of the inverse of the historical effective population size: s ≥ 1/N (14). Among these alleles, only relatively recessive ones can be selected more efficiently in smaller populations than in larger ones. Given the above, and the recency of the small size of most of our sampled populations, purging can reduce their burden below that of a larger population only for considerably recessive deleterious alleles that are relatively severe. For example, out of the cases considered in Fig. 4, only s = 0.1 leads to N > 1 in DON during its recent isolation period, while drift is expected to be the leading evolutionary force in the remaining cases. However, deleterious mutations with s ≥ 0.1 are expected to contribute a very small fraction of the derived count (75–77). Therefore, it would have been rather unexpected to detect a burden reduction in DON relative to AND. In contrast, the detection of purging at the species level suggests that many mutations in the putatively more deleterious categories have deleterious effects below 10%, and that a possibly small but still relevant fraction of them are substantially recessive (h < 0.25). Our predictions suggest that these effects are generally on the order of 1% or below, although effects between 1% and 10% could also be relevant, particularly if fragmentation during the Iberian lynx demographic decline played a role in promoting purging at the species level. In addition, the absence of differential purging at the population level reinforces the conclusion that the genome-wide purging observed in Iberian lynx took place before its two remnant populations became isolated, and can be attributed to its smaller long-term effective size due to the historical species-wide bottlenecks (26).

Despite the lack of derived count differences within each species, the smaller populations did show a more uniform distribution of frequencies, less diversity, higher homozygosity, and increased fixation compared to the larger ones. Although differences were much smaller for deleterious than for putatively neutral categories, these patterns, in particular, the increased burden in homozygosis and the presence of locally fixed deleterious alleles, might be indicative of some additional fitness deterioration in these smaller populations. In fact, although strongly deleterious mutations make a small contribution to the derived count, they can make a substantial contribution to the fitness inbreeding depression and inbreeding load (78). Therefore, a tiny increase in the number of deleterious alleles in homozygosis can produce an important fitness decline. In connection with this, there is evidence of substantial decay of fitness traits in the DON population during the recent period of isolation in the form of a decrease in survival and litter size with time, concomitant with an increase of homozygosity (36), which is suggestive of inbreeding depression. Since no such signs have been observed in the AND population, its decline may have been slow enough to have allowed additional (albeit undetectable) purging against severely deleterious alleles in the initial period after the isolation of both populations, which would be consistent with their reconstructed demography (31). At the same time, as noted above, a large share of the burden historically constrained by selection in the Iberian lynx has likely experienced drift-driven evolution during the very recent past, with severe population declines. Thus, it may have induced some level of undetectable species-wide fitness inbreeding depression when compared to the Eurasian lynx.

Together, these processes that are expected to occur at the population and species level imply the existence of a gap between the prediction of the fitness inbreeding load (and inbreeding depression) and the burden of deleterious alleles assessed in genomic analyses, even after their partition into categories according to their putative severity. It has recently been proposed that, since neutral genetic variation has a nonstraightforward association with adaptive genetic variation, its use in conservation management should be replaced with that of mapped functional genomic variants (ref. 79, but see refs. 80 and 81). However, the aforementioned gap illustrates the possible consequences of our inability to assess the magnitude of fitness effects for putatively deleterious alleles. While this problem should be progressively reduced at least for large-effect mutations [for example, results from knockout experiments using model organisms may provide rough estimates of the deleterious effects of LoF mutations (82)], all in all, the existence of such a gap calls for caution when using genomic information to assist conservation decisions, and advocates an integrative approach that incorporates careful quantitative fitness analyses together with the assessment of both the neutral genetic diversity (including the reconstruction of the demographic history) and the burden of deleterious alleles (81).

Implications for Conservation.

Following the recent description of the retention of key functional variation in the Iberian lynx by interlocus balancing selection (83), the genetic purging of deleterious variation likely constitutes yet another important mechanism to explain its long-term persistence, but several signs of extensive genetic erosion are still present in the species. These include values of overall nucleotide diversity well below the already low ones of the Eurasian lynx, which are consistent with previously reported measures of genome-wide diversity (26, 40), and larger homozygosity and fixation rates. In connection with this, the positive correlation found in the Iberian lynx between the proportion of abnormal sperm and the degree of individual inbreeding (39) indicates that the inbreeding load remaining in the species is still a concern. Moreover, several deleterious traits with a likely genetic basis have been observed in recent years (84–86). These signs of inbreeding depression emphasize the importance of the current conservation measures aimed at reducing inbreeding and increasing the effective size of the Iberian lynx. In particular, the admixture in captivity (34) of both remnant populations, AND and DON, should help alleviate the accumulation of inbreeding and facilitate the purging of deleterious alleles which may have become fixed in one of the two populations in recent times. Additionally, translocation in the wild (33) is not expected to introduce a dramatic amount of load in any of the two populations considering their recent demography (31). In any case, reconnection in the wild through either translocation or, ideally, restoration of natural connectivity should improve the species survival (87). Notwithstanding, deleterious alleles fixed at the species level can no longer be removed from these populations by increasing their size, which leaves some concern over their long-term fitness, and highlights that maintaining large population sizes should always be a conservation priority.

The larger derived count observed for the Eurasian lynx populations (compared to Iberian lynx) in deleterious categories is consistent with the higher fitness inbreeding load expected for larger outbred populations (44, 88). Although this inbreeding load should not be expressed in these populations, it alerts us to the possible risk of future inbreeding depression in a potential program using individuals from the largest populations to rescue endangered isolated conspecific populations (where a few of the most severely deleterious mutations could have been purged in the past with imperceptible impact on the overall deleterious burden). This risk derived from introducing load during rescue is conveniently illustrated by the well-documented case of extinction of the Isle Royale wolves, ascribed to the inbreeding load introduced by a single continental migrant male (89–91). The consequences of introducing inbreeding load by prioritizing large nonpurged donor populations with high genetic diversity in rescue programs have been noted by Kyriazis et al. (ref. 92, but see ref. 93 for criticisms). In fact, risks can be even more severe than warned by these authors (due to the unrealistic distribution of the deleterious effects assumed in their simulations), particularly under scenarios where the size of the rescued population allowed for purging in the past but is presently too small to enable the purging of the introduced load (87). Prior studies have called upon the need to establish integrative conservation actions for several European populations of the Eurasian lynx (40, 94, 95). On the basis of their declining demographic trends and their random differentiation following anthropogenic fragmentation and isolation (40), we argue that connectivity should be preventively favored between small, likely purged populations, and translocations from large populations should be initiated only when reasonable prospects for rapid growth exist.

Beyond the quantitative analyses of deleterious burden and patterns of deleterious variation in the Iberian and the Eurasian lynx, this work provides a genome-wide database of over 4 million derived mutations, of which several thousand are potentially deleterious and likely deserve further research. Incidentally, some deleterious traits believed to have a genetic basis, such as idiopathic epilepsy, membranous glomerulonephritis, and cryptorchidism, are currently relatively common in the Iberian lynx (84, 86), even in the captive population, despite its high level of admixture (34). Further investigation of these traits through genome-wide association studies and candidate gene approaches may ultimately lead to the identification of variants of major effect, whose selection could be eventually integrated into the ongoing genetic management.

Materials and Methods

Sampling and Genomic Datasets.

Our sample included 31 Iberian lynx individuals from the last two remnant populations of AND (n = 19) and DON (n = 12), in southern Spain, before they were reconnected through translocations and reintroductions, and 29 Eurasian lynx from KIR (n = 13), NOR (n = 8), and POL (n = 8). For details on the demography of these populations, please refer to . Additionally, we included one individual bobcat (Lynx rufus), from the Jerez Zoo (Spain). A summary of the sample is provided in . We arranged samples in four datasets for depth-filtering purposes. Our main analyses were carried out on 48 samples (12 each from AND and KIR, and eight each from DON, NOR, and POL) resequenced at a depth of ∼5.75× (range: 5× to 7×) using Illumina HiSeq2000 v3 and v4. See for further details on each dataset.

DNA Extraction, Library Preparation, and Sequencing.

DNA was extracted from good-quality tissue or blood samples. Paired-end libraries were prepared, and then sequenced in a fraction of a sequencing lane of HiSeq2000 flowcell v3 (Illumina Inc.), as detailed in .

Read Processing and Mapping.

Read processing and mapping were mostly performed as described elsewhere (40), and are further detailed in . Trimmed reads were mapped to the Iberian lynx reference genome using BWA-MEM (96) with default parameters.

Variant Calling, Polarization, and Filtering.

GATK v3.7 HaplotypeCaller (97) was used to generate variation data for the pool of 60 Iberian and Eurasian lynx samples. Repetitive, low-complexity, and low-mappability regions, which together account for 46.24% of the genome (26), were excluded from the calling. To polarize variants into the ancestral and derived states, we applied parsimony criteria using information from three outgroup species: bobcat, domestic cat (Felis catus), and tiger (Panthera tigris). For the bobcat, we took the available whole-genome mapping against the Iberian lynx reference genome (40), while, for the latter two species, the information for our final set of filtered variants was retrieved from the available Iberian lynx–cat–tiger synteny (26). Finally, a stepwise filtering framework was designed to keep only high-confidence biallelic SNPs. See for additional details.

We used SnpEff v4.3i (45) to annotate variants based on a custom annotation file (). This resource enabled the classification of variants into intergenic, intronic, CDS, and UCNE categories. The ability of SnpEff to predict the functional effects of variants allowed us to define three further mutation categories within CDS: synonymous, missense (i.e., nonsynonymous), and LoF. The latter group comprises variants whose predicted impact on function is classified by SnpEff as high effect, including stop-gained, start-gained, start lost, splice acceptor, and splice donor variants. Because LoF mutations are usually enriched in errors relative to other mutation types (64), we also identified high-confidence mutations to validate our results for this category (). In addition, we split missense variants on the basis of their degree of conservation, the biochemical distance between the encoded amino acids, and the sequence context, as inferred by Provean v1.1.5 (98), into “missense tolerated” (Provean score > −2.5) and “missense deleterious” (Provean score ≤ −2.5). Similarly, variants in regions annotated as UCNE were categorized according to their presumed evolutionary constraint, estimated through GERP scores (46) that were remapped from human to lynx, into “UCNE moderate effect” (2 < GERP ≤ 5) and “UCNE high effect” (GERP > 5). See for further details.

Genomic Variables Based on Individual Counts.

The alternative allele count was recorded across sites and annotation categories to obtain the individual counts of derived alleles (i.e., the derived count), heterozygous derived alleles, and homozygous derived alleles. Similarly, individual counts of fixed derived alleles were tallied to compare the fixation rates between the two species. For each of these four variables, we then calculated the population arithmetic mean of the individual counts, and the species average of the corresponding population means weighted by their corresponding sampling sizes. We also recorded the derived count per individual for two subsets of (purportedly) functionally relevant genes, as well as the autosomal subtelomeric and pericentromeric regions, and the X chromosome (). Average per-site estimates of heterozygosity were also estimated ().

Statistical Error and Significance Testing.

In order to estimate the statistical errors of the population and species averages for the different genomic variables, which allow testing of whether the observed differences between them are significant, we developed a resampling simulation procedure which takes into account both the sampling variance within each population or species and the evolutionary stochasticity ().

Here we consider the model developed by García-Dorado (13, 44) that predicts the inbreeding load at the MSD balance and its evolution after a reduction in population size, and we extend it to obtain predictions for the expected value and evolution of the per gamete (i.e., haploid) derived count of deleterious alleles (see for the extended description of the method). It assumes an ancestral population with a relatively large effective size (N) and predicts the overall count of segregating and fixed deleterious alleles through generations after the effective size reduces to N, including both the deleterious alleles that were segregating in the ancestral population and those arising due to continuous mutation as the population approaches a new MSD balance. This method involves several approximations (as explained in ) but provides very illustrative predictions. To improve our understanding of the contribution of deleterious mutations of different effects to the evolution of the derived count of the different putatively deleterious categories, we computed these predictions for a set of illustrative s,h values representing mutations of different severity (s = 0.001, 0.01, and 0.1) and dominance (h). Since h ≈ 0 can lead to huge amounts of hidden burden for very large populations and, therefore, to very large contributions of purging to the evolution of the overall derived count, we considered the more conservative h = 0.05 value and, for the sake of symmetry, h = 0.25 and 0.45.