The attenuation of myxoma virus (MYXV) following its introduction as a biological control into the European rabbit populations of Australia and Europe is the canonical study of the evolution of virulence. However, the evolutionary genetics of this profound change in host-pathogen relationship is unknown. We describe the genome-scale evolution of MYXV covering a range of virulence grades sampled over 49 years from the parallel Australian and European epidemics, including the high-virulence progenitor strains released in the early 1950s. MYXV evolved rapidly over the sampling period, exhibiting one of the highest nucleotide substitution rates ever reported for a double-stranded DNA virus, and indicative of a relatively high mutation rate and/or a continually changing selective environment. Our comparative sequence data reveal that changes in virulence involved multiple genes, likely losses of gene function due to insertion-deletion events, and no mutations common to specific virulence grades. Hence, despite the similarity in selection pressures there are multiple genetic routes to attain either highly virulent or attenuated phenotypes in MYXV, resulting in convergence for phenotype but not genotype.
The text-book example of the evolution of virulence is the attenuation of myxoma virus (MYXV) following its introduction as a biological control into the European rabbit populations of Australia and Europe in the 1950s. However, the key work on this topic, most notably by Frank Fenner and his colleagues, occurred before the availability of genome sequence data. The evolutionary genetic basis to the major changes in virulence in both the Australian and European epidemics is therefore largely unknown. We provide, for the first time, key details on the genome-wide changes that underpin this landmark example of pathogen emergence and virulence evolution. By sequencing and comparing MYXV genomes, including the original strains released in the 1950s, we show that (i) MYXV evolved rapidly in both Australia and Europe, producing one of the highest rates of evolutionary change ever recorded for a DNA virus, (ii) that changes in virulence were caused by mutations in multiple genes, often involving losses of gene function due to insertions and deletions, and that (iii) strains of the same virulence were defined by different mutations, such that both attenuated and virulent MYXV strains are produced by a variety genetic pathways, and generating convergent evolution for phenotype but not genotype.
Citation: Kerr PJ, Ghedin E, DePasse JV, Fitch A, Cattadori IM, Hudson PJ, et al. (2012) Evolutionary History and Attenuation of Myxoma Virus on Two Continents. PLoS Pathog 8(10): e1002950. https://doi.org/10.1371/journal.ppat.1002950
Editor: Bruce R. Levin, Emory University, United States of America
Received: April 26, 2012; Accepted: August 22, 2012; Published: October 4, 2012
Copyright: © Kerr et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Funding: This work was funded in part by grant R01 AI093804-01A1 from the National Institute of Allergy and Infectious Diseases, National Institutes of Health (http://www.niaid.nih.gov/Pages/default.aspx). We thank members of the Research and Policy in Infectious Disease Dynamics Program of the Science and Technology Directorate, Department of Homeland Security (http://www.dhs.gov/index.shtm) and the Fogarty International Center, National Institutes of Health (http://www.fic.nih.gov/Pages/Default.aspx)for support. We also thank the Huck Institutes of Life Sciences, The Pennsylvania State University (http://www.huck.psu.edu/), for seed funding. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
The classic model of pathogen evolution following a species jump is the introduction of the lethal myxoma virus (Poxviridae, genus Leporipoxvirus) into the European rabbit (Oryctolagus cuniculus) populations of Australia as a biological control . MYXV causes a largely innocuous cutaneous fibroma in its natural host, the South American tapeti (Sylvilagus brasiliensis). However, in European rabbits the same virus causes the systemic lethal disease myxomatosis. A strain of this highly virulent virus – Standard Laboratory Strain (SLS) – was successfully introduced into wild European rabbits in Australia in 1950 and over the next five years spread across the entire continental range of the rabbit. With outstanding foresight, Fenner and colleagues tracked the phenotypic evolution of virulence in the Australian epidemic by testing the lethality of viral isolates taken from the field over the subsequent decades in laboratory rabbits (of the same species as those infected in the wild).
SLS had a case fatality rate (CFR) estimated at 99.8% . However, within two years of the introduction of SLS, and despite the ongoing release of virulent viruses, slightly attenuated MYXV strains came to dominate field populations. Although they still killed 90–99% of infected rabbits, these lower virulence strains allowed infected rabbits to survive for longer, increasing the probability of transmission from skin lesions by mosquito vectors , . MYXV virus does not replicate in the vector , , so that transmission is dependent on high titres in cutaneous lesions (≥107 infectious particles/g) and the duration of survival of the infected rabbit. Highly attenuated viruses could be quickly controlled by the rabbit immune system and so tended to have relatively low transmissibility . For the next 30 years the majority of MYXV were of intermediate virulence with CFRs of 70–95% in laboratory rabbits . At the same time, natural selection acted on the wild rabbit population, resulting in the emergence of animals resistant to myxomatosis , , likely through an enhanced innate immune response , . This effectively reduced the virulence of field strains of MYXV in wild compared to laboratory rabbits. Uniquely, this evolutionary experiment was repeated in Europe following the release of a separate South American strain of MYXV (Lausanne strain; Lu) into France in 1952. The evolutionary outcomes were strikingly similar to those in Australia, with the emergence of attenuated virus strains and selection for resistant rabbits .
These remarkable natural experiments in biological control have become the bedrock of the mathematical theory of virulence evolution –, revealing much about the relationship between virulence and transmissibility. The hypothesis is that a reduction in virulence was selectively favoured because, by killing hosts so rapidly, highly virulent viruses had shorter infectious periods and hence lower fitness. The emergence of attenuated viruses may also have aided the selection of resistant rabbits , , in turn driving the evolution of more virulent MYXV in subsequent decades in an archetypal host-pathogen arms race , .
Despite the canonical position of MYXV in studies of virulence evolution, the critical body of work occurred in the pre-genomic era. Today, we know that at least 20 genes encoding proteins with immune evasion or host range functions have a demonstrated role in MYXV virulence based on gene knock-out studies, and a further 20 proteins are predicted to have similar functions , . However, both the pattern of molecular evolution and the genetic basis to the changes in virulence in both the Australian and European epidemics is largely unknown. Here we use comparative genomics to illuminate these uniquely important evolutionary events.
Evolutionary history and dynamics of MYXV
The genome of the Lu strain of MYXV comprises 161.8 kb of double-strand (ds) DNA with terminal inverted repeats (TIRs) of 11,577 bp. MYXV contains 158 unique open reading frames, 12 of which are duplicated in the TIRs. The central region of the genome encodes genes involved in viral replication that are strongly conserved among poxviruses. Genes located towards the termini and within the TIRs are less conserved, have roles in immune evasion and host-range, and have been shown to be critical for virulence in European rabbits , .
We obtained complete or near complete genome sequences of 22 isolates of MYXV (Table 1), with the Lu strain  also sequenced as a quality control. These data were combined with the published sequence of a single Spanish isolate (MYXV/6918) sampled in 1995 . A maximum likelihood phylogenetic analysis of these 24 complete MYXV genomes revealed a division between those viruses released in Australia (SLS progenitor) and Europe (Lu progenitor), a strong temporal structure, such that those viruses sampled earlier in time (1940s and 1950s) fall closer to the root of the tree, and frequent changes in virulence grade over the sampling period (Figure 1a).
A: Maximum likelihood (ML) phylogeny of 24 complete genome sequences of MYXV. Viruses are color-coded according to virulence grade (defined in Table 1); grade 1 = red, grades 3–4 = green, grade 5 = blue. The tree is rooted between the two oldest strains in the collection – Lausanne (1949) and SLS (1950) – which were used to seed the European and Australian epidemics, respectively (marked by asterisks). The same root position is obtained under the mid-point method and from the Maximum Clade Credibility (MCC) tree in BEAST . All horizontal branch lengths are drawn to a scale of nucleotide substitutions per site, and bootstrap values >70% are shown. B: rate of the nucleotide substitution (range of 95% HPD values) for the real MYXV genomic data set (on right in red) compared to 20 sequence data sets in which the year of sampling has been randomized among the taxa. Note that the lower 95% HPD value for the randomized data tends to a zero rate. C: regression of root-to-tip genetic distances against year of sampling for the 24 complete genomes of MYXV and inferred from the ML phylogenetic tree. The correlation coefficient is 0.98 and the slope of the line, which is an estimator of the mean substitution rate, is 1.0×10−5 subs/site/year.
A parsimony analysis reconstructed 482 mutations across the MYXV phylogeny (mean pairwise distance among isolates = 0.05%), of which 207 were synonymous, 249 nonsynonymous, and 26 occurred in non-coding or intergenic regions. These data were also characterized by a high consistency index (0.89) indicative of limited homoplasy, and no recombination was observed. A Bayesian Markov chain Monte Carlo analysis  recovered a mean evolutionary rate of 9.6×10−6 nucleotide substitutions per site, per year (95% HPD values, 8.3–10.9×10−6 subs/site/year) under a strict molecular clock, with almost identical estimates obtained using a relaxed clock. Because both clock models have essentially identical likelihoods, and the coefficient of rate variation includes zero, we conclude that MYXV has evolved in a strongly clock-like manner. In addition, significant temporal structure was observed using both a phylogenetic randomization test, in which MYXV sequences with randomized times of sampling produced a distribution of substitution rates significantly different to that observed in the real MYXV data (Figure 1b), and in a regression analysis of root-to-tip genetic distances against sampling times (Figure 1c; correlation coefficient = 0.98, although this value will be inflated to some extent by phylogenetic non-independence in the data). In sum, there is measurable, and therefore relatively rapid, viral evolution over the sampling period. Accordingly, those Australian MYXV isolates sampled from the 1990s share a common ancestor that existed between 1964–1971.
Evolution of MYXV virulence in Australia
Fenner and colleagues categorised virulence on a scale from grade 1 (extreme lethality) to grade 5 (case fatality rate <50%) (Table 2). Nine of the 19 Australian MYXV isolates sequenced here have been phenotyped for virulence (Figure 1a; Table 1). These cover almost 50 years of viral evolution and the full range of virulence, including viruses that are as lethal as the progenitor SLS (grade 1) strain, as well as those with attenuated phenotypes and prolonged survival times (grade 5). Three of the phenotyped isolates were sampled during the initial epidemic; the Glenfield strain (Dubbo, Feb. 1951; grade 1), KM13 (Corowa, Dec. 1952; grade 3) and Uriarra (Feb. 1953; originally characterised as grade 4, but grade 5 in more recent testing ). Glenfield was more virulent than SLS when tested in resistant rabbits demonstrating that virulence can increase as well as decrease . The remaining six phenotyped MYXV viruses were isolated between 1991 and 1999: Bendigo, Wellington, BRK (grade 1), SWH, Gung (grade 4) and Meby (grade 5) . Among all nine phenotyped viruses there were 109 nonsynonymous mutations and 18 insertion-deletions (indels; ranging from 1–92 bp in length) in 67 open reading frames compared to the ancestral SLS sequence (Table S1). Strikingly, only two nonsynonymous (amino changes P76H in M140R and L204S in M153R) and two synonymous mutations are shared by all 19 Australian viruses compared to SLS. The genomic locations of the mutations in all 19 Australian MYXV isolates are shown in Figure 2. A total of 23 nonsynonymous, nine synonymous and one intergenic mutation fall on the branch separating the Australian viruses sampled during the 1990s from their 1950s relatives (Table S1). Interestingly, two indels are also reconstructed to fall on the branch separating the 1950s and 1990s Australian viruses; an A insert in M009L which disrupts this ORF (and that reverts in Bendigo), and a G insert in M083L that corrects a reading frame disruption in SLS and the other early Australian viruses. As MYXV will have clearly been subject to a variety of selection pressures over our sampling period, including adaptation to the European rabbit following the species jump from the tapeti, it is likely that only a subset of these mutations would be responsible for virulence evolution.
Open reading frames and their direction of transcription are shown as arrows, and gene identities are also indicated. Light arrows above the sequence indicate the terminal inverted repeats of 11577 bp at each end of the genome. Vertical lines show the location of mutations in the Australian isolates of MYXV compared to the SLS progenitor strain that initiated the Australian epizootic of myxomatosis.
Some mutations, and particularly indels, have a clear impact on virulence. For example, the highly attenuated Uriarra strain contains a C nucleotide insertion in M005L/R that disrupts the reading frame. M005 is an E3 Ub ligase that manipulates cell cycle progression and inhibits cell death; as it is known to be critical for virulence, this indel is likely the main mutation responsible for the attenuation of Uriarra , . The grade 5 Meby strain has a 73 bp deletion in the conserved region of M153R ,  which leads to read-through of the normal stop codon and an altered C-terminal protein sequence which is no longer highly acidic. M153R is a key virulence gene , encoding an E3 Ub ligase with an N-terminal RING-CH domain responsible for downregulation of MHC-1, CD4, ALCAM/CD166 and Fas/CD95 on the membrane of infected cells –. Hence, this mutation is predicted to impair the capacity of MYXV to interfere with host destruction of infected cells, resulting in an attenuated phenotype.
Aside from these particular indels, the mutational events responsible for other changes in virulence are less obvious. To frame our analysis of virulence evolution, we consider the two simplest hypotheses to account for this changing MYXV phenotype in Australia. The most parsimonious explanation is that the mutations that caused attenuation in the early 1950s were maintained over the subsequent half century; hence, the phenotypic diversity in the 1990s was due to the co-circulation of virulence determinants present in the ancestral strain and attenuation mutations favoured by natural selection in the period immediately after viral release (the ‘attenuation-only’ model). A competing hypothesis is that the early MYXV strains fixed mutations that resulted in an attenuated phenotype, and that the appearance of the highly virulent grade 1 viruses in the 1990s was due to mutations which restored virulence. Under this ‘attenuation-restoration’ model patterns of phenotypic evolution over the half century following release would be due to changes in the relative frequency of strains with ancestral virulence, attenuation mutations that occurred soon after release, and strains with de novo mutations that restored high virulence.
Under either scenario, our genomic data can account for the observed pattern of virulence evolution only if phenotypic convergence has been achieved through multiple genetic routes. Attenuation is never associated with convergence at the genotypic level: among coding mutations, none uniquely define our single grade 3 virus, none are uniquely shared by both grade 4 viruses, and none are uniquely shared by both grade 5 viruses (Table S1). For example, the two grade 5 viruses, Meby and Uriarra, only share those nonsynonymous mutations that are fixed in all Australian strains, including those of high virulence. Similarly, no mutations are uniquely shared by the three grade 1 viruses isolated in the 1990s. Thus, if novel mutations did restore virulence to these viruses, these are specific to each. A picture of multiple genetic pathways is also observed at the genic level: none of the virulence grades can be defined by uniquely mutated genes.
Although it is clear that there are multiple routes to attenuation/virulence, mapping virulence grade on to our phylogeny using a parsimony procedure (Figure 1a) results in ambiguous state characteristics for most nodes, so that it impossible to determine the exact direction of virulence evolution. We note, however, that the ‘attenuation-only’ model cannot easily account for the phenotype of KM13, a grade 3 virus recovered late in 1952. Although KM13 differs from SLS at five amino acids, all are also found in at least one grade 1 virus. In addition, KM13 falls in a phylogenetic position that is directly ancestral to the grade 1 strains sampled in the 1990s (Figure 1a), suggesting that there was global fixation of at least some attenuating mutations as predicted by the ‘attenuation-restoration’ model.
Evolution of MYXV in Europe
To further explore the evolution of MYXV we sequenced two British strains that have been phenotyped for virulence. Cornwall is a grade 1 virulence virus isolated in April 1954 soon after the release of MYXV in Britain , while Nottingham attenuated (Nott), a grade 5 virulence virus, was obtained from an infected rabbit in April 1955 . We compared these viruses to both the parental Lu strain (Table S2) and to MYXV/6918, a more recent grade 5 virus from Spain.
The most likely explanation for the profound attenuation of Nott is the insertion of a TG dinucleotide in M150R which leads to a premature stop codon after amino acid 196. M150R encodes an inhibitor of NFκB, and a virus in which this gene is insertionally inactivated is highly attenuated , . Thus, as with the M153 mutation in Meby, attenuation has apparently been achieved by the virus reducing its capacity to manipulate host immunity. This TG insertion is not present in any of the Australian viruses, and although six of the phenotyped Australian strains share a single substitution in this gene, they span the full virulence spectrum (Table S1). Similarly, single nucleotide insertions in M135R and M148R are the most likely cause of the attenuation of MYXV/6918 , none of which are seen in the Australian viruses. Hence, despite the parallel evolution of MYXV attenuation in Australia and Europe, different mutations in different genes were responsible for the reduction in virulence in both cases. Indeed, of the 482 mutations in the MYXV phylogeny, only two occur in parallel between the Australian and European epidemics. Both are nonsynonymous – a A37V amino acid replacement in M003.1 and a P173S replacement in M150 – but they are found in Australian viruses spanning the virulence spectrum and occur on the long branches separating the 1990s from the 1950s viruses.
Despite their importance as agents of disease, far less is known about the evolutionary dynamics of large DNA viruses than of RNA viruses, and what research has been done often relies on the assumption that the viruses in question have co-diverged with their hosts over long periods of evolutionary time . Our sampling of 49 years of MYXV evolution therefore enables a unique snap-shot of the pace of viral evolution in real-time. Accordingly, our genome scale data provides strong evidence for the rapid evolution of MYXV in both Australia and Europe. The rate of nucleotide substitution we estimate here (∼1×10−5 subs/site/year) is one of the highest ever reported for a dsDNA virus, but similar to that documented in another poxvirus – Variola (VARV) virus – the agent of smallpox . These rates are approximately three orders of magnitude higher than those estimated in the dsDNA vertebrate herpesviruses assuming virus-host co-divergence , and where there is a marked lack of temporal structure in the short-term . It is therefore possible that poxviruses are characterized by relatively high background mutation rates, which may in part explain the clock constancy. However, direct estimates of mutation rate are only available for a limited number of dsDNA viruses, with none from the Poxviridae . Alternatively, as MYXV was likely subject to strong selection for enhanced transmissibility and a variety of other phenotypic traits during our sampling period, it is possible that its substitution rate has been elevated by adaptive evolution. Nonsynonymous mutations were commonplace in the MYXV phylogeny, and particularly on the branch separating the 1990s Australian viruses from those sampled forty years earlier. However, their frequency of occurrence (∼55% of all mutations in coding regions) is still compatible with a predominantly neutral evolutionary process. Indeed, the small number of point mutations in individual MYXV genes precluded an analysis of gene and site-specific selection pressures.
These data are also notable in that they show that the attenuation of the grade 1 virulence of SLS in Australia was achieved by unique mutations affecting a variety of different genes. This may be a consequence of the initially wide geographic spread of MYXV which, in combination with population devastation and seasonal population bottlenecks, would have caused regionally successful virulence mutants to arise and be lost stochastically. As rabbits subsequently evolved resistance, it seems likely that a variety of virulence-restoring mutations would be locally selected among attenuated strains to maintain relative transmissibility, which was maximal for grade 4 phenotypes when measured experimentally in non-resistant laboratory rabbits . Such a scenario is compatible with our MYXV phylogeny (Figure 1a) in which relatively few mutations are fixed on internal nodes of the tree, particularly during the 1950s, with more mutations located on the branches leading to individual variants and indicative of localized evolution.
More broadly, our genome scale analysis of MYXV evolution on two continents reveals that there are multiple routes to viral virulence or attenuation, of which we have likely sampled only a small number. Hence, there is convergence for phenotype but not genotype. It is likely that many of the mutational changes associated with virulence have subtle effects and are subject to complex, but as yet undefined, epistatic interactions. Such complexity could in part reflect the evolutionary arms race involving selection for resistance in rabbits and counter-selection for virulence in the virus to maintain transmissibility and competitiveness, and which led to the reversal of attenuating mutations, mutations that compensate for attenuating mutations, or those that increase virulence by novel pathways. Indeed, it is striking that both highly virulent and highly attenuated viruses appear to be successful at least under some local conditions, although the BRK and Bendigo viruses that possessed grade 1 virulence in laboratory rabbits were of grade 4 and 5 virulence, respectively, when tested in wild resistant rabbits , and the progenitor SLS strain was also effectively a grade 5 virus in resistant rabbits . These laboratory tests were performed under cage or pen conditions and it is likely that virulence would be higher in the field due to a variety of biotic and abiotic stress factors, and which may in part explain the apparent success of attenuated viruses such as SWH or Gung. In transmission trials using viruses of virulence grades 1, 3 and 5, rabbits that died quickly following infection and those that survived infection were both poor sources of viral transmission by fleas; most transmission occurred from those rabbits that had prolonged survival times but ultimately died . Stress due to nutrition, cold, and other diseases may have increased the effective virulence of a field strain of MYXV in individual rabbits, thereby enhancing the transmission of attenuated viruses and decreasing the transmission of more virulent strains. In the fragmented rabbit populations of today which are still considerably reduced from those in 1950, it may be that both local ecological conditions and virulence plasticity are critical for virus success. Clearly, multiple genetic variants of MYXV may appear in relatively small geographic areas, with the predominant viral type shifting from year to year .
Although in vitro studies of virulence determinants in viruses are commonplace, MYXV occupies a unique position in studies of virulence evolution because this evolutionary process occurred in natural populations and virulence was experimentally assayed in a controlled and relevant common garden species (laboratory rabbits of the same species as those infected in the wild). Our comparative genome sequence data suggests that the large and complex genome of MYXV provides the plasticity for multiple routes to attenuation, and multiple and complex routes back to virulence. Such genetic flexibility sits in contrast to RNA viruses such as West Nile , influenza A , and HIV-1 –, in which it has been suggested that a limited number of mutations may control virulence. Key differences between large dsDNA viruses like MYXV and these RNA viruses are that the latter possess a limited number of genes that commonly encode multiple functions, and experience infrequent gene duplication and lateral gene transfer, all of which will limit the number of viable evolutionary pathways . In light of such evolutionary complexity, it will be necessary to employ a systematic reverse genetic approach to precisely determine the contribution of individual mutations to the virulence and attenuation of this important model of disease emergence. More generally, the complex relationship between genotype and phenotype observed here suggests that it will often be difficult to predict the future course of virulence evolution in emerging pathogens from comprehensive genome sequence data.
Materials and Methods
Virus growth and DNA extraction
All virus isolates analyzed here, as well as their virulence, are described in Table S1. Viruses were passaged twice in RK13 cells to produce working stocks. Concentrated virions were prepared from infected RK13 cells and DNA extracted as described previously, with the exception that DNA was ethanol precipitated rather than dialysed .
Sequencing and genome assembly
We obtained and assembled complete or almost complete genome sequences for 23 isolates of MYXV using the Roche GS-FLX sequencing platform. Whole genome shotgun libraries were created for each sample by fragmenting the genomic DNA using a nebulization technique. The genomic libraries were sequenced in multiple runs on the GS-FLX using the Titanium chemistry and a 16-well gasket to separate each sample on the picotiter plate (PTP). Approximately 50–60% of all the sequence reads generated corresponded to host genomic DNA so that the majority of the samples had to be sequenced twice leading to an average sequence coverage of 30× (range of 11–60X). Sequence reads were assembled into contigs using the comparative assembler AMOScmp  which uses a reference genome (Lu; GenBank accession NC_001132) as a guide for the assembly of related genomes, as well as with gsMapper (Roche). We also re-sequenced a Lu strain as a control of our sequencing quality. In general, the assembly produced two virus contigs per genome; one at ∼11,500 bp corresponding to the collapsed terminal inverted repeats and another, at ∼140,000 bp, corresponding to the core portion of the genome. For some genomes lighter coverage of various regions led to gaps, or suspect areas, that were closed by targeted Sanger sequencing. All high-quality reported single nucleotide polymorphisms (SNPs) were also confirmed by Sanger sequencing. SeqMan Pro (Lasergene) was used to incorporate these finishing reads into the larger assembly, followed by visual confirmation of coverage. All sequence data generated here has been submitted to GenBank and assigned accession numbers JX565562–JX565584.
We first screened for recombination in the 24 complete genome sequences of MYXV (alignment length of 163,538 bp) using the RDP, GENECOV, BOOTSCAN programs available within the RDP3 package . As no recombination was observed, a phylogenetic tree was then inferred using the maximum likelihood (ML) method available in the PAUP* package , and employing Tree Bisection and Reconnection (TBR) branch-swapping. Because of the very small number of nucleotide substitutions across the MYXV genome, such that any multiple substitution will be very limited, we utilized the HKY85 model of nucleotide substitution with default parameters. A bootstrap resampling analysis (1000 replications) was used to determine the support for key nodes. Finally, we employed a parsimony protocol to reconstruct the nucleotide changes along each branch of the ML tree, again utilizing the PAUP* package.
To infer rates of nucleotide substitution we employed the Bayesian Markov Chain Monte Carlo (MCMC) method available in the BEAST package . This analysis employed both strict (i.e. constant) and relaxed (uncorrelated lognormal) molecular clocks and utilized the year of sampling of each viral isolate. We again employed the HKY85 model of nucleotide substitution as well as a Bayesian skyline coalescent prior. MCMC chains of 100 million generations (with 10% burn-in) were run multiple times until convergence was achieved for all parameters. Statistical uncertainty was reflected in values for the 95% Highest Probability Density (HPD). To determine whether there was sufficient temporal structure in the data to accurately infer evolutionary rates, we randomized sampling times 20 times across the MYXV alignment and repeated the BEAST analysis described above. In addition, we plotted root-to-tip genetic distances determined from the ML tree against year of sampling using the Path-O-Gen program (http://tree.bio.ed.ac.uk/software/pathogen/). For this analysis the ML tree is rooted between the oldest Lu and SLS strains, and which is also the deepest split in both BEAST and mid-point rooted trees.
For the BEAST analysis we used two possible dates for the SLS sequence; although the first description of this virus was published in 1911 , it was subjected to regular rabbit passage until 1950 when it was released into Australia. Therefore, the ‘correct’ evolutionary date for SLS falls at some point between 1911 and 1950. Although estimates of the rate of nucleotide substitution were similar under both the 1911 (95% HPD = 7.6–9.0×10−6 subs/site/year; strict clock) and 1950 (95% HPD = 8.3–10.9×10−6 subs/site/year; strict clock) dates, we based our analysis on the 1950 date as it provided a better fit to a strict molecular clock, both with respect to the coefficient of variation in BEAST (values of 0.418–1.076 and 0.00008–0.388 for the 1911 and 1950 dates, respectively) and in the root-to-tip regression (correlation coefficient = 0.92 and 0.98 for the 1911 and 1950 dates, respectively). However, the uncertainty over the correct evolutionary date for SLS does compromise estimates of when this virus shared a common ancestor with Lu.
Sequence differences between SLS and Australian isolates of MYXV. Mutations in viruses with known phenotype are shaded.
Sequence differences between Lu and the European isolates of MYXV sequenced here.
This paper is dedicated to the memory of Frank Fenner (1914–2010), without whose foresight this unique experiment would have gone unexplored. We thank Dr. Andrew Rambaut (University of Edinburgh) for his assistance with the Path-O-Gen analysis.
Conceived and designed the experiments: PJK EG IMC PJH AFR ECH. Performed the experiments: PJK EG JVD AF. Analyzed the data: PJK EG JVD AF ECH. Contributed reagents/materials/analysis tools: PJK EG JVD AF DCT. Wrote the paper: PJK EG IMC PJH DCT AFR ECH.
Fenner F, Ratcliffe FN (1965) Myxomatosis. Cambridge: Cambridge University Press.
- 2. Fenner F, Day MF, Woodroofe GM (1956) Epidemiological consequences of the mechanical transmission of myxomatosis by mosquitoes. J Hyg (Camb) 54: 284–303.
- 3. Fenner F, Marshall ID (1957) A comparison of the virulence for European rabbits (Oryctolagus cuniculus) of strains of myxoma virus recovered in the field in Australia, Europe and America. J Hyg (Camb) 55: 149–191.
- 4. Day MF, Fenner F, Woodroofe GM, McIntyre GA (1956) Further studies on the mechanism of mosquito transmission of myxomatosis in the European rabbit. J Hyg (Camb) 54: 258–283.
- 5. Fenner F (1983) Biological control as exemplified by smallpox eradication and myxomatosis. Proc Royal Soc Lond B 218: 259–285.
- 6. Marshall ID, Douglas GW (1961) Studies in the epidemiology of infectious myxomatosis of rabbits. VIII. Further observations on changes in the innate resistance of Australian wild rabbits exposed to myxomatosis. J Hyg (Camb) 59: 117–122.
- 7. Marshall ID, Fenner F (1959) Studies in the epidemiology of infectious myxomatosis of rabbits. V. Changes in the innate resistance of wild rabbits between 1951 and 1959. J Hyg (Camb) 56: 288–302.
- 8. Best SM, Kerr PJ (2000) Coevolution of host and virus: the pathogenesis of virulent and attenuated strains of myxoma virus in resistant and susceptible European rabbits. Virology 267: 36–48.
- 9. Kerr PJ, McFadden G (2002) The immune response to myxoma virus. Viral Immunol 15: 229–246.
- 10. Anderson RM, May RM (1982) Coevolution of hosts and parasites. Parasitol 85: 411–426.
- 11. Bull JJ (1994) Virulence. Evolution 48: 1423–1437.
- 12. Mackinnon MJ, Gandon S, Read AF (2008) Virulence evolution in response to vaccination: The case of malaria. Vaccine 26: C42–C52.
- 13. Read AF (1994) The evolution of virulence. Trends Microbiol 2: 73–76.
- 14. Dwyer G, Levin SA, Buttel L (1990) A simulation model of the population dynamics and evolution of myxomatosis. Ecol Monog 60: 423–447.
- 15. Boots M, Hudson PJ, Sasaki A (2004) Large shifts in pathogen virulence relate to host population structure. Science 303: 842–845.
- 16. Cameron C, Hota-Mitchell S, Chen L, Barrett J, Cao J-X, et al. (1999) The complete DNA sequence of myxoma virus. Virology 264: 298–318.
- 17. Stanford MM, Werden SJ, McFadden G (2007) Myxoma virus in the European rabbit: interactions between the virus and its susceptible host. Vet Res 38: 299–318.
- 18. Morales M, Ramirez MA, Cano MJ, Parraga M, Castilla J, et al. (2009) Genome comparison of a nonpathogenic myxoma virus field strain with its ancestor, the virulent Lausanne strain. J Virol 83: 2397–2403.
- 19. Drummond AJ, Rambaut A (2007) BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol Biol 7: 214–221.
- 20. Kerr PJ, Jackson RJ, Robinson AJ, Swan J, Silvers L, et al. (1999) Infertility in female rabbits (Oryctolagus cuniculus) alloimmunized with the rabbit zona pellucida protein ZPB either as a purified recombinant protein or expressed by recombinant myxoma virus. Biol Reprod 61: 606–613.
- 21. Douglas GW (1962) The Glenfield strain of myxoma virus. Its use in Victoria. J Agri (Victoria) 60: 511–515.
- 22. Saint KM, French N, Kerr PJ (2001) Genetic variation in Australian isolates of myxoma virus: an evolutionary and epidemiological study. Arch Virol 146: 1105–1123.
- 23. Werden SJ, McFadden G (2008) The role of cell signaling in poxvirus tropism: the case of the M-T5 host range protein of myxoma virus. Biochim Biophys Acta 1784: 228–237.
- 24. Mossman K, Fong Lee S, Barry M, Boshkov L, McFadden G (1996) Disruption of M-T5, a novel myxoma virus gene member of the poxvirus host range superfamily, results in dramatic attenuation of myxomatosis in infected European rabbits. J Virol 70: 4394–4410.
- 25. Collin N, Guérin J-L, Drexler I, Blanié S, Gelfi J, et al. (2005) The poxviral scrapin MV-LAP requires a myxoma viral infection context to efficiently downregulate MHC-I molecules. Virology 343: 171–178.
- 26. Guérin J-L, Gelfi J, Bouillier S, Delverdier M, Bellanger F-A, et al. (2002) Myxoma virus leukemia-associated protein is responsible for major histocompatibility complex class 1 and Fas-CD95 down-regulation and defines scrapins, a new group of surface cellular receptor abductor proteins. J Virol 76: 2912–2923.
- 27. Bartee E, McCormack A, Früh H (2006) Quantitative membrane proteomics reveals new cellular targets of viral immune modulators. PLoS Pathog 2: E107.
- 28. Mansouri M, Bartee E, Gouveia K, Hovey Nerenberg BT, et al. (2003) The PHD/LAP-domain protein M153R of myxoma virus is a ubiquitin ligase that induces the rapid internalization and lysosomal destruction of CD4. J Virol 77: 1427–1440.
- 29. Blanié S, Gelfi J, Bertagnoli S, Camus-Bouclainville C (2010) MNF, an ankyrin repeat protein of myxoma virus, is part of a native cellular SCF complex during viral infection. Virol J 7: 56.
- 30. Camus-Bouclainville C, Fiette L, Bouchiha S, Pignolet B, Counor D, et al. (2004) A virulence factor of myxoma virus colocalizes with NF-kappaB in the nucleus and interferes with inflammation. J Virol 78: 2510–2516.
- 31. Duffy S, Shackelton LA, Holmes EC (2008) Rates of evolutionary change in viruses: Patterns and determinants. Nat Rev Genet 9: 267–276.
- 32. Firth C, Kitchen A, Shapiro B, Suchard MA, Holmes EC, et al. (2010) Using time-structured data to estimate evolutionary rates of double-stranded DNA viruses. Mol Biol Evol 27: 2038–2051.
- 33. McGeoch DJ, Gatherer D (2005) Integrating reptilian herpesviruses into the family Herpesviridae. J Virol 79: 725–731.
- 34. Sanjuán R, Nebot MR, Chirico N, Mansky LM, Belshaw R (2010) Viral mutation rates. J Virol 84: 9733–9748.
- 35. Kerr PJ, Merchant JC, Silvers L, Hood G, Robinson AJ (2003) Monitoring the spread of myxoma virus in rabbit populations in the southern tablelands of New South Wales, Australia. II. Selection of a virus strain that was transmissible and could be monitored by polymerase chain reaction. Epidemiol Infect 130: 123–133.
- 36. Kerr PJ, Perkins HD, Inglis B, Stagg R, McLaughlin E, et al. (2004) Expression of rabbit IL-4 by recombinant myxoma viruses enhances virulence and overcomes genetic resistance to myxomatosis. Virology 324: 117–128.
- 37. Mead-Briggs AR, Vaughan JA (1975) The differential transmissibility of myxoma virus strains of differing virulence grades by the rabbit flea Spilopsyllus cuniculi (Dale). J Hyg (Camb) 75: 237–247.
- 38. Kerr PJ, Hone J, Perrin L, French N, Williams CK (2010) Molecular and serological analysis of the epidemiology of myxoma virus in rabbits. Vet Microbiol 143: 167–178.
- 39. Brault AC, Huang C Y-H, Langevin SA, Kinney RM, Bowen RA, et al. (2007) A single positively selected West Nile viral mutation confers increased avian virogenesis in American crows. Nat Genet 39: 1162–1166.
- 40. Baigent SJ, McCauley JW (2003) Influenza type A in humans, mammals and birds: determinants of virus virulence, host-range and interspecies transmission. Bioessays 25: 657–671.
- 41. Fenyö EM, Albert J, Asjö B (1989) Replicative capacity, cytopathic effect and cell tropism of HIV. AIDS 3: S5–12.
- 42. Connor RI, Ho DD (1994) Human immunodeficiency virus type 1 variants with increased replicative capacity develop during the asymptomatic stage before disease progression. J Virol 68: 4400–4408.
- 43. Kimata JT, Kuller L, Anderson DB, Dailey P, Overbaugh J (1999) Emerging cytopathic and antigenic simian immunodeficiency virus variants influence AIDS progression. Nat Med 5: 535–541.
Holmes EC (2009) The Evolution and Emergence of RNA viruses. Oxford: Oxford University Press.
- 45. Pop M, Phillippy A, Delcher AL, Salzberg SL (2004) Comparative genome assembly. Briefings Bioinformat 5: 237–248.
- 46. Martin DP, Lemey P, Lott M, Moulton V, Posada D, et al. (2010) RDP3: a flexible and fast computer program for analyzing recombination. Bioinformatics 26: 2462.
Swofford DL (2003) PAUP*. Phylogenetic Analysis Using Parsimony (*and other methods). Version 4. Sunderland, Mass: Sinauer Associates.
- 48. Moses A (1911) O virus do myxoma dos coelhos. Mem Inst Osw Cruz 3: 46–53.
- 49. Russell RJ, Robbins SJ (1989) Cloning and molecular characterization of the myxoma virus genome. Virology 170: 147–159.
- 50. Berman D, Kerr PJ, Stagg R. van Leeuwen BH, Gonzalez T (2006) Should the 40-year-old practice of releasing virulent myxoma virus to control rabbits (Oryctolagus cuniculus) be continued? Wildl Res 33: 549–556.