Genomic signatures of extensive inbreeding in Isle Royale wolves, a population on the threshold of extinction

See allHide authors and affiliations

Science Advances  29 May 2019:
Vol. 5, no. 5, eaau0757
DOI: 10.1126/sciadv.aau0757


The observation that small isolated populations often suffer reduced fitness from inbreeding depression has guided conservation theory and practice for decades. However, investigating the genome-wide dynamics associated with inbreeding depression in natural populations is only now feasible with relatively inexpensive sequencing technology and annotated reference genomes. To characterize the genome-wide effects of intense inbreeding and isolation, we performed whole-genome sequencing and morphological analysis of an iconic inbred population, the gray wolves (Canis lupus) of Isle Royale. Through population genetic simulations and comparison with wolf genomes from a variety of demographic histories, we find evidence that severe inbreeding depression in this population is due to increased homozygosity of strongly deleterious recessive mutations. Our results have particular relevance in light of the recent translocation of wolves from the mainland to Isle Royale, as well as broader implications for management of genetic variation in the fragmented landscape of the modern world.


Under increasing human population pressure, many species with once continuous ranges have been reduced to small fragmented populations (1). Higher levels of inbreeding in these small populations may elevate the risk of extinction through inbreeding depression [e.g., (2, 3)]. Large carnivores are particularly susceptible, since they typically have lower population densities relative to the herbivores they prey on and they often require extensive natural areas to persist. Further, they are frequently persecuted because of real or perceived threats to humans and livestock (4). Some well-known examples of inbreeding depression in the wild have been observed in large carnivores, such as Florida panthers [Puma concolor (5)] and gray wolves (68). Studies of inbred wolves in the wild and in captivity have found elevated rates of blindness, cryptorchidism, heart and kidney defects, dental anomalies, and vertebral malformations, as well as decreased reproduction, litter size, body weight, and longevity (6, 913). Phenotypic defects associated with inbreeding are not limited to carnivores, however, and have been observed in numerous plant and animal species (14). Quantifying and maintaining genetic diversity to minimize the risk of inbreeding depression therefore remain a fundamental goal of conservation biology.

Studying the genetic basis of inbreeding depression remains challenging, particularly in wild populations. Previous studies largely support the hypothesis that inbreeding depression is caused by the increased homozygosity of recessive strongly deleterious alleles [reviewed in (15)]. However, the adverse consequences of small population size have been debated, in part, because theory predicts that strongly deleterious alleles are more likely to be purged from small populations (16). With the ever-decreasing costs of whole-genome sequencing, it is now feasible to estimate the genome-wide burden of deleterious variants (genetic load). Thus far, genomic studies of deleterious variation have primarily dealt with the effects of ancient bottlenecks and long-term reduced population size [e.g., (1721)] or population expansion (22), rather than inbreeding in small populations. In general, small increases in the number of deleterious alleles per genome (additive genetic load) have been observed in these studies, but the effects of recent inbreeding within the past few generations are likely to be distinct.

In this study, we present results from complete genome sequencing of a small, highly inbred population of wolves on Isle Royale, MI, USA. The island was colonized by two to three mainland wolves that crossed frozen Lake Superior in the 1940s, establishing a population that, at its peak, included approximately 50 individuals (23). Following a disease outbreak in the early 1980s, the population crashed to 14 individuals and failed to rebound for approximately 15 years. The population exhibited a substantial but short-lived improvement in numbers after a single wolf migrated from the mainland in 1997 before falling to even lower numbers by 2010 (2325). Severe inbreeding in Isle Royale wolves has resulted in conspicuous morphologic defects, especially malformed vertebrae (8). Other abnormalities have also been observed in Isle Royale wolves, including syndactyly (8), probable cataracts, an unusual “rope tail,” and anomalous fur phenotypes. A previous analysis revealed that highly inbred wolves on Isle Royale had reduced survival and reproduction (25) and the population has fallen from 30 individuals to only two over the past 12 years (26). Meanwhile, the number of moose on Isle Royale, the wolves’ main prey, has swelled from 450 to 1600 individuals over the same period (26), implying that the demise of the wolf population cannot be attributed to the lack of available prey. A pair of individuals, both father-daughter and half-siblings, descended from a legacy of repeated close inbreeding events were the only wolves that remained by early 2018, by which time they were either in or approaching senescence. No successful reproduction has occurred since 2014 (26), prompting the National Park Service to initiate a plan to introduce 20 to 30 wolves to the island over the next few years.

Although previous investigations of Isle Royale wolves focused on inbreeding using genetic assays (24, 25, 27) or morphological assessment (8), our study combines both approaches and uses complete genome sequence data. Recent genomic studies of isolated European wolf populations have found that inbreeding is associated with the formation of large runs of homozygosity (ROH) (28, 29), which are chromosomal segments inherited identically by descent from a recent common ancestor through both the maternal and paternal lineages. However, the extent to which these large ROH are associated with inbreeding depression or deleterious variation remains unclear. Analysis of putatively damaging mutations within ROH in inbred populations suffering from reduced fitness has the potential to reveal the underlying basis of inbreeding depression. Here, we analyzed patterns of genome-wide diversity, deleterious variation, and skeletal phenotypes in a sample of Isle Royale wolves to gain insight into the severe inbreeding depression and collapse of this population. Through comparison with a sample of wolf genomes worldwide and through population genetic simulations, we find that long-term demographic history is the key factor driving patterns of heterozygosity and deleterious variation across the genome. These results have implications for the effective management of small populations recently derived from large outbred populations, as is the case for many species threatened by habitat changes resulting from recent human impacts on the landscape.


Genome-wide patterns of variation shaped by demographic history

We obtained genetic samples from 11 Isle Royale wolves collected between 1988 and 2012 for whole-genome sequencing and analysis (table S1). The population size on the island was estimated to number 8 to 30 individuals during this period (26). A pedigree, adapted from Hedrick et al. (25), shows the relationships of the sequenced wolves (where known) and inbreeding events between close relatives (Fig. 1A). We supplemented the Isle Royale wolf genomes with publicly available and newly generated sequences to include six mainland wolves from nearby Minnesota, nine gray wolves from elsewhere in North America, six Eurasian gray wolves, and a single genome from each of the following species: red wolf (Canis rufus), coyote (Canis latrans), and Ethiopian wolf (Canis simensis) (Fig. 1B and table S1) (30, 31). Our dataset spans the Holarctic range of the gray wolf and contains individuals derived from a variety of demographic histories that feature recent inbreeding, long-term small population size, isolation, and admixture (table S2). To our knowledge, Isle Royale and Mexican wolves are the only populations we sampled with documented evidence of inbreeding depression (8, 11). All genomes were aligned, genotyped, and annotated with respect to the domestic dog reference genome (canFam3.1), yielding mean genome-wide coverage values of 12× to 49× after read filtering (table S1). We determined the ancestral allele at each single-nucleotide polymorphism (SNP) by including two outgroup genomes from an African golden wolf (Canis anthus) and a gray fox (Urocyon cinereoargenteus) (see Materials and Methods) (20, 32). A cladogram depicting the phylogenetic relationships between wolf populations and sister taxa is shown in fig. S1.

Fig. 1 Isle Royale wolf pedigree and geographic origins of genomes in this study.

(A) Pedigree of Isle Royale wolves sequenced in this study (numbered individuals), adapted from Hedrick et al. (25). Circles represent females, and squares represent males. Relationships were inferred from genotypes at 18 microsatellite loci. Shaded individuals were examined for the presence of vertebral abnormalities for this study (see Table 1). The ancestries of F25, F55, M61, and F67 are unknown, but F67 is known to be the mother of F65 and M175. M was a wolf that migrated from the mainland in 1997 (24). The “A” individuals were the last two wolves alive on Isle Royale in 2018. (B) Map showing approximate origins of all individuals analyzed in this study. BI, Baffin Island; EI, Ellesmere Island; NU, Nunavut; VI, Victoria Island. See table S1 for further sample information.

To assess how demographic history has shaped patterns of genomic diversity in wolves, we calculated per-site heterozygosity in nonoverlapping 1-Mb windows across the genome within each individual (Fig. 2 and fig. S2). Qualitatively, we observed three distinct patterns of genome-wide heterozygosity: (i) genomes with high heterozygosity throughout (Fig. 2, A and B), (ii) genomes with low heterozygosity throughout (Fig. 2, C and D), and (iii) genomes with a sawtooth-like pattern characterized by regions of high heterozygosity interspersed by long ROH devoid of variation (Fig. 2, E and F). High heterozygosity across the entire genome was present in individuals derived from outbred populations with large long-term effective population sizes, such as lowland Chinese (Xinjiang) and Minnesota wolves (31, 3335). Conversely, low heterozygosity across the entire genome was associated with a history of isolation and long-term low effective population size, such as in Ethiopian and Tibetan wolves (30, 31, 36, 37). Last, the sawtooth-like pattern characterized by long ROH interspersed with regions of high heterozygosity was observed in individuals with a history of recent inbreeding, as in the Isle Royale and Mexican wolves (11, 25). In some cases, ROH spanned entire chromosomes, consistent with predictions from a simulation-based study of the Isle Royale wolf pedigree (38). ROH spanning entire or near-entire chromosomes were also observed in recent studies of highly inbred Scandinavian and Iberian wolf genomes (28, 29). In total, we found that 98.2% of the autosomal genome falls within ROH in at least one Isle Royale wolf and the other 1.8% contains genes with significantly reduced (or no) sequence coverage (P = 4.16 × 10−13; fig. S3).

Fig. 2 Distributions of heterozygosity across the genome under different demographic histories.

In each panel: (left) Example barplots showing per-site heterozygosity in nonoverlapping 1-Mb windows across the autosomal genome; (right) histograms of per-window heterozygosity. (A and B) The Xinjiang and Minnesota wolves represent large outbred populations. (C and D) The Tibetan and Ethiopian wolves represent small isolated populations. (E and F) The Mexican and Isle Royale wolves represent populations with recent inbreeding. See fig. S2 for plots of all individuals. Het., heterozygotes.

Overall, genome-wide heterozygosity in Isle Royale wolves (0.96 to 1.47 heterozygotes/kb) is 11 to 42% lower than the mean heterozygosity of wolves from nearby mainland Minnesota (1.65 heterozygotes/kb) (Fig. 3A). Although mean genome-wide heterozygosity correlates with the total amount of the genome within ROH, it does not capture the distribution of ROH lengths, which is shaped by demographic history (Fig. 3A). Specifically, long ROH result from recent close inbreeding, whereas shorter ROH reflect ancient shared ancestry between parental lineages. Isle Royale wolves contain 12 to 39 long ROH (>10 Mb) per individual, whereas Minnesota wolves contain just 1 to 10 long ROH per individual (fig. S4). Long ROH are also prevalent in the genomes of the inbred Mexican wolf and the Ellesmere Island wolf, the latter also suspected to be inbred (Fig. 3A and fig S4) (39). The Mexican wolf genome derives from the highly inbred Ghost Ranch lineage before measures were taken to decrease inbreeding within the Mexican wolf captive breeding program. In contrast, the Tibetan and Ethiopian wolf genomes contain many shorter ROH but few long ROH, reflecting the historical isolation and small effective sizes in these populations (30, 31, 36, 37). Notably, these wolves have lower genome-wide heterozygosity overall than the Isle Royale wolves (Fig. 3A). The increased prevalence of long rather than shorter ROH in Isle Royale wolf genomes is consistent with their recent descent from mainland wolves, which bear few ROH of any length and have high heterozygosity across their entire genomes, followed by severe inbreeding.

Fig. 3 ROH and the correspondence with genome-wide heterozygosity and FPED.

(A) Left: Per-site autosomal heterozygosity across the autosomal genome. Samples are ordered by decreasing heterozygosity from top to bottom. Right: Summed lengths of short (0.1 Mb ≤ ROH < 1 Mb), medium (1 Mb ≤ ROH < 10 Mb), and long (10 Mb ≤ ROH < 100 Mb) ROH per individual. (B) FROH is the proportion of the genome contained within ROH of ≥100 kb. FPED values were calculated from the pedigree of Hedrick et al. (25) (Fig. 1A). The gray dashed line shows the diagonal y = x. (C) FROH in Isle Royale wolves over time (by birth year; see table S1) reflects the effects of increasing inbreeding over time before and after the arrival of the migrant individual in 1997 (gray dashed line; see Fig. 1A). Gray bars indicate the census population size on Isle Royale. (B and C) Open circles denote individuals not related to the migrant, and closed circles denote individuals descended from the migrant (according to the pedigree).

Pedigree-based inbreeding coefficients (FPED) provide the expected proportion of the genome that is contained within ROH (FROH). We calculated that 23 to 48% of the Isle Royale wolf genomes were contained within ROH of 100 kb or longer, which is significantly greater than FROH in Minnesota wolves who carry 12 to 24% of their genomes within ROH [one-tailed Mann-Whitney U (MWU) test, P = 1.62 × 10−4]. Among the pedigreed Isle Royale wolves in our dataset (n = 7), FPED ranged from 0 to 37.5%. Since the pedigree only includes Isle Royale wolves genotyped between 1998 and 2013, there are several individuals for which we have no genealogical records (F25, F55, M61, and F67). We found that FPED (where known) consistently underestimated FROH (one-tailed paired Wilcoxon signed-rank test, P = 7.81 × 10−3; Fig. 3B). We repeated the analysis with ROH of 1 Mb or longer, since these reflect recent rather than ancient inbreeding and can be unambiguously identified with high confidence, and obtained the same result (P = 7.81 × 10−3). Thus, recent severe inbreeding in Isle Royale wolves has sharply increased homozygosity across the genome, beyond the expected values suggested by the known pedigree.

Last, we evaluated FROH in Isle Royale wolves over time to identify patterns consistent with the known history of the population (Fig. 3C). The Isle Royale wolves in our study were born between 1983 and 2009. In the early 1980s, the Isle Royale population crashed from an all-time high of 50 individuals to just 14 within a span of 2 years due to an outbreak of canine parvovirus (40). The population remained small until the mid-1990s, when it began to gradually increase. We found an overall trend of increasing FROH consistent with high levels of inbreeding during this time of stagnant population growth (Fig. 3C). In 1997, a male wolf naturally immigrated to Isle Royale (individual “M” in Fig. 1A) and became a highly successful breeder (24). As shown previously, this event led to an immediate decrease in the mean level of inbreeding and an increase in heterozygosity in the population (24). However, the migrant was so successful that his genome effectively swamped the population within ~2 generations and the level of inbreeding rapidly rose again (24, 25). We found that patterns of FROH in the migrant’s descendants reflect this series of events, with an initial drop in FROH directly after the migrant’s arrival, followed by a rapid increase in premigrant levels by the late 2000s, by which time the population had entered its final period of decline (Fig. 3C). Unfortunately, the migrant individual was never directly sampled; therefore, we were unable to sequence his genome for this study. Nonetheless, we found that overall patterns of ROH in Isle Royale wolf genomes show trends consistent with notable events that occurred during this period of the population’s history.

Phenotypic evidence of inbreeding depression in Isle Royale wolves

Isle Royale wolves have previously been found to suffer from a high incidence of congenital vertebral anomalies, including extra vertebrae, severely asymmetrical vertebrae, and transitional vertebrae (vertebrae exhibiting characteristics of two different types of vertebra, e.g., thoracolumbar, lumbosacral, and sacrococcygeal transitional vertebrae) (8). We examined skeletal remains collected post-mortem from 6 of the 11 wolves sequenced for this study and a litter of 8 pups to identify congenital malformations, specifically in the spine and rib cage. Only one of these specimens (F65) was part of the 2009 study. The incidence of observed anatomical defects was high: Only one individual was free of aberrations (M61), one had a minor transitional variation (F65), and the remaining individuals (F75, M152, M175, and F189) had two to four abnormalities each, including transitional vertebrae, extra vertebrae, and extra ribs (Table 1). One wolf, F75, had three transitional vertebrae and an extra pair of ribs. This wolf was the product of a brother-sister mating and died while giving birth to a litter of pups thought to have been sired by her father, M62 (25). The pedigree-based inbreeding coefficient of this litter, in which all pups showed vertebral changes and all but one had extra ribs, was 0.375 (25). The pups died within hours, and our examination revealed that all eight had extra vertebrae and all but one had extra ribs.

Table 1 Vertebral phenotypes of Isle Royale wolves.

Six individuals were examined for vertebral defects, and their phenotypes are listed. The skeletal remains of the other Isle Royale wolves sequenced in this study (F25, F55, M62, F67, and M141) were not recovered and were therefore unavailable for examination. LSTV, lumbosacral transitional vertebra; SCTV, sacrococcygeal transitional vertebra; TLTV, thoracolumbar transitional vertebra; N/A, not applicable.

View this table:

The fitness effects of defects observed in Isle Royale wolves are uncertain, but lumbosacral transitional vertebrae (LSTV) in dogs are reportedly a predisposing factor of a potentially debilitating condition known as cauda equina syndrome (41). Wolf F65 exhibited a minor aberration that presumably would have no impact on fitness, but two other wolves with LSTV (F75 and M175) showed severe changes with a combination of malformed features. The incidence of LSTV in outbred wolf populations in Finland and historic Scandinavia is reportedly 0 to 1%, compared to 10% in modern inbred Scandinavian wolves, and 33% in Isle Royale wolves before 2007 (7, 8). In addition, three of the six specimens (F75, M152, and M175) exhibited spondylosis deformans, a condition in which the vertebrae have varying degrees of bone spurs (osteophytes) and bony bridges across the vertebral disc space. These lesions are not congenital but are a degenerative change commonly related to age (42) that may also form as a result of traumatic fractures or abnormal vertebral anatomy and alignment (43). Two specimens (F75 and M152) had severe spondylosis that may have affected their spinal mobility. In summary, wolves on Isle Royale suffer from a suite of vertebral changes, including types that can impair mobility and fitness, that are rarely observed in outbred wolves.

The genetic basis of inbreeding depression in Isle Royale wolves

To investigate the genetic basis of inbreeding depression in Isle Royale wolves, we analyzed mutations within protein-coding regions of the genome, which are more likely to directly affect fitness and are also more amenable to functional interpretation. Variants in coding regions were annotated with respect to their predicted impact on the encoded amino acid (e.g., synonymous, nonsynonymous, etc.), and nonsynonymous SNPs were further classified as putatively deleterious or tolerated on the basis of amino acid conservation across taxa (mutations at sites highly conserved over evolutionary time are assumed to be deleterious) (44). We then grouped variants as putatively damaging or benign; here, synonymous and tolerated missense mutations comprise the benign group, whereas deleterious missense mutations and mutations that disrupt splice sites or start or stop codons comprise the damaging group. Studies in humans have revealed that long ROH are enriched for homozygous deleterious variants (45). We did not observe enrichment of deleterious variants within ROH in our dataset (fig. S5). However, the comparison of our results with those obtained in humans is complicated by biological differences (e.g., highly divergent demographic histories of genomes within our study) and technical differences between studies (e.g., possible errors in the classification of variants as damaging or benign in wolves and different methods for identifying and classifying ROH).

To control for background demographic history and gauge specifically how recent inbreeding has affected Isle Royale wolf genomes, we focused on comparing Isle Royale to mainland Minnesota wolves, which have a shared history until the recent founding of the Isle Royale population. Overall, both damaging and benign variants in Isle Royale wolves have shifted from the heterozygous to the homozygous state (Fig. 4, A and B). The proportion of damaging homozygous genotypes was 38.4% higher in Isle Royale compared to the mainland (MWU test, P = 3.23 × 10−4) and 28.4% higher for tolerated homozygotes (MWU test, P = 6.46 × 10−4). However, the total proportion of derived (i.e., nonancestral) alleles per genome was unchanged between the Minnesota and Isle Royale wolves for both damaging and benign mutations (MWU test, P > 0.18; Fig. 4C), consistent with population genetic theory, as inbreeding distorts genotype frequencies rather than allele frequencies. Further, because the founding of the Isle Royale population was very recent [~16 generations ago, assuming 4.2 years per generation; (40)] and its numbers have remained low, an accumulation of new deleterious variants entering the population through mutation or extensive drift of existing weakly deleterious alleles would not be expected. In other words, inbreeding depression in Isle Royale wolves cannot be explained by an increase in the number of derived deleterious variants per genome; rather, it must be accounted for by the increased homozygosity of deleterious variants.

Fig. 4 Genotype and allele frequencies in Isle Royale versus mainland Minnesota wolves.

(A) Inbred Isle Royale wolves contain significantly fewer heterozygotes and (B) significantly more homozygotes than Minnesota wolves, for both damaging and benign SNPs. (C) The total number of derived alleles per individual is unaffected by recent inbreeding. ***P < 0.001. NS, not significant. (D to G) Two-dimensional allele frequency spectra showing the correlation in derived allele frequencies (AF) between outbred North American wolves (Quebec, Yellowstone, and Canadian Arctic excluding Ellesmere Island) and Isle Royale or Minnesota wolves for variants present within Isle Royale or Minnesota wolves. Derived nonsynonymous variants occur at low frequencies due to negative selection in large outbred populations, (D and F) but nonsynonymous variants segregating in Isle Royale wolves occur at much higher frequencies due to drift and high relatedness among individuals (E and G). This effect is apparent for putatively damaging and benign variants. All sites were randomly down-sampled to include exactly five individuals in each group. Color indicates the density of points (see legends). The dashed line represents the diagonal y = x, and the solid line represents the regression line (see table S3). (D) 4387 sites, (E) 3109 sites, (F) 47,784 sites, and (G) 35,979 sites.

Even severely deleterious recessive alleles are expected to be present in standing genetic variation, segregating at low frequencies in large populations where drift is minimal and inbreeding is rare. In contrast, variants that are at high frequency in large populations of outbred wolves are not likely to be strongly deleterious. In the absence of appreciable gene flow with the mainland, we predicted that strongly deleterious recessive alleles carried in the founder genomes could have attained high frequency within the Isle Royale population. We compared the frequencies of segregating variants in Isle Royale wolves and mainland Minnesota wolves to those of outbred North American wolves by constructing two-dimensional allele frequency spectra and performing linear regression to assess correlations between populations (Fig. 4, D to G, and table S3). For this analysis, all groups were randomly down-sampled to five individuals each. We found that variants with low frequency in noninbred North American wolves are also typically at low frequency in Minnesota wolves, consistent with weak drift and efficient selection in the large outbred mainland population (Fig. 4, D and F). In contrast, the frequencies of nonsynonymous variants within Isle Royale wolves were elevated, demonstrating the effects of isolation and high relatedness among individuals (Fig. 4, E and G). Further, derived allele frequencies in Isle Royale wolves were less correlated with allele frequencies in outbred wolves compared to Minnesota wolves, consistent with our prediction that the founder effect and inbreeding allowed damaging variants to drift and attain high frequencies in Isle Royale wolves (P < 2.22 × 10−16; table S3).

We reasoned that these damaging variants might account for the high incidence of vertebral anomalies observed in Isle Royale wolves. We used the following criteria to identify candidate deleterious variants underlying the phenotypes of Isle Royale wolves in our dataset: (i) homozygous and located within ROH in the affected individuals (F65, F75, M152, M175, and F189) but heterozygous or absent in the unaffected individual (M61) and (ii) low frequency (<10%) among other gray wolves (n = 21). Two hundred sixty-four genes containing these mutations were found in at least one affected individual, and two of these genes contained mutations in all five affected individuals: RTTN and SCUBE2. Only RTTN (rotatin) is known to be associated with abnormal phenotypes and notably affects vertebral development (46). The five Isle Royale wolves with vertebral malformations are homozygous for a C-to-T transition that converts a leucine residue to phenylalanine in exon 11 of RTTN. The SIFT (Sorting Intolerant From Tolerant) score for this mutation is 0, indicating very strong conservation at this site and that this mutation is therefore predicted to have a deleterious effect. Two other Isle Royale wolves were also homozygous for this mutation, and three were heterozygous. The mutation was absent from all other individuals except for two Minnesota wolves, which were heterozygotes. Of the 264 candidate genes we identified, 10 are associated with the Human Phenotype Ontology (HPO) term “abnormality of the vertebral column” (HP:0000925) but were not shared across all five affected individuals: ABCC6, CAPN1, ELN, ERCC1, MLXIPL, PSAT1, TCTN2, TERT, ENSCAFG00000001588 (ortholog of SLC52A2), and ENSCAFG00000006532 (ortholog of DCHS1) (fig. S6). The heterogeneity of abnormal phenotypes in our sample of Isle Royale wolves suggests that multiple genes may be involved.

Testing models for the mechanistic basis of inbreeding depression through simulation

Our genomic results suggest that the sharp decline of fitness in Isle Royale wolves is due to increased homozygosity of strongly deleterious recessive alleles rather than a change in the number of deleterious alleles per genome, contrasting with a number of studies showing reduced diversity and an elevated burden of deleterious alleles in historically small or bottlenecked populations [e.g., (17–21)]. However, current genome annotation methods cannot reliably distinguish between strongly and weakly deleterious mutations in nonmodel organisms nor determine whether mutations are additive or recessive. We therefore used simulations to investigate expected patterns of deleterious variation in populations with historically large versus small sizes to infer the likely impact of inbreeding and the mechanism of inbreeding depression in these populations. We performed forward simulations in SLiM (47) under a two-population model incorporating estimates of the long-term effective population sizes (31) of outbred North American (Ne = 17,350) and Tibetan wolves (Ne = 2500) (Fig. 5A) to compare the number of weakly (0 < Nes ≤ 10), moderately (10 < Nes ≤ 100), and strongly (Nes > 100) deleterious mutations. Here, Ne is the ancestral population size before the North American and Tibetan populations diverged (Ne = 45,000). We used the MWU test to evaluate statistical significance in comparisons between the two populations. We conducted one set of simulations in which all mutations were additive (h = 0.5) and one in which all were recessive (h = 0) to explore the effects of dominance. The simulated North American wolf population represents the source population for Isle Royale wolves, which only became isolated fewer than 100 years ago, in contrast to the Tibetan wolf, which has had a small effective population size relative to other gray wolves for more than 10,000 years (31).

Fig. 5 Model and results from simulations of deleterious variation.

(A) Demographic model used to simulate the expected number and age of mutations in a large population (North America, 17,350 individuals) versus a small population (Tibet, 2500 individuals). Both populations split from a large ancestral population (45,000 individuals) 12,500 years before the present (3-year generation time). Population sizes and split time from Fan et al. (31). Model not drawn to scale. (B and C) Results from simulations were grouped according to dominance and selection coefficients. Additive, h = 0.5; recessive, h = 0; strongly deleterious, Nes > 100; moderately deleterious, 100 ≥ Nes > 10; weakly deleterious, 10 ≥ Nes > 0; neutral, Nes = 0. (B) The average number of homozygous derived genotypes and total derived alleles per individual (ind.). (C) The average ages of segregating mutations in each population.

Our simulations confirmed that the number and the homozygosity of deleterious alleles per individual are expected to differ significantly between the larger North American and the smaller Tibetan populations (Fig. 5B). Except in the case of strongly deleterious mutations, homozygosity was consistently higher in the Tibetan population (P < 1.31 × 10−32). Regarding the number of deleterious alleles per genome, however, different patterns emerged depending on whether mutations were additive or recessive. In simulations with additive mutations, the overall number of deleterious alleles per individual was slightly higher in the smaller Tibetan population relative to the larger North American population (+1.87%; P = 1.84 × 10−2) due to the accumulation of moderately deleterious alleles (+34%; P = 2.49 × 10−11). There was no significant difference in the numbers of strongly or weakly deleterious additive alleles between the two simulated populations. Thus, in an additive model, selection against strongly deleterious alleles is not hindered by drift, but long-term small population size is predicted to result in an overall increase in the burden of deleterious variants due to the accumulation of moderately deleterious alleles.

In simulations with recessive mutations, however, a very different pattern emerged. Here, the total number of deleterious alleles per simulated genome was higher in the larger North American population (+3.58%; P = 1.42 × 10−6) due to sharply elevated numbers of strongly (+59.2%; P = 2.56 × 10−34) and moderately (+23.4%; P = 3.04 × 10−24) deleterious recessive alleles relative to the Tibetan population. Increased rates of homozygosity in the simulated Tibetan population caused by its smaller effective size allow recessive deleterious alleles to be removed more efficiently than in the larger population. Furthermore, the simulations showed that the mean age of segregating strongly deleterious recessive mutations was 2.22-fold greater in the North American relative to the Tibetan population (P = 2.56 × 10−34; Fig. 5C), indicating that these mutations can persist over longer time periods while remaining hidden from selection as heterozygotes in the larger population. By contrast, segregating mutations in all other categories tended to be younger in the North American population, as a consequence of more new mutations arising each generation in the larger population. On the basis of these simulations, inbreeding would be expected to produce more strongly deleterious homozygous mutations in individuals drawn from a historically large population, such as the mainland source population for Isle Royale wolves, compared to those drawn from a historically smaller population.


The persistence of wolves on Isle Royale was once used to support the claim that a very small population in isolation may persist, and even thrive, without succumbing to genetic deterioration (48). During the past few decades, however, Isle Royale wolves have experienced a precipitous decline following generations of inbreeding and physical degeneration (8, 25). Although inbreeding depression was not the sole determinant, it undoubtedly played a role in the collapse of the population, along with stochastic demographic and environmental events, such as periodic disease outbreaks, severe winters, and the drowning of three wolves in a flooded abandoned mine shaft in 2011 (25). The genomes of Isle Royale wolves bear the hallmarks of their extreme demographic history characterized by extensive ROH, in some cases, spanning whole chromosomes, leading to a marked increase in the homozygosity of deleterious variants.

We found no difference in the total number of derived putatively deleterious alleles between Isle Royale and mainland Minnesota wolf genomes. Instead, reduced fitness in Isle Royale wolves can be explained by increased homozygosity of recessive mutations, rather than an increase in the additive genetic load. Alternatively, the additive genetic load may be higher in Isle Royale wolves, but we are unable to detect this increased load due to the difficulties of determining which genetic variants negatively affect fitness. Currently, predicting and calculating the burden of strongly deleterious recessive alleles within a genome is challenging, particularly in nonmodel species. Here, we focused on SNPs within protein-coding regions and assumed that amino acid changing mutations at highly conserved sites and mutations predicted to markedly affect protein structure were deleterious, while others were relatively benign. These assumptions may not always hold, and this framework also does not account for other potentially deleterious changes (e.g., structural variants, mutations in regulatory regions, etc.) or for beneficial mutations. Nonetheless, our findings suggest that populations with a similar total number of putatively deleterious alleles may differ markedly in terms of fitness. With genomic studies of small populations and deleterious variation in small populations accumulating [e.g., (1721, 49, 50)], we propose that additional metrics should be used to quantify the load in populations (51).

Further, our simulated results affirm that patterns of deleterious variation within a genome (i.e., total number of alleles, dominance, deleteriousness, and zygosity) are strongly influenced by demographic history. We found that individuals derived from populations with historically small sizes have a higher burden of deleterious alleles overall, consistent with population genetic theory and the results of recent empirical studies of genomes from populations with long-term reduced effective size [e.g., (1721)]. However, we also found that individuals from these smaller populations may carry fewer strongly deleterious recessive alleles in the heterozygous state, reducing the risk of severe fitness declines due to inbreeding. The Ethiopian and Tibetan wolves may represent examples of this phenomenon; both exist in small isolated populations (30, 31, 36, 37) but are not known to suffer from inbreeding depression. Although the Tibetan wolf has not been well studied, the Ethiopian wolf has been studied for decades [e.g., (36, 37, 52, 53)], and notable signs of inbreeding depression, such as those in Isle Royale wolves, have not been reported. Further studies are needed to test this possibility, but purging of strongly deleterious recessive alleles may be one reason why some small isolated populations have existed at small sizes for many generations (20, 21, 49, 50). In contrast, individuals derived from historically large populations, such as the mainland source population for Isle Royale wolves, carry a greater burden of strongly deleterious recessive mutations in the heterozygous state relative to individuals from historically smaller populations. Strongly deleterious recessive alleles carried as heterozygotes within the genomes of founders or immigrants quickly become homozygous through inbreeding in small isolated populations, resulting in inbreeding depression. Similar phenomena have been noted in maize, which were domesticated from historically large populations. Here, the initial inbred lines exhibited severe inbreeding depression and reduced yield, creating the need for hybrid lines to mask the effects of recessive deleterious alleles (54).

The collapse of the Isle Royale wolf population occurred despite a reported genetic rescue and evidence of earlier sporadic migration events from the mainland. Previous genetic analysis revealed that undetected migration from the mainland may have occurred during winters that were cold enough for an ice bridge to form (24, 25). However, warmer winters over the past several decades have resulted in a marked reduction in the formation of ice bridges, a trend that is likely to continue in a warming climate (25). The reported genetic rescue of the Isle Royale wolf population by a single migrant from the mainland in 1997 also appears to have been short-lived (25, 55). This male wolf was such a successful breeder that within 2.5 generations, every individual in the Isle Royale population was related to him, leading to intense inbreeding among his descendants (Fig. 1A) (24, 25). A similar episode occurred in the inbred Scandinavian wolf population following the arrival of an immigrant wolf in 1991, but the population subsequently entered a period of growth in the wake of additional immigration (13). A higher rate of gene flow between the mainland and the island following the colonization of Isle Royale may have mitigated or postponed the onset of inbreeding depression, but the effective rate of naturally occurring gene flow was insufficient to sustain the population.

Consequently, periodic human-assisted gene flow to offset inbreeding may be necessary for the long-term persistence of a reintroduced population on Isle Royale. Alternatively, our simulation results suggest that selecting founders from a historically isolated population may reduce the risk of future inbreeding depression, since these founder individuals would be expected to carry a reduced load of recessive deleterious alleles in the heterozygous state. This alternative strategy has drawbacks, however; for one, its success requires the absence of gene flow with large populations nearby that could introduce strongly deleterious recessive alleles back into the smaller population. Given the intermittent migration that may occur with the mainland, sourcing from a small population may not be advisable for the Isle Royale reintroduction. Moreover, individuals from historically smaller populations are predicted to have a higher burden of weakly deleterious additive alleles, which are not purged effectively in smaller populations, and they may also have a reduced capacity for adapting to new environmental conditions in the reintroduction site due to lower genetic variability. Nonetheless, the idea that founders from historically isolated populations should be selected to mitigate the risk of future inbreeding depression stands in contrast with the conventional wisdom of selecting individuals with the most diversity from nearby populations. However, this alternative should be considered as a potentially more successful strategy for preserving species consisting of small and isolated populations, which is the inevitable future for many species.

Last, life history traits must be considered when determining the best course of action to reduce the risk of inbreeding. For example, wolves typically avoid mating with close relatives through the exchange of individuals between different packs, which are usually familial units that consist of a breeding pair and its offspring (56). Previously, Isle Royale sustained three to four wolf packs; hence, the estimate that the long-term effective population size of wolves on Isle Royale was a mere 3.8 individuals (40). In simulations with models incorporating ecological and demographic stochasticity, the mean time to extinction for social organisms, specifically wolves, was found to be strongly tied to the number of social groups rather than the number of individuals (57). Thus, the complex social behavior of wolves further reduces their ratio of effective population size relative to census population size and likely contributed to the extreme inbreeding that occurred on Isle Royale. A minimum number of wolves to sustain multiple packs, maximizing the number of breeding pairs within an isolated population such as Isle Royale, is therefore essential for population persistence.

The rapid deterioration of fitness in Isle Royale wolves within just a few generations demonstrates the substantial risks of inbreeding in very small isolated populations. In the absence of recurring immigration, whether managed or not, the fate of a restored population on Isle Royale is grim. Nevertheless, wolf predation is an important top-down influence on Isle Royale, and its absence threatens the stability of the island ecosystem (23). Active management of the population to reduce inbreeding may improve its long-term viability. We also recommend genomic screening of all wolves reintroduced to Isle Royale and their eventual descendants for the purposes of cataloging and monitoring ROH and deleterious alleles over time. Although Isle Royale represents an extreme scenario, evidence of inbreeding (and, in some cases, inbreeding depression) has now been documented in multiple isolated wolf populations in Europe and North America [e.g., (68, 1113, 28, 29)]. The demise of the iconic Isle Royale wolf population provides lessons for its future restoration and guidance for the management of other species or populations to minimize the consequences of inbreeding.


Experimental design

Samples and sequencing. DNA from Isle Royale wolves was extracted from blood samples archived at Michigan Technological University. Before this study, the Isle Royale population was not actively managed, and human interventions (e.g., trapping and fitting radio collars) were kept to a minimum. Thus, many wolf samples have been collected opportunistically over the years from urine, scat, and carcasses, and complete lifetime survival and reproduction data are generally not available. Of the DNA samples present in the archives, those with the highest concentration and quality were selected for sequencing to ensure high-quality sequence data and high coverage. DNA from Minnesota and Canadian Arctic wolves was extracted from blood and tissue samples from the archive of R.K.W. that were used in previous studies (34, 58). Although our sample of Minnesota wolves may not be directly derived from the ancestral population for Isle Royale wolves, a previous analysis of 3 Isle Royale and 11 Minnesota wolves at >44,000 SNP loci suggests that the two populations have very high genetic similarity (34). Whole-genome sequencing was performed on an Illumina HiSeq 4000 at the Vincent J. Coates Genomics Sequencing Laboratory at University of California, Berkeley. Previously sequenced genomes with high coverage were downloaded from the National Center for Biotechnology Information (NCBI) Sequence Read Archive (see table S1).

Statistical analyses

Read processing and alignment. A pipeline adapted from the Genome Analysis Toolkit (GATK) (59) Best Practices Guide was used to process raw reads before genotype calling. Briefly, paired-end raw sequence reads, 150 bases in length, were aligned to the domestic dog reference genome, canFam3.1, using BWA-MEM (60) before the removal of polymerase chain reaction duplicates and low-quality reads. Lacking a database of known variants, the bootstrapping method of base quality score recalibration as recommended by GATK was performed by calling raw genotypes with GATK UnifiedGenotyper (minimum base quality Phred score of 20) and using these variants as input for recalibration with BaseRecalibrator. This process was repeated three times to reach convergence between reported and empirical quality scores. The full pipeline used for read processing and alignment, with example commands, is available at

Genotype calling and filtering

Joint genotype calling at all sites across the reference genome (including invariant positions) was performed with GATK HaplotypeCaller. Genotypes were filtered for quality and depth, leaving only high-quality biallelic SNPs and monomorphic sites. Only genotypes with at least six supporting reads and high quality (minimum Phred score of 20) were included. A per-individual excess depth filter, set at the 99th percentile of depth for each sample (to control for differences in coverage between samples), was also used. Each site was then filtered on the following criteria: Variants failing the recommended GATK hard filters were excluded (QD < 2.0 || FS > 60.0 || MQ < 40.0 || MQRankSum < −12.5 || ReadPosRankSum < −8.0 || SOR > 3.0), as well as sites with excess depth (>99th percentile for total depth across all samples), low Phred score (QUAL of <30), more than 20% missing data, excess heterozygosity (>50% of individuals heterozygous), or sites found within repeat regions and CpG islands [coordinates from (19)]. The pipeline used for genotype calling and filtering, with example commands, is available at

Variant annotation and polarization of derived versus ancestral alleles

We annotated variant sites using the Ensembl Variant Effect Predictor (VEP) [version 87; (61)] with SIFT enabled (44). SIFT determines whether a nonsynonymous mutation is likely to be damaging or benign on the basis of phylogenetic constraint. The SIFT scores derive from a precomputed database included with VEP, which includes scores for all possible amino acid changes in the dog reference genome. According to the VEP manual, the scores were generated with SIFT (version 5.2.2) from alignments of homologous sequences in the UniRef90 (release 2014_11) protein database, according to the protocol recommended by the SIFT developers. SIFT scores were therefore not calculated from the variation present within our dataset. We grouped protein-coding variants into “damaging” and “benign” classes as follows: Damaging variants included “deleterious” nonsynonymous variants (SIFT score of <0.05) and variants that disrupted splice sites, start codons, or stop codons; benign variants included “tolerated” nonsynonymous variants (SIFT score of ≥0.05) and synonymous mutations. Alleles were polarized as derived or ancestral with respect to the gray fox and African golden wolf genomes. Specifically, when both outgroup genomes contained an allele that did not match the reference, the polarity was reversed (i.e. the nonreference allele was assumed to be ancestral) with probability equal to the frequency of the nonreference allele in the outgroup genomes. For example, the polarity was changed with a probability of 50% if the nonreference allele was present at a frequency of 0.5 in the gray fox and African golden wolf genomes.

Phylogenetic analysis and cladogram construction

A cladogram representing the relationships between 20 genomes in this study was constructed using SNPhylo (62). Where multiple individuals were available from a single population, one individual was chosen at random for inclusion in the tree (Minnesota, RKW119; Isle Royale, CL141; and Yellowstone, 569F). SNPs were pruned for linkage disequilibrium (threshold of 0.25), and SNPs with a minor allele frequency of <0.1 or with missingness above 10% were excluded. The tree was constructed with the 38,408 remaining SNPs. One thousand bootstrap replicates were performed.

Morphological analysis

Skeletons from 6 of the 11 Isle Royale wolves included in this study were retrieved from storage at Michigan Technological University, photographed, and assessed for morphological anomalies as described previously (7, 8). Skeletons were obtained from carcasses collected post-mortem. The vertebral column and rib variations were evaluated. In some cases, specimens were missing a few vertebrae or ribs; the list of examined wolves is noted in Table 1. A litter of eight newborn pups from F75 was also examined through radiographic analysis.

Calculation of genome-wide heterozygosity

In this study, we calculated heterozygosity as the number of heterozygous genotypes divided by the total number of called genotypes within a single individual. For each individual, we calculated heterozygosity for the entire autosomal genome and in nonoverlapping 1-Mb windows across the autosomes. Windows where more than 80% of sites failed filters or were missing were excluded.

Identification and analysis of ROH

ROH were identified using VCFtools (63). ROH spanning regions smaller than 100 kb were ignored. Highly similar patterns of ROH were also inferred using PLINK v1.90b6.7 (64) with the following options: --homozyg-window-snp 200 --homozyg-window-het 3 --homozyg-window-missing 50. Overall, we found that PLINK identified a greater number of short ROH (100 kb to 1 Mb) relative to VCFtools and also inferred fewer medium ROH (1 to 10 Mb) in the Ethiopian wolf genome (fig. S7). We used ROH identified by VCFtools in further analyses.

The proportion of the genome within ROH (FROH) was calculated as the total length of ROH of 100 kb or longer within an individual divided by the length of the genome (2,203,764,842 nucleotides). The X chromosome was included in the calculation for females only (total genome length = 2,327,633,984 nucleotides). To test for enrichment of homozygous deleterious variants within ROH, we followed the method of Szpiech et al. (45). The fraction of the genome within ROH and the fraction of damaging and benign homozygous genotypes inside ROH were calculated in each individual. Following Equation 10 of Szpiech et al. (45), we fit linear models to test whether variant impact (benign versus deleterious) and FROH were significant predictors (β2 and β3, respectively) for the proportion of nonreference homozygotes within ROH.

Identification of candidate genes underlying Isle Royale phenotypes

We identified candidate mutations underlying the abnormal phenotypes observed within our sample of Isle Royale wolves as those satisfying the following criteria: (i) Mutations had to be classified as damaging (missense mutations classified as deleterious by SIFT and mutations disrupting splice sites or start/stop codons) and passing all quality control filters described above. (ii) Mutations had to be homozygous and contained within ROH of 100 kb or more in the affected wolves (F65, F75, M152, M175, and F189) but not homozygous for the derived allele in the unaffected wolf (M61). (iii) Only mutations with low frequency (<10%) in 21 other gray wolves [wolves from Minnesota (n = 6), Canadian Arctic (n = 4), Yellowstone (n = 3), Quebec, Mexico, Portugal, Spain, Italy, Iran, Tibet, and Xinjiang] were considered. Alleles were polarized as ancestral or derived with respect to the gray fox and African wolf outgroups as described above. Genes containing mutations satisfying all criteria were extracted, and HPO terms associated with the genes were obtained with gProfileR [r1741_e90_eg37 release; (65)].

Simulations of neutral and deleterious variation

We performed forward genetic simulations in SLiM [version 2.4.2; (47)] under a divergence model with parameters estimated by Fan et al. (31). Each simulated individual consisted of a diploid 1-Mb genome, with a simple architecture of 1000 “genes” carried on 38 chromosomes proportional to chromosome lengths in the dog genome. Each gene consisted of a contiguous 1-kb sequence that accumulated mutations at a rate of 1 × 10−8 per site per generation. Selection coefficients for deleterious mutations were drawn from the distribution of fitness effects inferred from a large sample of humans by Kim et al. (66). Seventy percent of mutations were deleterious, and the remaining 30% were neutral (s = 0). Each simulation began with a burn-in period of 450,000 (10 × Ne) generations to allow the ancestral population to reach equilibrium. Recombination was permitted at single base positions between each gene at a rate of 1 × 10−3 per site per generation to simulate the effective rate of crossing over that would occur in 100-kb noncoding regions between each gene, assuming a recombination rate of 1 × 10−8 per site per generation. At the end of each simulation, the average number of alleles per individual and the average age of segregating mutations were calculated for weakly (0 < Nes ≤ 10), moderately (10 < Nes ≤ 100), and strongly (Nes > 100) deleterious and neutral (Nes = 0) mutations. We performed 100 replicates in which mutations were additive (h = 0.5) and 100 replicates in which mutations were completely recessive (h = 0.0) to examine the effects of dominance. Simulation scripts are available at


Supplementary material for this article is available at

Fig. S1. Cladogram of genome sequences in this study.

Fig. S2. Genome-wide distributions of heterozygosity in all individuals.

Fig. S3. ROH sharing among Isle Royale wolves across the autosomal genome.

Fig. S4. Number of ROH per individual for different ROH length categories.

Fig. S5. Proportion of derived homozygotes in ROH as a function of the proportion of the genome within ROH.

Fig. S6. Candidate genes underlying Isle Royale phenotypes associated with HPO terms related to skeletal anatomy.

Fig. S7. Comparison of ROH identified in VCFtools and PLINK.

Table S1. Sample information for sequences included in this study.

Table S2. Notes and information from the literature about the demographic history of populations included in this study.

Table S3. Parameters of linear regression models for two-dimensional allele frequency spectra.

References (6776)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank P. W. Hedrick for helpful discussions in preparation of this manuscript. We acknowledge N. C. Nelson from Michigan State University for providing radiographs of the wolf pups. Funding: This work used the Vincent J. Coates Genomics Sequencing Laboratory at UC Berkeley, supported by NIH S10 OD018174 Instrumentation Grant. This work was funded by NIH (grant R35GM119856 to K.E.L.), NSF (DEB-1453041 to J.A.V.), Isle Royale National Park (CESU task agreement no. P16AC00004, under Master Cooperative Agreement no. P12AC31164), the Robbins Chair in Sustainable Management of the Environment to R.O.P. at Michigan Technological University, and McIntyre-Stennis Grant (USDA-NIFA-1014575). Author contributions: J.A.R., K.E.L., R.K.W., and R.O.P. conceived and designed the study. L.M.V., J.A.V., and R.O.P. performed the sample acquisition. J.R. carried out morphological analyses. J.A.R. carried out genomic analyses and simulations and wrote the manuscript. All authors contributed to manuscript revision and approved the final version. R.K.W. and K.E.L. jointly supervised this work. Competing interests: The authors declare that they have no competing interests. Data and materials availability: Newly generated genome sequence data were deposited in the NCBI Sequence Read Archive under BioProject PRJNA512209. All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article