ABSTRACT
Myxoma virus (MYXV) has been evolving in a novel host species—European rabbits—in Australia since 1950. Previous studies of viruses sampled from 1950 to 1999 revealed a remarkably clock-like evolutionary process across all Australian lineages of MYXV. Through an analysis of 49 newly generated MYXV genome sequences isolated in Australia between 2008 and 2017, we show that MYXV evolution in Australia can be characterized by three lineages, one of which exhibited a greatly elevated rate of evolutionary change and a dramatic breakdown of temporal structure. Phylogenetic analysis revealed that this apparently punctuated evolutionary event occurred between 1996 and 2012. The branch leading to the rapidly evolving lineage contained a relatively high number of nonsynonymous substitutions, and viruses in this lineage reversed a mutation found in the progenitor standard laboratory strain (SLS) and all previous sequences that disrupts the reading frame of the M005L/R gene. Analysis of genes encoding proteins involved in DNA synthesis or RNA transcription did not reveal any mutations likely to cause rapid evolution. Although there was some evidence for recombination across the MYXV phylogeny, this was not associated with the increase in the evolutionary rate. The period from 1996 to 2012 saw significant declines in wild rabbit numbers, due to the introduction of rabbit hemorrhagic disease and prolonged drought in southeastern Australia, followed by the partial recovery of populations. It is therefore possible that a rapidly changing environment for virus transmission changed the selection pressures faced by MYXV, altering the course and pace of virus evolution.
IMPORTANCE The coevolution of myxoma virus (MYXV) and European rabbits in Australia is one of the most important natural experiments in evolutionary biology, providing insights into virus adaptation to new hosts and the evolution of virulence. Previous studies of MYXV evolution have also shown that the virus evolves both relatively rapidly and in a strongly clock-like manner. Using newly acquired MYXV genome sequences from Australia, we show that the virus has experienced a dramatic change in evolutionary behavior over the last 20 years, with a breakdown in clock-like structure, the appearance of a rapidly evolving virus lineage, and the accumulation of multiple nonsynonymous and indel mutations. We suggest that this punctuated evolutionary event may reflect a change in selection pressures as rabbit numbers declined following the introduction of rabbit hemorrhagic disease virus and drought in the geographic regions inhabited by rabbits.
INTRODUCTION
An experimental field study in Australia in 1950 to test whether myxoma virus (MYXV) could be used as a biological control for European rabbits (Oryctolagus cuniculus) accidentally initiated one of the largest and best-known natural experiments in host-pathogen coevolution. Detailed studies on the coevolution of the virus with its new host over the subsequent 60 years have done much to inform thinking on pathogen evolution in general and the evolution of virulence in particular (1–5).
MYXV is a poxvirus of the Leporipoxvirus genus (family Poxviridae; subfamily Chordopoxvirinae). The natural host of the virus released in Australia is the South American tapeti, or forest rabbit (Sylvilagus minensis [6], also called Sylvilagus brasiliensis [7, 8]). In the tapeti, MYXV causes a localized, cutaneous fibroma at the inoculation site (6). Virus is passively transmitted on the mouthparts of biting insects, such as mosquitoes or fleas, probing through the virus-rich fibroma. In general, MYXV does not appear to cause significant clinical disease in the natural host. However, when European rabbits (Oryctolagus cuniculus) are infected, the virus spreads systemically, overwhelming the immune system, to cause the lethal, disseminated disease myxomatosis (9). Related viruses, rabbit fibroma virus and Californian MYXV, occur in separate Sylvilagus spp. in North America.
Following its experimental release, MYXV spread across large swathes of rabbit-infested southeastern Australia (10). The released virus, originally isolated in Brazil and later termed the standard laboratory strain (SLS), had a case fatality rate (CFR) estimated at 99.8% (11). However, there was rapid natural selection for slightly attenuated viruses which, by allowing longer survival of the infected rabbit, were more efficiently transmitted by the mosquito vectors (12, 13).
To study the evolution of virulence, field isolates of MYXV were classified into five virulence grades based on the CFR and average survival time (AST) in small groups of laboratory rabbits which were the same species and had essentially the same outcomes following infection as wild European rabbits. Moderately attenuated grade 3 strains with CFRs of 70% to 95% predominated in the field in the decades following the release. More attenuated grade 4 viruses were less common, while highly virulent grade 1 and 2 viruses became relatively rare, despite ongoing releases of virulent virus. Highly attenuated grade 5 viruses with CFRs of <50% were never common (2, 13–16). At the same time, there was a strong selection for resistance to myxomatosis in the wild rabbit population. The short generation interval of rabbits meant that this evolutionary process could be measured in real time by challenging successive generations of rabbits with virus of known virulence. At one study site, CFRs dropped from 90% to 26% over a 7-year period (17, 18).
The deliberate release in France of a separate Brazilian isolate of MYXV in 1952 led to the establishment and spread of MYXV in Europe. Later termed the Lausanne (Lu) strain, this virus was shown to be more virulent than SLS when tested in genetically resistant rabbits (2, 19). However, although the starting virus, ecological conditions, and vectors were different, the evolutionary outcomes in Europe at the phenotypic level were remarkably similar to those in Australia, with the emergence of moderately attenuated strains of virus and selection for genetic resistance in the wild rabbit population (16). The Lu sequence is the type sequence for MYXV. This virus was originally isolated in 1949 and possesses a genome of 161,777 bp of linear double-stranded DNA (dsDNA) with inverted terminal repeats (TIRs) of 11,577 bp and closed single-strand hairpin loops at the termini. The genome encodes 158 unique open reading frames (ORFs), 12 of which are duplicated in the TIRs (20, 21).
The virus originally released in Australia, SLS, is believed to be derived from an isolate made in Brazil in about 1910 (22). A sample of this virus was later transferred to the Rockefeller Institute (23), where it was maintained by passage in laboratory rabbits. MYXV used for testing in Australia and Britain was obtained from the Rockefeller Institute in the 1930s (24). SLS has a genome of 161,763 bp with the same organization as that of Lu, but with single-base indels disrupting two ORFs, M005L/R in the TIRs and M083L (4). Although the Lu strain of MYXV was widely released in Australia for tactical rabbit control from the 1970s to 1990s, there is no evidence that this virus has established (25, 26), and all of the viruses that we have sequenced are derived from SLS.
Detailed molecular studies, in particular gene knockout experiments using the Lu strain, have revealed much about the way in which MYXV subverts and overwhelms the immune system of the European rabbit (9, 27). In many cases, deletion of a single virulence gene is sufficient to substantially attenuate the virus. Sequencing studies of field isolates of MYXV from both the Australian and European radiations have shown that reading frame disruptions in key virulence genes are not uncommon (4, 28–30). However, characterization of some of these viruses has revealed that disruption to virulence genes can be compatible with high virulence and that recent viruses from both radiations exhibit a novel immune-suppressive phenotype in laboratory rabbits (5, 30). In addition, reverse genetics to correct potentially attenuating mutations or insert a potentially attenuating mutation did not result in phenotypic changes (31). This has made it difficult to define the molecular basis of the initial attenuation or to define the molecular drivers of the ongoing virus evolution in wild rabbit populations, potentially in competition with continuing selection for resistance to myxomatosis.
Phylogenetic analyses of almost 50 full-genome sequences from Australian and European viruses have demonstrated that despite the similarity in evolutionary outcomes on the two continents, there were almost no shared mutations across the Australian and UK sequences and, hence, no convergent evolution at the genomic level (30). However, the evolutionary rates for the two radiations were both remarkably similar, at approximately 1 × 10−5 nucleotide substitutions per site per year. Virus evolution on both continents was also strongly clock-like, manifest in a consistent relationship between the time of sampling and the genetic distance (30).
Previous genomic analyses for the Australian viruses only included viruses sampled from 1950 to 1999, with most samples being obtained between 1990 and 1996. However, the ecological environment facing rabbits in Australia changed dramatically in the mid-1990s with the release of a second biocontrol agent, rabbit hemorrhagic disease virus (RHDV), that greatly reduced the size of the rabbit population (32–35). To determine whether MYXV evolution differed markedly before and after the emergence of RHDV, we sequenced a large set MYXV isolates sampled between 2008 and 2016, extending previous studies by nearly 20 years and considerably expanding the geographic coverage. In addition, to investigate possible variation in the progenitor virus, we sequenced SLS stocks prepared in 1949 and 1956.
(This article was submitted to an online preprint archive [36].)
RESULTS
Phylogenetic analysis of Australian MYXV.We sequenced and analyzed 50 field strains of MYXV sampled in southeastern Australia. Of these, 49 were collected during the period from 2008 to 2016. One archival sample from 1990 was also sequenced (Table 1). In addition, we resequenced the SLS progenitor and an SLS-like virus from archival samples and resequenced regions of the Glenfield virus (GV) strain (originally sampled in February 1951 during the first epizootic). The sequencing coverage for each virus ranged from 360 to 1,540 times. The geographic locations of the viruses sampled in this study and previously, as well as their lineages (see below), are shown in Fig. 1. The complete genome of each virus was obtained and aligned with the Lu reference genome covering nucleotide positions 1 to 161,777.
Australian field isolates of MYXV sequenced as part of this study
Sampling locations in Australia of the MYXV isolates sequenced in this study and previously. The sampling locations are color coded by phylogenetic lineage (see Results and Fig. 2). The different Australian states and territories are marked.
Our phylogenetic analysis of all complete MYXV genomes from Australia (n = 78) revealed that those viruses sampled between 1990 and 2016 can be placed into three main lineages, denoted here as lineages (clades) A, B, and C, although only lineage C received strong bootstrap support (Fig. 2). A single virus, Meby, isolated in Tasmania in 1991, could not be placed in these three lineages and fell in a more basal position. Lineage A contains viruses sampled between 1990 and 2016, including a small cluster of five viruses isolated between 2013 and 2016 (Fig. 2). Similarly, lineage B contains viruses sampled from 1993 to 2016 and from a variety of geographic and climatic areas, including semiarid South Australia (Port Augusta, Olympic Dam, Flinders Ranges, Wilpena, ABC Range) and more temperate, higher-rainfall areas in South Australia, New South Wales (NSW), and close to Canberra in the Australian Capital Territory (Fig. 1 and 2). This lineage also includes previously sequenced viruses isolated from the Canberra region and western Queensland in the 1990s. Interestingly, all but four viruses in clade B, including all the recent isolates, have a 1,635-nucleotide (nt) duplication encompassing the M156R and M154L genes and part of the M153R genes normally found at the right-hand (RH) TIR boundary that has been duplicated at the left-hand (LH) TIR boundary; this is effectively an expansion of the TIRs, with the loss of 923 nt of the 3′ end of the M009L gene at the LH TIR boundary (Fig. 2). This duplication was first seen in viruses isolated from the Canberra district in 1995 (26, 28).
Evolutionary history of MYXV in Australia. The tree depicts the phylogenetic relationships among all 78 complete genome sequences of MYXV sequenced to date. All viruses sampled from 1990 have been color coded and, with exception of the Meby strain sample in Tasmania in 1991, fall into three main lineages (denoted A, B, and C), with their range of sampling times shown. Two branches—denoted x and y—associated with the emergence of the rapidly evolving lineage C are indicated. All viruses within lineages B and C, aside from the four marked, possess a duplication in M156R, M154L, and part of M153R. The tree is rooted by the European Lausanne strain, and all branches are scaled according to the number of nucleotide substitutions per site (subs/site). All bootstrap values of >70% are shown. (Inset) Regression of root-to-tip genetic distances against the year of sampling for all available Australian isolates of MYXV. Each point is color coded according to the phylogenetic lineage. The recently sampled lineage C, which has experienced an elevation in evolutionary rate, is marked.
Finally, lineage C includes viruses isolated from different geographic and climatic areas of Victoria (Hattah, Hoppers Crossing, Wonga Park), southern NSW (Deniliquin), and a broad geographic and climatic region of South Australia from the far north desert country (Pukatja) to the higher-rainfall Adelaide Hills and Mt. Gambier in the south. This lineage is notable in that it only contains viruses sampled from 2012 onwards and that all viruses contain the TIR expansion seen in clade B. More dramatically, individual tree branches within lineage C are longer than those in the other two lineages and exhibit a marked lack of temporal structure; that is, many of the earlier sampled viruses (from 2012) fall no closer to the root of the tree than those sampled in 2016. This includes viruses sampled between 2012 and 2016 at the same site: Turretfield in South Australia. That viruses sampled in 2012 occupy phylogenetic positions both at the root and toward the tips of this lineage suggests that most of the divergence in this lineage, and hence the rapid evolutionary event, occurred prior to 2012.
To assess the extent of temporal structure in the Australian MYXV sequences in a more quantitative manner, we performed a regression of root-to-tip genetic distances against the year of sampling. This revealed that all MYXV sequences, with the exception of those assigned to lineage C, exhibited a relatively strong rate constancy, as previously observed in MYXV in both Australia and Europe (30). In contrast, most of those viruses assigned to lineage C had elevated rates, such that they consistently fell above the regression line (Fig. 2, inset). Across the Australian viruses as a whole, the R2 value of the relationship between root-to-tip genetic distances against the year of sampling was 0.58, whereas values of between 0.93 and 0.98 were observed in previous studies that did not include lineage C. Importantly, the viruses assigned to lineage C were sequenced at the same time as those from the other lineages, indicating that their divergent position is not a laboratory artifact.
Mutations associated with the rapidly evolving MYXV lineage C.To determine what processes might be associated with the rate elevation in lineage C, we mapped all nucleotide and amino acid substitutions onto the maximum likelihood (ML) phylogenetic tree and examined in detail those that occurred early in the history of lineage C. Specifically, we focused on (i) the branch directly leading to lineage C immediately prior to the divergence of the Hoppers Crossing virus (designated branch x in Fig. 2) and (ii) the first (i.e., next) internal branch in lineage C (designated branch y in Fig. 2).
A total of 10 nucleotide substitutions were reconstructed as falling on branch x, 9 of which were nonsynonymous. Such a relatively high frequency of nonsynonymous substitutions is compatible with either a relaxation of selective constraints or, perhaps more likely, the occurrence of positive selection on this lineage. These nonsynonymous substitutions included an amino acid change in the M010L gene (A66V). M010 is a homologue of cellular epidermal growth factor/transforming growth factor alpha and is critical for virulence (37, 38). Interestingly, this is the only amino acid substitution in M010 that we have seen in Australian isolates, although a valine is also found at this position in the Californian MYXV strain MSW and rabbit fibroma virus. Two nonsynonymous mutations were also seen in M036L, which encodes a protein orthologous to the orthopoxvirus vaccinia virus (VACV) O1. O1 potentiates extracellular signal-regulated kinase 1/2 (ERK1/2) signaling stimulated by the VACV epidermal growth factor homologue (39). There are no data on the role (if any) of M036 in MYXV infection and whether it potentiates the activity of M010. However, mutations in the M036L gene are relatively common, and disruption of the gene in field isolates of MYXV is compatible with high virulence (5, 25, 30).
In addition to the novel mutations found on this branch, there were two reversals of nonsynonymous mutations in M014L and M134R that occurred on the root branch leading to all the sequences from 1990 onwards. The role of M014 is undefined, but it likely functions as part of an E3 ubiquitin (Ub) ligase complex (40). M134 is an orthologue of variola virus (VARV) B22, which has been implicated in suppression of T cells in orthopoxvirus infections (41).
A total of 19 substitutions were constructed to fall on branch y (Fig. 2); 10 of these were nonsynonymous, and 6 represented a reversal of mutations occurring on branches leading to all the post-1990 sequences. Perhaps of most interest is an additional single nucleotide deletion of a G nucleotide at nt 31 in the 5′ end of the M005L ORF (and in the corresponding position in M005R) that corrects a G insert at this position found in SLS and in all subsequent Australian virus sequences previously sequenced. The G insert disrupts the M005L/R ORF with stop codons after codon 78 and an altered amino acid sequence after codon 11 (Fig. 3). Targeted disruption of the M005L/R gene in Lu caused severe attenuation of the virus, with case fatality rates being reduced from 100% to 0 (42). All of the viruses within rapidly evolving lineage C have an intact M005L/R ORF, with the exception of the Hoppers Crossing virus.
Reading frame correction in M005L/R. (A) Partial nucleotide sequences of M005R for Lu, SLS, and representative group A, group B, and group C viruses are shown. The C insertion at nt 31 in SLS and all other Australian viruses in group A (Lameroo) and group B (Acton 08-2014 [Acton 8-14]) (shaded) is corrected by a C deletion in the group C viruses (Pukatja), except for the Hoppers Crossing virus. The early stop codon at nt 235 to 237 is underlined. Acton 08-2014 has a T > C mutation at nt 133 that inserts an earlier in-frame stop codon (underlined and shaded). There is an alternate ATG (underlined) in SLS upstream of the C insertion which would correct the reading frame. However, we have no evidence that this is used. (B) Partial amino acid sequences of M005 for the same viruses showing an altered amino acid sequence after amino acid 11 (underlined), an early stop after amino acid 78 for SLS, and Lameroo and an early stop after amino acid 44 for Acton 08-2014. The full-length M005 protein is 483 amino acid residues. Virus names have been abbreviated for ease of presentation.
Recombination in Australian MYXV.Recombination may be relatively common in poxviruses (43) and has been observed between MYXV and rabbit fibroma virus (44, 45) and between South American and Californian MYXV (46, 47). Despite this, we previously found no evidence for recombination in either the Australian or European MYXV (4, 28, 30).
An analysis of possible recombination events in our expanded data for Australian MYXV revealed a complex pattern of phylogenetic movement compatible with the action of recombination, although the location of precise recombination breakpoints is difficult to determine. An analysis of recombination events using a suite of methods within the RDP package failed to identify either consistent (i.e., statistically significant) recombination events or recombination breakpoints. For example, although some analyses suggest that the Hattah (2012), Coomandook (2013), and Adelaide Hills (2012) viruses are recombinant, this was not confirmed in more detailed phylogenetic analyses because there was insufficient bootstrap support to confirm phylogenetic incongruence (data not shown).
In an additional attempt to reveal recombinants, we inferred phylogenetic trees for nonoverlapping 25-kb regions of the MYXV genome (Fig. 4). This revealed clear, although complex, patterns of topological movement, and in most cases it was impossible to statistically prove recombination (i.e., prove with strong bootstrap support for incongruent trees) because of a lack of phylogenetic resolution due to low levels of sequence diversity. Indeed, the only case in which the differences in tree topology were supported by high levels of bootstrap support concerned a cluster of three viruses from Euchareena sampled in NSW in 2012 and 2013 and a singleton comprising a virus from Port Augusta in South Australia sampled in 2014. These viruses fell into lineage B in phylogenies of the whole-virus genome, as well as that of the region from 25 to 50 kb but clustered with some lineage C viruses in the region from 75 to 100 kb, particularly four viruses sampled at Turretfield (South Australia) in 2016 (Fig. 4).
Phylogenetic evidence for recombination in Australian MYXV. (A) The tree for the aligned genomic region from 25 to 50 kb. (B) The tree for the region from 75 to 100 kb. All taxa are color coded according to their whole-genome lineage (Fig. 2). Putative recombinant strains are shaded, with the relevant sequence names shown in bold. All nodes with bootstrap support values of >70% are shown. All trees are drawn to a scale of the number of nucleotide substitution per site and rooted by the Lausanne sequence.
In addition, as noted below, some individual mutations exhibited phylogenetic distributions that were compatible with the occurrence of recombination, although this could not be demonstrated in a statistically robust manner.
Gene-specific mutations in Australian MYXV analysis.Across all Australian MYXV isolates sequenced to date, 13 genes had no mutations and 16 had only synonymous mutations. These conserved genes largely encode structural proteins, fusion/entry proteins, or those involved in genome replication or transcription. The exceptions include M128L, which encodes an immunomodulatory protein important for virulence (48), and M123R, an orthologue of VACV A35R that encodes an immunomodulatory protein (49, 50). Multiple genes had in-frame indels, reading frame-disrupting indels, or stop codons (Table 2).
Insertions, deletions, disruptions, and repairs in genes from new MYXV sequencesa
Mutations in proteins involved in DNA replication and transcription.We examined whether high evolutionary rates could have been due to a mutator phenotype or possibly more rapid replication. The MYXV DNA polymerase (DNA pol) complex consists of four proteins: M034 (DNA pol; VACV E9), M079 (DNA processivity/uracil DNA glycosylase; VACV D4), M080 (helicase/primase; VACV D5), and M111 (DNA pol processivity; VACV A20) (51). M034 and M079 had no amino acid substitutions in any virus. A single amino acid substitution was observed in M080 in one virus (Bluegums), and another amino acid change, R430C in M111, was identified in some of the rapidly evolving viruses from lineage C (Wonga Park, Adelaide Hills, Turretfield 2012 isolates, Mt. Gambier, Pukatja, Quorn, Sandy Creek). Published mutational studies of VACV A20 do not provide any information on the phenotypic impact of this substitution (52–54), although the cysteine at this position is also found in rabbit fibroma virus.
In addition to the polymerase complex, M112 (Holliday junction resolvase; VACV A22), M133 (DNA ligase; VACV A50), and M074 (DNA topoisomerase; VACV H6R) are required for DNA replication. None of the genes for these proteins contained any nonsynonymous substitutions. Proteins involved in generating nucleotide precursors, M061 (thymidine kinase), M015 (ribonucleotide reductase small subunit), and M012 (dUTPase), had no consistent mutations in lineage C viruses.
Of 27 genes encoding proteins involved in virus transcription, most had only synonymous mutations, mutations that occurred only in single viruses, or mutations common to most of the Australian viruses. However, there were two mutations in M044, the viral NPH II protein (nucleoside triphosphatase/RNA helicase), that occurred uniquely in most lineage C viruses (R118C, P478S). Of note was the finding that a number of viruses (Wonga Park, Hattah, Adelaide Hills, Turretfield 2012 and 2016 isolates, Coomandook, Deniliquin, Mt. Gambier, Pukatja, and Sandy Creek) had P478S, although the proline at this site (within the helicase C-terminal domain) was otherwise completely conserved across the chordopoxviruses.
Interestingly, there was a widespread mutation in the C terminus of the RNA polymerase subunit encoded by M114R (VACV A24) in Australian MYXV isolates (P1147H). Many of the viruses with this M114 mutation also shared a mutation in M156 (L71P), the MYXV member of the poxvirus K3 family of protein kinase R (PKR) inhibitors, that abolishes PKR inhibition (55). A T1120M mutation in the VACV RNA polymerase subunit A24 was able to confer resistance to PKR in an artificial evolution system (50).
Analysis of MYXV promoters.We examined mutations up to 100 nt 5′ to ATG start codons to determine whether there were lineage-specific differences in promoter sequences. Since only a small number of MYXV genes have had their promoters mapped, promoters have necessarily been inferred by homology with conserved chordopoxvirus early, intermediate, and late promoter sequences (56). Recent studies with VACV have demonstrated that transcriptional regulation may be complex (57–59). The majority (but not all) of the viruses in lineage C had a T deletion at the 3′ end of the M156R gene (Table 2), which forms part of the upstream T-rich region of the M008.1L/R late promoter. This mutation was also seen in three lineage B viruses and might be expected to decrease transcription from this promoter, based on exhaustive mutational studies of poxvirus late promoters (60). An additional A at the 5′ end of the M040L E promoter sequence was present in most recent Australian isolates in lineages A and B but not lineage C, which might be expected to have a minor increase in transcription efficacy, based on mutational studies of early promoters (61). Most other mutations were considered unlikely to have an effect on promoters or occurred only in single sequences.
Mutations causing loss of the M153R ORF are common and persist in the field.The M153R gene encodes a 206-amino-acid RING-CH E3 ubiquitin ligase that downregulates major histocompatibility class I, CD4, Fas/CD95, and ALCAM/CD166 from the surface of infected cells in vitro (62–65). Targeted deletion of this gene in Lu causes marked attenuation of the virus (62). However, mutations that disrupt the ORF and that are predicted to lead to a loss of function were common in Australian MXYV isolates (Fig. 5).
Mutations in M153R leading to reading frame disruptions. (A) Partial nucleotide sequence of M153R from SLS, Bluegums (BGM), Monarto Zoo (MON), and ABC Range (ABC). The early stop codon formed by 3 nucleotide substitutions in Bluegums is underlined and shaded. The A deletion in Monarto Zoo is shaded, and the G insertion in ABC Range is shaded. (B) Amino acid sequences of M153 for the same viruses described in the legend to panel A. The RING-CH E2 binding site at the N terminus is shaded, the two transmembrane domains are underlined, and the C-terminal conserved acidic domain is shaded. Lowercase amino acids in Monarto Zoo and ABC Range indicate an altered amino acid sequence following the indel. Virus names have been abbreviated for ease of presentation.
Our previous sampling of field isolates of Australian MYXV did not allow us to determine whether any of the viruses with disrupted M153R genes left descendants. However, with increased sampling we can now identify multiple viruses with M153R disruptions and their descendants. The first of these was initially seen in a lineage B virus isolated in Canberra in 2008 (Acton 07-2008) and in later isolates sampled from the same site (Acton 04-2014, Acton 08-2014), approximately 35 km to the north (Murrumbateman 2014 and Bluegums 2015), and 10 km to the southeast (Symonston 2016). This is a complex mutation in which three nucleotide substitutions within a 4-nucleotide span, AAAC > GTAG, form an in-frame TAG stop codon after codon 11 (Fig. 5). The next in-frame ATG is at nt 275. The identical mutation was found in a lineage C virus isolated at Deniliquin in 2015, some 400 km west of Canberra. The complex nature of this mutation makes it unlikely to have occurred more than once, and its presence in viruses in different lineages may therefore reflect the occurrence of recombination.
In the second example, a G insertion at nt 330, causing disruption of the M153R ORF with an early stop after amino acid 124, was found in viruses isolated from Wilpena and the ABC Range in the Flinders Ranges of outback South Australia in November and December 2014 (Fig. 5). The same mutation was present in related viruses isolated 12 months later in the Flinders Ranges and in August 2015 at Olympic Dam, some 100 km to the west. While we cannot exclude the possibility that this insertion occurred multiple times independently, it seems more likely that these closely related viruses were descendants of a single common ancestral strain that contained the mutation.
Sequences of the Australian progenitor SLS virus.The original epizootic of myxomatosis in 1950 stemmed from wild rabbits inoculated with virus later termed the standard laboratory strain (SLS). In the following years, many releases of this virus were made by inoculating rabbits. It is possible that polymorphisms in stocks contributed to diversity in the field. For example, the nucleotides at positions 2576 and 149127 in the stock of SLS that we originally sequenced were not found in any other Australian viruses. This suggests that mutations at these positions may have occurred during passage in rabbits or cell culture (4) or existed as polymorphisms in stocks of SLS used for release. To test this, we grew and sequenced virus from stocks originally obtained from Frank Fenner and labeled by him as “SLS (1956)” and “Myxoma virus 1949.”
There are a number of sequence differences between our original SLS and the 1956 stock and 1949 viruses (Table 3). In the case of nt positions 2756 and 149127, the SLS 1956 stock had the same sequence as the field isolates, suggesting that these were the original sequences at these positions. In addition, the 1956 virus had two single nucleotide insertions which were seen in some early isolates: the G insertion at position 122397 was present in the GV strain (which was isolated in February 1951 but which was more virulent than SLS, despite the gene disruptions), and the A insertion at position 131587 was present in the Ur strain (1953). Hence, the originally released SLS may have comprised a mixture of sequences. The provenance of the 1949 virus is not known. It is clearly closely related to SLS and distinct from Lu. Fenner obtained and tested virus from the Rockefeller Institute (13) from which SLS was originally derived. If this is that virus, then it appears that SLS has diverged during passage. In addition, using a new stock of virus, we confirmed the gene disruptions in M014L, M130R, and M153R originally identified in the high-virulence GV strain (4).
Sequence variation in SLS stocks
DISCUSSION
Previous studies of MYXV evolution had shown that the evolutionary rate in this virus since its release in Australia and Europe was comparatively high for a DNA virus at ∼1 × 10−5 nucleotide substitutions per site per year (4, 30) and extremely clock-like in both continents (30). While this still appears to be the case for more recent isolates assigned to lineages A and B, it is clear that the isolates in a third and large group of viruses, which were designated lineage C viruses and which were collected from diverse locations across southeastern Australia, have a much greater genetic diversity and a breakdown of the clock-like structure, compatible with a greatly elevated rate of evolutionary change. The temporal structure in lineage C, with the earliest sampled viruses from 2012 falling at diverse phylogenetic locations, including both close to the root and near the tips, suggests that the full diversity in this lineage, and hence the rate elevation, were established at some point between 1996 (i.e., the most recent date of the majority of our last sample of sequences) and 2012.
What, then, is the most likely cause of this punctuated evolutionary event? Sequence comparison did not reveal any convincing evidence for mutations in genes involved in DNA replication or transcription that might lead to a more error-prone replication or more cycles of replication, which could have conceivably led to an increase in the rate, although this clearly merits detailed experimental analysis. Similarly, while there is some evidence for recombination in the recent isolates of MYXV, including some viruses within lineage C, its impact appears to be localized to specific viruses and particular genomic locations. We therefore suspect that recombination is unlikely to be responsible for the punctuated evolution of lineage C.
A change in selection pressures during the period from 1996 to 2012 may be the most likely explanation for the emergence of the divergent viruses associated with lineage C. Two major changes in rabbit populations occurred during this period, and these could conceivably have influenced the transmission of MYXV and, hence, altered selection pressures. From 1995, a novel viral pathogen, RHDV, spread through Australia, causing rabbit population crashes in many areas and reducing the number of MYXV-susceptible hosts (32–35, 66–69). In addition, a major, prolonged drought occurred over southeastern Australia between 1996 and 2009, putting additional pressure on rabbit populations. A possible third influence is the establishment of the Spanish rabbit flea (Xenopsylla cunicularis) as a vector for MYXV in the arid areas of southeastern Australia, although there is little information on its impact. Decreased mosquito numbers due to dry conditions may have meant that fleas became the major vectors for the transmission of MYXV, in turn making transmission much more localized or making it completely fail in regions where fleas were absent. Changes in the timing of epizootics of myxomatosis were also seen (67). Smaller, fragmented rabbit populations may have favored viruses with longer average persistence in the individual rabbit, and local extinctions of MYXV were likely commonplace. As rabbit populations recovered from both the initial impact of rabbit hemorrhagic disease (70) and the drought, any viruses better able to persist, perhaps by adaptation to more efficient transmission, would have had a selective advantage.
Within lineage C, it is tempting to speculate that the reversion of the M005L/R gene mutation to generate an intact ORF provided the phenotype conferring that selective advantage. The M005 protein has been shown to inhibit apoptosis in MYXV-infected lymphocytes in European rabbits and to have host range functions in human cancer cells and is a key virulence gene in Lu (42, 71, 72). Mutational disruption of M005 has been shown to reduce virulence in the Lu genetic background with no development of secondary cutaneous lesions (42). Previous Australian isolates do not appear to have expressed M005, despite many having high virulence. The role of virulence in transmission in low-density rabbit populations with high genetic resistance is likely complex, and it is possible that mutations acting on other phenotypic characteristics, such as lesion morphology, may be important. The genetic determinants of virulence in the field are clearly complex and likely multifactorial, and wild-type virulence is not easily explained by the so-called virulence determinants identified in gene knockout studies.
Aside from the potential selection pressures acting on lineage C, there is the surprising success of viruses with the duplication of M154L and M156R. Whether there is a selective advantage in duplication of two potential immunomodulatory genes or loss of the M009L gene is not clear. M154 has homology to VACV M2, which inhibits NF-κB activity (73), while M156 inhibits PKR (55, 74). In artificial evolution experiments using VACV, duplication of PKR inhibitors, such as K3L, occurred quickly and provided a selective advantage (75, 76). Interestingly, a mutation (L71P) in M156R that was present in all of the viruses in lineage B and in the duplicated M156R genes has been demonstrated to inactivate PKR inhibition (55). Despite this loss of function, some viruses with this mutation have high virulence in laboratory rabbits (5, 25). In addition, reverse engineering a virulent virus to have the 71P mutation did not alter virulence (31). However, most viruses in lineage C, all of which had the duplication, did not have the L71P mutation, indicating a reversion at this position. The majority of these viruses also had a T deletion in the 3′ end of the M156R gene, causing readthrough of the stop codon and the addition of two amino acids (E76 and G77) at the C terminus.
The M009L gene encodes a 509-amino-acid protein that is likely part of an E3 ubiquitin ligase complex (40). M009L has suffered multiple disruptions during the evolution of MYXV in Australia (Table 2) (4). The M009L ORF was also lost in the Californian MSW strain of MYXV and in the Kazza strain of rabbit fibroma virus (47, 56). Only one group of viruses sequenced here had an intact M009L ORF—Kangarilla, Brooklyn Park, Lameroo, Echunga, and Monarto in lineage A—as did the previously sequenced Bendigo virus from 1992 (also lineage A).
In sum, during the period from 1996 to 2012, viruses that had the M154L/M156R gene duplication and that comprised the newly described lineages B and C became dominant in our sampling over a wide geographic area of southeastern Australia. Within this group is a highly divergent and rapidly evolving virus lineage in which the normally constant rate of molecular evolution has been disrupted. Although we are uncertain what evolutionary forces are responsible for this punctuated evolutionary event, we suggest that it may reflect a combination of the combined impact of RHDV and a large-scale drought, both of which would have severely reduced rabbit numbers and, in turn, changed the selective environment facing MYXV. Similarly, while it is currently unclear what mutations may have been most impacted by these selective events, the mutations in M005L/R and M156R, which are of plausible phenotypic importance, should clearly be a starting point in further experimental investigation.
MATERIALS AND METHODS
Field samples for virus isolation.Samples for virus isolation were submitted either as tissues (usually eyelids) or as conjunctival swabs from rabbits with clinical signs of myxomatosis found dead or collected for other purposes. All samples were processed and viruses were isolated on RK13 cell monolayers as previously described (30).
Archival samples of SLS and GV.Small glass straws containing freeze-dried aliquots of SLS produced in rabbits were originally obtained from the late Frank Fenner (Australian National University) in 1988. A vial labeled “myxoma virus 1949” and containing freeze-dried material was also present from the same source. Virus from these vials was amplified in RK13 cells. Vials of freeze-dried Glenfield virus (GV) were produced by the Commonwealth Serum Laboratories in 1962 as homogenates of tissue from infected rabbits.
DNA preparation.Virus was concentrated from infected RK13 cells by osmotically bursting the cells, centrifugation to remove nuclei and other cell debris, and DNase and RNase digestion, followed by polyethylene glycol precipitation (30, 77, 78). Viral DNA was prepared using a MasterPure DNA extraction kit (Epicorp).
PCR analysis and sequencing of disrupted genes in GV.The primer sequences used were as follows: for M153R, forward primer AATGGAGGTGTTATAAACGCGACTGCCACG-3′ and reverse primer CGATCTTTGTTATAGACAAACGAGATACCT; for M014L, forward primer GCCAGGTCTATTCTGTCGATC and reverse primer CGTCATTTACGTCTTGGGAGGAGTCTCGTAC; and for M130R, forward primer TATATTATACACGGCCTTATTCGACGGCG and reverse primer CGTCGTACGGAAGGTGACTGTCTACGTTAA.
The sequences of both strands of the PCR-amplified DNA were determined by Sanger sequencing.
MYXV genome sequencing.Purified viral DNA was first quantified using a Qubit dsDNA assay (Invitrogen), before being diluted to 0.2 ng/µl with molecular-grade water. Sequence libraries were then prepared from the diluted viral DNA using a Nextera XT sample preparation kit (Illumina). Individual libraries were quantified using a Kapa library quantification kit (Kapa Biosystems) before manual normalization. The pooled libraries were then sequenced on a 500-cycle (v2) Illumina MiSeq sequencer (250-nt paired end reads) at the Ramaciotti Centre for Genomics, Sydney, Australia, with between 12 and 24 samples being indexed in each run.
Genome assembly.Sequence quality control and genome assembly were performed using CLC Genomics (v8) software (Qiagen). First, the demultiplexed sequence reads were imported and quality trimmed to ensure removal of adapter sequences and low-quality bases (phred score < 20). Next, overlapping sequence reads were merged. The full set of trimmed sequence reads was then assembled de novo with a word (k-mer) size of 35. This typically yielded a large contig corresponding to the core genome (∼138 kb) and a small contig corresponding to the TIR (∼11 kb). The TIR contig was duplicated, and one copy was reverse complemented before manually assembling each onto the end of the core genome (where sequence overlap was present) to generate a draft genome. Following this procedure, the full set of trimmed sequence reads was remapped back to check for polymorphisms, misalignments, and the evenness of coverage, with the final majority consensus taken as the final genome assembly.
Gene annotation.Genes were annotated as in the Lu reference (Brazil/Campinas-1949; GenBank accession number NC_001132) (20) and SLS (Brazil-1910; GenBank accession number JX565574) (4) sequences. Gene numbers are italicized, with the direction of transcription indicated as left (L) or right (R) (for example, M010L). Genes in the terminal inverted repeats are indicated by L/R. Proteins are indicated by the gene number without the transcription direction (for example, M010). In the case of the M156R gene, we revised the annotation to begin at a start codon downstream to the original annotation, such that the protein is 75 instead of 102 amino acids (47, 55).
Evolutionary analysis.The 50 MYXV genome sequences determined here were combined with the sequences of 25 previously sequenced viruses from the Australian outbreak. Details on these previously sequenced viruses, including their place and date of sampling, can be found elsewhere (4, 28). The European Lausanne sequence was used as an outgroup to root the tree. These sequences were initially aligned using the MAFTT program (79) and adjusted manually, resulting in a final sequence alignment data set of 78 MYXV sequences 163,711 bp in length and covering the years 1949 to 2016. The final alignment is available from the authors on request.
A phylogenetic tree of the aligned sequences was estimated using the maximum likelihood (ML) procedure available in the PhyML package (80). This analysis utilized the GTR+Γ4 model of nucleotide substitution and NNI+SPR (nearest-neighbor interchange and subtree pruning and regrafting) branch swapping. A bootstrap resampling procedure (1,000 replications) was used to assess the statistical robustness of individual nodes on the MYXV phylogeny. To determine and visualize the extent of clock-like structures in these MYXV data, we performed a regression of root-to-tip genetic distances against the year of sampling on the ML tree using the TempEst software tool (81).
To test for the presence of recombination in the Australian MYXV data, we first performed a screen for recombination using the RDP, Genecov, and Bootscan methods (with default settings) available within the RDP4 package (82). In addition, we estimated ML phylogenetic trees for nonoverlapping 25-kb sections of all the Australian MXYV sequences sampled here using the same procedure described above. A visual inspection for phylogenetic incongruence supported by high numbers of bootstrap replications (>70%) was then used to identify putative recombination events.
Accession number(s).All consensus MYXV genome sequences generated here have been submitted to GenBank and assigned accession numbers MK388093 to MK388144. The raw sequence reads have been submitted to SRA under BioProject accession number PRJNA513218. The BioSample accession numbers for each virus are SRR8402032 to SRR8402083.
ACKNOWLEDGMENTS
Rabbit samples were kindly provided by Tarnya Cox and Pam Whitely. Karen Glover provided technical support.
E.C.H. is funded by an ARC Australian Laureate Fellowship (FL170100022). This work was funded in part by the National Institute of Allergy and Infectious Diseases (grant R01AI093804).
FOOTNOTES
- Received 7 November 2018.
- Accepted 23 January 2019.
- Accepted manuscript posted online 6 February 2019.
- Copyright © 2019 American Society for Microbiology.
REFERENCES
- 1.↵
- 2.↵
- 3.↵
- 4.↵
- 5.↵
- 6.↵
- 7.↵
- 8.↵
- 9.↵
- 10.↵
- 11.↵
- 12.↵
- 13.↵
- 14.↵
- 15.↵
- 16.↵
- 17.↵
- 18.↵
- 19.↵
- 20.↵
- 21.↵
- 22.↵
- 23.↵
- 24.↵
- 25.↵
- 26.↵
- 27.↵
- 28.↵
- 29.↵
- 30.↵
- 31.↵
- 32.↵
- 33.↵
- 34.↵
- 35.↵
- 36.↵
- 37.↵
- 38.↵
- 39.↵
- 40.↵
- 41.↵
- 42.↵
- 43.↵
- 44.↵
- 45.↵
- 46.↵
- 47.↵
- 48.↵
- 49.↵
- 50.↵
- 51.↵
- 52.↵
- 53.↵
- 54.↵
- 55.↵
- 56.↵
- 57.↵
- 58.↵
- 59.↵
- 60.↵
- 61.↵
- 62.↵
- 63.↵
- 64.↵
- 65.↵
- 66.↵
- 67.↵
- 68.↵
- 69.↵
- 70.↵
- 71.↵
- 72.↵
- 73.↵
- 74.↵
- 75.↵
- 76.↵
- 77.↵
- 78.↵
- 79.↵
- 80.↵
- 81.↵
- 82.↵