Research Article

Parallel adaptation of rabbit populations to myxoma virus

See allHide authors and affiliations

Science  22 Mar 2019:
Vol. 363, Issue 6433, pp. 1319-1326
DOI: 10.1126/science.aau7285

Locating myxomatosis resistance

Myxomatosis is a viral infection that was deliberately introduced from American cottontail rabbits into European rabbit populations to control their population. Over the past 60 years or so, similar resistance variants have emerged in parallel in the United Kingdom, France, and Australia. Alves et al. discovered that the basis for this resistance is polygenic, with selection converging on several host immunity and proviral alleles (see the Perspective by Miller and Metcalf). Interestingly, it now seems that the virus is counterevolving immune suppressive traits.

Science, this issue p. 1319; see also p. 1277


In the 1950s the myxoma virus was released into European rabbit populations in Australia and Europe, decimating populations and resulting in the rapid evolution of resistance. We investigated the genetic basis of resistance by comparing the exomes of rabbits collected before and after the pandemic. We found a strong pattern of parallel evolution, with selection on standing genetic variation favoring the same alleles in Australia, France, and the United Kingdom. Many of these changes occurred in immunity-related genes, supporting a polygenic basis of resistance. We experimentally validated the role of several genes in viral replication and showed that selection acting on an interferon protein has increased the protein’s antiviral effect.

The emergence of new infectious diseases can result in intense selective pressures on populations and cause rapid evolutionary change in both host and parasite. Whereas pathogens must adapt to a new ecology and cellular environment, hosts can rapidly evolve resistance in a continuous arms race. One of the most emblematic examples of this coevolutionary process arose when the wild European rabbit (Oryctolagus cuniculus) populations in Australia and Europe were exposed for the first time to the myxoma virus (MYXV) (genus Leporipoxvirus, family Poxviridae). MYXV circulates naturally in American cottontail rabbits (Sylvilagus spp.), where it causes benign skin tumors, but in European rabbits it causes the highly lethal systemic disease myxomatosis (1).

Rabbits were initially introduced into Australia by European settlers, resulting in extensive ecological and economic damage (2). In an attempt to control the rabbit populations, MYXV was released in 1950 in Australia, after which it spread rapidly across the country, causing massive population reductions (1). In 1952 it was released in France, and in 1953 it reached the United Kingdom, with similar outcomes (2). In a series of classic experiments, the evolution of the virus and rabbits was tracked, and in all three countries, substantial declines in case fatality rates in wild rabbits were observed, due both to the evolution of less virulent viral phenotypes and to increased resistance in rabbit populations (35). Considered “one of the greatest natural experiments in evolution,” these observations ultimately became a textbook example of host-parasite coevolution (2).

Sixty-nine years have passed since MYXV was first released in Australia, and today the virus continues to evolve in an ongoing arms race against the rabbit immune system (6, 7). Despite much research on the genetics of MYXV, little is known about the genetic basis of resistance to myxomatosis. The intense selective pressure exerted by the disease in rabbits and the parallel phenotypic evolution across multiple populations provide an exceptional framework to study the evolution of resistance to a highly lethal pandemic in natural host populations. To understand the genetic basis of these changes, we focused on three countries where genetic resistance independently emerged: Australia, France, and the United Kingdom. For each country, we compared modern rabbits with historical specimens, obtained from museums, that were collected before or soon after the release of the virus (between 1856 and 1956) (Fig. 1, fig. S1, and data file S1).

Fig. 1 Rabbit origins and sampling locations.

Historical (circles) and modern (triangles) sampling locations. Dates in red inside the maps are the dates of the first myxomatosis outbreaks in the respective countries. Orange dashed arrows and dates reflect historical and archaeological records of the colonization of European rabbits from France to the United Kingdom and Australia.

Colonization route of rabbits

To obtain genome-wide polymorphism data, we used oligonucleotide probes to capture 32.10 Mb (1.17%) of the rabbit genome, including the exome (19,293 genes), the mitochondrial genome, and three genomic regions that contain the major histocompatibility complex (MHC) region, encompassing 1.75 Mb. We sequenced a total of 152 rabbits from Australia (historical, n = 17 rabbits; modern, n = 26 rabbits), France (historical, n = 29; modern, n = 26), and the United Kingdom (historical, n = 29; modern, n = 25). The mean sequence coverage of on-target reads per individual after filtering was 33×. After filtering, the number of biallelic single-nucleotide polymorphisms (SNPs) was 757,333.

Historical records support the introduction of rabbits to the island of Great Britain from continental Europe by the 13th century (8), and most Australian rabbits are thought to be derived from an introduction in 1859 from the United Kingdom (9) (Fig. 1). Our genetic data reflect these historical records. Both principal components analysis (PCA) (Fig. 2A) and structure analyses (Fig. 2B) (K = 3, where K is the number of ancestral populations) reveal three clusters composed of individuals from the same country. The sequential colonization from France to the United Kingdom and then Australia is reflected in the levels of genetic differentiation, which are greatest between France and Australia and lowest between the United Kingdom and Australia (table S1). This pattern is repeated with the Australian populations clustering together with the U.K. populations in the structure analysis (Fig. 2B) (K = 2). Population bottlenecks during colonization can reduce genetic diversity and increase linkage disequilibrium (LD). Both aspects are reflected in the successive increases in LD (Fig. 2C) and decreases in heterozygosity (Fig. 2D) as rabbits moved from continental Europe to the United Kingdom and then Australia.

Fig. 2 Genetic structure and diversity in historical and modern populations of France, the United Kingdom, and Australia.

(A) PCA. Ellipses show 95% confidence intervals. (B) Ancestry fractions inferred with Ohana structure analysis for K = 2 (top) and K = 3 (bottom). Each bar shows the inferred ancestry fraction for an individual. Black lines between bars separate populations. Labels above bars identify the country, and labels below bars identify the temporal set and sample size. Individuals are ordered geographically within each population. (C) Decay of LD for each population. Dots represent the mean pairwise r2 values (squared correlation coefficients in the allelic state) between pairs of SNPs in nonoverlapping 500-bp windows. Colors represent different populations. FR, France; AU, Australia. (D) Expected heterozygosity for each population. The bars are the means across chromosome arms, and error bars are 95% bootstrap confidence intervals from resampling of chromosome arms. Colors represent different populations. (E) The correlation between the frequencies of the alternative allele in historical and modern populations for France, the United Kingdom, and Australia. Colors reflect the relative density of points according to the scale in the bottom right of each plot, from darker (more density) to lighter (less density). R2, coefficient of determination.

Genetic variation in historical and modern populations

The genome-wide allele frequencies in our modern and historical samples are similar, an important consideration in the detection of changes in allele frequency due to natural selection. Although changes in allele frequency are an expected signal of selection, artifactual variation generated by DNA degradation can generate similar signals between modern and historical DNA. To mitigate this, we sequenced samples to a high coverage, corrected for the effects of DNA damage patterns, and used a stringent set of filters at read and variant levels. Genome-wide differences in allele frequency between our modern and historical samples may also arise because of population substructure. To minimize this effect, we sampled modern rabbits from locations near where the historical specimens had been sampled (Fig. 1 and fig. S1).

We consistently found that historical and modern populations from the same country exhibit similar patterns of genetic structure and diversity. In the PCA, the 95% confidence ellipses for historical and modern samples from the same country are interspersed (Fig. 2A). In the structure analysis, the patterns again reflect geography rather than the time of sampling (Fig. 2B) (K = 3), and increasing the number of ancestral populations (K) reveals finer population substructure instead of a split between the different time points (fig. S2). More generally, across all SNPs, the allele frequencies of historical and modern populations are highly correlated in the three countries (Fig. 2E).

The collapse of populations because of myxomatosis has not increased levels of LD or decreased genetic diversity, as both these parameters are similar between historical and modern populations (Fig. 2, C and D). This is to be expected, as the effective population size (Ne) of rabbits is ~106 (10), which is large for a vertebrate species. Therefore, even if an ~99% reduction in population size [as seen in some populations immediately after the first epidemic (6)] had been sustained for the subsequent 65 generations, the expected reduction in heterozygosity would have been only 0.35% [per generation reduction in heterozygosity = 1 − 1/(2Ne)]. Similarly, bottlenecks of this size are expected to have little effect on LD (11). Together, these results demonstrate that historical and modern samples are drawn from genetically similar populations. This strong correlation in allele frequencies between the two time points within each country is conducive to the identification of unusual shifts in allele frequency resulting from natural selection.

Parallel changes in allele frequency

With a 99.8% case fatality rate for the originally released strain of MYXV, the myxomatosis pandemic imposed intense selection on rabbits, resulting in rapid and parallel evolution of increased resistance in exposed populations (6). To investigate whether parallel genetic changes occurred in Australia, France, and the United Kingdom, we calculated the fixation index (FST) between the historical and modern samples for each country and identified the 1000 SNPs that show the highest FST. By intersecting these three lists, we found that more alleles have changed in frequency across two or three countries than expected by chance (Fig. 3A and data file S2). Furthermore, considering the 92 SNPs in the intersections of Fig. 3A, in 87 cases it was the same allele that had increased in frequency in the countries involved. The SNPs that are among the top 1000 most differentiated in any two populations tend to show elevated FST in the third population (Fig. 3B). These results demonstrate the occurrence of parallel changes in allele frequency across the three populations.

Fig. 3 Parallel changes in allele frequency across three countries.

(A) Venn diagram showing the overlap of the 1000 SNPs with the greatest changes in allele frequency between modern and historical samples (FST) in France, the United Kingdom, and Australia. Numbers in black are the observed numbers of SNPs, and numbers in red show the expected overlap after 1000 random permutations of modern and historical samples within each country. (B) Scaled histogram of the FST values in the three countries. Bars with dark colors reflect SNPs that are in the top 1000 in both of the other two countries. Bars with light colors represent SNPs that are in the top 1000 for only one of the other countries. Gray bars represent all the remaining SNPs. (C) Genome-wide selection scan on the basis of allele frequency changes after the introduction of myxomatosis (supplementary methods and eq. S5) (the strength of selection in each population is allowed to vary independently). (D) Selection scan testing whether selection has acted in all three populations (positive values) or just one population (negative values) (supplementary methods and eq. S6). [(C) and (D)] The y axis shows the likelihood ratio statistic of each model. The orange dotted lines show the genome-wide 95% significance threshold from permuting sample collection dates within each country. Shaded boxes show SNPs located in the highlighted genes. Different shades of blue represent chromosomes. (E) Mean of the posterior distribution of the derived allele frequency as a function of time for the IFN-α21A and FCRL3 loci from the Bayesian selection analysis (additive model). Ninety-five percent credible intervals are shaded. Triangles across the bottom represent years of samples (only samples post-1940 are shown). Dotted lines show the dates of the first reports of MYXV and RHDV. A list of the top 1000 SNPs for all figures is available in data files S2, S3, and S4.

To locate the putative targets of selection in the genome, we searched for SNPs that showed large changes in allele frequency since the release of the MYXV. We accounted for the effects of population structure by scanning for outliers where the difference in allele frequency between modern and historical populations was larger than that expected given the covariance matrix describing the joint allelic frequencies among populations inferred from the structure analysis (Fig. 2B) (K =3 ). In each population, we assumed that selection started in the year the virus was introduced. For each SNP, we then compared the likelihood of a model where the change in allele frequency followed that predicted from the genome-wide amount of genetic drift across all historical and modern samples with the likelihood of a model that allowed additional changes in allele frequency through time (fig. S3). Combining data across the three populations and allowing the strength of selection to vary independently in the three countries, we identified 193 SNPs in 98 genes and seven intergenic regions that showed significant changes in allele frequency (genome-wide significance, P < 0.05) (Fig. 3C and data file S3). The results were similar when we assumed the same strength of selection in all countries (fig. S4 and data file S3).

To explicitly test for parallel evolution and identify variants that significantly changed in frequency in all three populations, we compared the likelihood of a model that assumed equal selection across the three populations with a model where we allowed the alleles to change in frequency as a result of selection in one population only. We identified 94 SNPs that had a significant signature of parallel selection (genome-wide significance, P < 0.05) (Fig. 3D, positive values).

It is common to find considerable standing genetic variation in susceptibility to infection within populations (12, 13), which may allow populations to rapidly evolve resistance when they encounter new pathogens (14). This can be explicitly tested with our experimental design that incorporated historical samples. We found that in the large majority of cases, the SNPs under selection were also present in at least one of our historical populations (93% and 84% of the SNPs with a genome-wide P value of <0.05 in Fig. 3, C and D, respectively). Therefore, parallel evolution has occurred because natural selection has been acting on standing genetic variation that was present >800 years ago in continental Europe and was carried with rabbits to the United Kingdom and Australia. The existence of this variation likely explains the rapid development of resistance to myxomatosis observed almost immediately after the first outbreaks.

Population-specific evolution

Despite the common selection pressure imposed by myxomatosis, the populations of France, the United Kingdom, and Australia will have experienced their distinct selection pressures, as resistance was evolving in different ecological and genetic backgrounds. We can identify SNPs that have experienced population-specific changes in allele frequency from our previous analysis when the model of selection in one population is preferred over the model of selection in all three populations (negative values in Fig. 3D).

To quantify the proportion of SNPs under parallel or population-specific selection, we used a Bayesian approach. First, we analyzed data from each population independently to identify variants under selection (fig. S5 and data file S4) (genome-wide P < 0.05). This minimizes the inherent bias toward detecting SNPs under parallel selection that occurs when the all populations are combined (Fig. 3, C and D). We then returned to the combined dataset, and for each of the SNPs that were significant in the single population analysis, we calculated a Bayes factor to compare models of population-specific versus parallel selection (fig. S5). For each gene, we retained only the most significant SNP (n = 43 SNPs) (data file S5). We found evidence that 20 SNPs were under population-specific and 20 were under parallel selection (for three SNPs, we could not distinguish the models), implying that a large component of the recent selection in the three populations has occurred on a common set of variants.

Because of population bottlenecks as rabbits colonized new areas (Fig. 2, C and D), alleles selected in one population may be rare or missing in other populations. This means that population-specific adaptation could result not only from differences in selection pressures but also from differences in the variants available for selection to act on. To test this, we examined the frequency of the 20 alleles under population-specific selection in our historical samples. There were no consistent differences in the ancestral allele frequencies between the populations where we detected the effects of selection and the populations where we did not (table S2). Therefore, we can conclude that population-specific changes in allele frequency result from population-specific selection pressures, perhaps due to differences in ecology, genetic background, or independent paths of coevolution with MYXV (15).

Changes in the strength of selection

The proportion of rabbits killed by myxomatosis has fallen since the 1950s because of declines in the virulence of the virus and increases in resistance (6). Therefore, the strength of selection on variants that confer resistance to MYXV is expected to have declined. However, in 1984 a new lethal viral pathogen was identified in rabbits, the rabbit hemorrhagic disease virus (RHDV) (genus Lagovirus, family Caliciviridae) (16). RHDV, which has a case fatality rate similar to that of MYXV, was first found in domestic rabbits in China, from which it spread to continental Europe in 1986, the United Kingdom in 1992, and Australia in 1995 (2, 17, 18). Like myxomatosis, RHDV infection caused high mortality across the two continents, which may have contributed to the observed changes in allele frequency.

To evaluate the role of RHDV in our selection signatures and understand how the strength of selection has changed through time, we obtained 70 rabbit samples that were collected before or soon after RHDV appeared, between 1985 and 1996, in the United Kingdom and Australia (Fig. 3E and data file S1). We selected four SNPs that had shown significant changes in allele frequency since the 1950s and that were located in or near genes with known immune functions (CD96, FCRL3, IFN-α21A, and MHC class I) and genotyped them in these samples by Sanger sequencing. Combining these genotypes with data from our exome sequences, we used a Bayesian approach (19) to estimate the selection coefficient acting on these SNPs. We allowed two phases of selection: the first from the time when MYXV was released to the appearance of RHDV, and the second from the appearance of RHDV to the present day (Fig. 3E) (the date when the strength of selection changes is fixed in the model).

We found that the strength of selection on FCRL3 and IFN-α21A has declined through time in both the United Kingdom and Australia (Fig. 3E and data file S6). This is consistent with selection being driven by MYXV, whose case fatality rate has decreased since its release, and does not support a role for RHDV. The data suggest that these variants have been negatively selected in recent decades (Fig. 3E), which is predicted if they have negative pleiotropic costs on other fitness components. For CD96 and the MHC, there was no significant change in the strength of selection through time (data file S6). Nonetheless, both genes in the United Kingdom and CD96 in Australia have a significant signal of positive selection before the release of RHDV.

We observed large changes in allele frequency across the two time periods (Fig. 3E). Averaging across populations during the first phase of selection, we estimate that the per-year selection coefficient was 0.16 for IFN-α21A, 0.13 for FCRL3, 0.08 for CD96, and 0.07 for the MHC (data file S6). These analyses assume that the alleles were additive, but estimates of the timing and strength of selection remained similar if we assumed a recessive or dominant mode of inheritance.

Selection on the immune system

The innate immune response provides the first line of defense against viruses. It is activated by the secretion of cytokines, including interferon-α (IFN-α), which binds to receptors on uninfected neighboring cells and induces an antiviral state. To circumvent this, MYXV encodes potent inhibitors of the IFN response, including the double-stranded RNA binding protein M029 (20, 21). Among the alleles with the most significant increases in frequency in our dataset, three nonsynonymous variants segregate as a haplotype in the IFN-α21A gene (IFN-α21A) (Fig. 3C). To evaluate the roles of these SNPs in the antiviral activity of IFN-α21A, we synthesized the two corresponding protein isoforms and tested their antiviral effect by measuring MYXV replication in a rabbit cell line. Neither isoform of IFN-α21A affected the replication of the wild-type MYXV (Lausanne strain) (fig. S6a); however, both significantly reduced the replication of an attenuated strain of MYXV (the M029 mutant) (22) (Fig. 4A and fig. S6b). Moreover, we found that the isoform of IFN-α21A that had been favored by selection (varIFN-α21A) more strongly inhibited the replication of M029 mutant MYXV. This indicates that natural selection has increased the potency of the IFN response in modern rabbit populations that have coevolved with MYXV. The selected allele also had an antiviral effect on vesicular stomatitis virus (VSV), an RNA virus (Fig. 4B). This suggests a general increase in the protein’s antiviral activity, and selection by MYXV may have increased resistance to other viruses, including RHDV. Although this effect was not apparent with wild-type MYXV, this may be a limitation of the cell culture–based experiment, as in vivo innate immune responses involve the coordinated action of multiple cytokines across many tissues. In such circumstances, it is possible that the isoform of IFN-α21A that has been favored by selection contributes to attenuating wild-type MYXV infection.

Fig. 4 The effect of IFN-α21A and VPS4 on viral titers.

(A and B) Wild-type IFN-α21A and variant IFN-α21A (varIFN-α21A) were added at different concentrations to cell cultures before infection with (A) the MYXV M029 knockout mutant (M029KO) and (B) VSV. Viral titer was measured 1 hour postinfection and 24 or 48 hours postinfection. Error bars show SEM. Statistical significance between wild and mutated IFN treatments was inferred with Student’s t test across three replicate assays (*P < 0.05; **P < 0.01). FFU, focus-forming units. (C) Human embryonic kidney (HEK) 293 cell lines stably expressing the wild-type isoform of human VPS4 or a dominant-negative VPS4 (EQ mutant) under the control of ponasterone A (PonA). The HEK293 nontransfected cell line (parental) was included as an additional control. Cells were either untreated (PonA−) or pretreated (PonA+) with 1 μM PonA for 20 to 24 hours and then infected with wild-type MYXV expressing a red fluorescent protein (vMyx-tdTomato) at a MOI of 10. The percentage of infected cells (tdTomato+) was assessed by flow cytometry. Error bars show SEM. (D) Fluorescence microscope images of VPS4 wild-type and VPS4 EQ mutant HEK293 cells pretreated with PonA (20 to 24 hours), 48 hours postinfection with vMyx-tdTomato (MOI of 10). The live-cell images were taken using an inverted fluorescence microscope at 10× magnification.

We found a strong population-specific signal of selection on a nonsynonymous variant in CD200-R, which encodes the receptor for the negative regulator of innate immune responses CD200. The selected allele was not observed in any of the historical populations, but it increased to a frequency of 56% in the modern French population (Fig. 3D). The selected residue is part of the CD200 binding interface (23) (fig. S7). Mutation of this site in the human protein does not prevent CD200 binding (24), but it may affect binding affinity in a quantitative way. Therefore, any effects of this variant may occur via modulation of the rabbit CD200–CD200-R interaction.

Alongside the innate immune response, hosts mount an adaptive antiviral immune response mediated by MHC class I proteins. Polymorphisms in the peptide-binding region of MHC proteins affect the repertoire of peptides they can present and are frequently associated with variation in susceptibility to infectious disease (25). We found multiple SNPs in the MHC region that have shown positive selection across the three populations (Fig. 3C). However, the top hits in this region were located in an intergenic region of 35.7 kb between two MHC class I genes (gray shaded box in Fig. 3C). We failed to sequence the regions encoding the peptide-binding regions of the proteins for most individuals (either the sequence capture or the read mapping failed). Therefore, the variants that have increased in frequency in the intergenic region may be in LD with variants in the protein coding sequence that affect susceptibility to MYXV.

Other immunity genes that exhibit parallel changes in frequency include four nonsynonymous variants in FCRL3 (Fc receptor-like 3), which encodes a cell surface receptor that inhibits the activation of B cells (26). Three highly significant nonsynonymous variants are in CD96, which encodes a transmembrane receptor of T and natural killer cells involved in modulating immune responses (27, 28). The four SNPs with the highest likelihood in the entire dataset were all intronic variants present in the MFSD1, which encodes a transmembrane protein in the major facilitator superfamily. Several paralogs of this gene have been identified in a genome-wide RNA interference (RNAi) screen for genes affecting the replication of vaccinia virus, a poxvirus related to MYXV (29). By using RNAi, we found that this gene also modulated MYXV titers in rabbit cells (fig. S6e).

Selection on proviral genes

Resistance to viruses can evolve not only through the improvement of antiviral defenses but also through changes in host proviral proteins that viruses hijack for their own benefit. Among our significant associations was an SNP in the 3′ untranslated region of VPS4 (Fig. 3C). This gene has no known role in the replication cycle of poxviruses but is required for the envelopment of herpes simplex virus (HSV) in the cytoplasm (30). To evaluate the role of VPS4 in MYXV replication, we took advantage of a human cell line expressing a dominant-negative form of the VPS4 protein (30). By comparing the effects of MYXV in wild-type and VPS4 dominant-negative human cell lines, we found that the latter strongly inhibited MYXV replication (Fig. 4, C and D, and fig. S6, c and d). The effect at a high multiplicity of infection (MOI) suggests that VPS4 may affect earlier replication steps in MYXV than in HSV. Therefore, the VPS4 protein is proviral, and it is possible that selection by MYXV altered its expression in the modern rabbit populations.

The proteasome is required both for poxvirus uncoating after cell entry (31) and for processing of antigenic peptides, and three variants in the gene PSMG3, which encodes proteasome assembly chaperone 3, increased 83% in frequency in France, with smaller changes in other populations. By using Sanger sequencing, we found that these SNPs were associated with a 7–base pair (bp) insertion in the first exon and a 50-bp deletion spanning the first intron and second exon. However, because genome annotation of this gene is incomplete, the importance of these changes is unclear.


The myxomatosis pandemic caused massive mortality in rabbit populations, leading to the evolution of genetic resistance to the disease in Australia, France, and the United Kingdom. We have found that over the last 60 years, the convergent phenotype of viral resistance in these populations was also accompanied by parallel genetic changes. This was a consequence of natural selection acting on standing genetic variation that was present in the ancestral rabbit populations in continental Europe and was retained in the subsequent colonization process of the United Kingdom and Australia. The presence of this variation likely explains the rapid development of resistance to myxomatosis observed in rabbit populations almost immediately after the first outbreaks and may frequently be critical to allow populations to respond to new pathogens. This is in marked contrast to the evolution of the virus, where parallel changes in MYXV virulence do not have a common genetic basis (15).

Despite rapid phenotypic evolution, only 2% of the alleles favored by selection have reached fixation in any of the modern populations (Fig. 3C) (3 out of 193 SNPs). Together with our estimates of selection coefficients and the moderate antiviral effect of the three IFN SNPs, this suggests that genetic resistance to myxomatosis is a polygenic trait resulting from the cumulative effect of multiple alleles shifting in frequency across the genome, as opposed to a few major-effect changes to the immune response. Such adaptation is likely to result in a gradual “quantitative” increase in resistance. When resistance reduces the within-host replication rate of the virus rather than completely preventing infection, selection will favor increases in viral virulence (32). Consistent with this prediction, it has recently been observed in wild populations of rabbits that the decline in virulence seen in the years immediately after the virus was released has been reversed and highly virulent viral genotypes have emerged (7).

The evolution of resistance to MYXV is associated with enhanced innate antiviral immunity (6). Homologs of CD96, FCRL, CD200-R, and IFN-α all play a regulatory role in the innate immune response, with effects on lymphocyte proliferation and activation (21, 26, 27, 33). The increased virulence seen in recent MYXV isolates is associated with the virus becoming highly immunosuppressive, causing the loss of cellular inflammatory responses, lymphocytes, and neutrophils (7). Therefore, viral evolution may have found ways to counter the effect of many of the genetic adaptations that we have observed. In conclusion, our results reveal how standing genetic variation in the immune system allowed populations to rapidly evolve resistance to a new and highly virulent pathogen and describe the molecular and genetic basis of an iconic example of natural selection.

Supplementary Materials

Materials and Methods

Figs. S1 to S7

Tables S1 to S4

Data Files S1 to S10

References (3560)

References and Notes

Acknowledgments: We thank the Australian Museum, American Museum of Natural History, Booth Museum of Natural History, Natural History Museum (London), Museum of Comparative Zoology (Harvard University), Musée des Confluences, Muséum National d’Histoire Naturelle (Paris), Museum Victoria, Queensland Museum, Museum of Zoology (University of Michigan), Smithsonian Institution National Museum of Natural History, and all the curators and museum technicians who generously sampled and provided historical samples; M.-D. Wandhammer from Musée Zoologique de la Ville de Strasbourg, who helped to track historical French samples; and the British Association for Shooting and Conservation (BASC), A. Holroyd, S. Whitehead, and all volunteers who contributed with modern rabbit samples. P. Kerr provided rabbit samples from NSW (Australia) from before the emergence of RHDV. P. Elsworth and W. Dobbie contributed Queensland post-RHDV samples. C. Crump provided human cell lines expressing VPS4. R. Fonseca provided valuable advice concerning ancient DNA bioinformatics. Funding: This work was funded by grants from the Programa Operacional Potencial Humano–Quadro de Referência Estratégica Nacional funds from the European Social Fund and Portuguese Ministério da Ciência, Tecnologia e Ensino Superior to M.C. (IF/00283/2014/CP1256/CT0012), to P.J.E. (IF/00376/2015), and to J.M.A. (SFRH/BD/72381/2010). J.M.A. was supported by a travel grant from the Middleton Fund (Cambridge) to undertake work in the Centre of GeoGenetics (Copenhagen). A.M., F.M.J., and L.L. were supported by the European Research Council (grants 647787-LocalAdaptation, 281668, and 339941-ADAPT, respectively). M.T.P.G. was funded by Lundbeck Foundation grant R52-A5062 (Pathogen Palaeogenomics). The McFadden lab is supported by National Institutes of Health (NIH) grant R01 AI080607. S.C.G. holds a Sir Henry Dale Fellowship, cofunded by the Wellcome Trust and the Royal Society (098406/Z/12/Z). Competing interests: None declared. Author contributions: J.M.A., M.C., P.J.E., N.F., and F.M.J. designed and led the study. J.M.A., P.F.C., N.W., S.A., J.P.D., and M.T.P.G. generated sequencing data. J.M.A., T.S., D.J.B., S.J.F., S.M., W.J.P., G.Q., and A.K.S. did fieldwork and extracted DNA. A.L.d.M., M.M.R., S.C.G., L.B., and G.M. performed the experiments in cell culture. J.M.A., M.C., J.Y.C., L.L., A.E., A.M., F.G.V., R.N., and F.M.J. analyzed the data. J.M.A., M.C., and F.M.J. wrote the paper with input from the other authors. All authors approved the manuscript before submission. Data and materials availability: Original sequence data are available in the Sequence Read Archive under BioProject PRJNA393806 (SRP118358). The variant calls are available in the Cambridge Data Repository (34). French modern samples are available from S.M. under a materials agreement with the Office National de la Chasse et de la Faune Sauvage.
View Abstract

Stay Connected to Science

Navigate This Article