Resolving the Origin of Rabbit Hemorrhagic Disease Virus: Insights from an Investigation of the Viral Stocks Released in Australia

Abstract

To resolve the evolutionary history of rabbit hemorrhagic disease virus (RHDV), we performed a genomic analysis of the viral stocks imported and released as a biocontrol measure in Australia, as well as a global phylogenetic analysis. Importantly, conflicts were identified between the sequences determined here and those previously published that may have affected evolutionary rate estimates. By removing likely erroneous sequences, we show that RHDV emerged only shortly before its initial description in China.

FIG 1

SNPs in the original stocks of RHDV introduced into Australia. Panel A shows the genome positions (nucleotides) and frequencies (percent) of the SNPs in the original RHDV stocks. Trends in SNP frequency are indicated following passaging from the master stock, through an intermediate working stock, to the release strain with SNPs of increasing and decreasing frequencies plotted above and below the x axis, respectively. The genome positions on a gray background are the SNPs located in the major capsid protein (VP60). Panel B maps the three nonsynonymous SNPs in the capsid gene onto the published crystal structure of the RHDV VP60 protruding (P) domain dimer (Protein Data Bank code 4EGT). A side view is shown with the surface of one subunit (colored gray) in the background, the ribbon structure of the other subunit in the foreground, and the P1 and P2 subdomains colored blue and yellow, respectively. The table provides the codon sequences of all of the SNPs located in the capsid, with the nucleotide (nt) of each SNP in bold and the coding amino acid (aa) in brackets. The prototype sequence corresponds to the majority consensus sequence of the original master stock. Across both panels, shades of green and red indicate trends of increasing and decreasing SNP frequency, respectively.

FIG 2

Phylogenetic analysis of the original stocks of RHDV from Australia and New Zealand. Shown is a phylogeny of the complete genome sequences of 28 RHDV strains, including six original stocks of RHDV from Australia and New Zealand (green), as well as the previously published sequences of the purported founder virus, RHDV V351 (red). The gray dashed line shows the sequences compared here with the number of mismatches (nucleotides [nt]) and overall genetic difference (percent) indicated. The phylogeny was rooted with an early European RHDV strain (outgroup), and branch lengths are scaled according to the number of nucleotide substitutions per site. Support for individual nodes was estimated from 1,000 bootstrap replicates with asterisks marking nodes with ≥70% bootstrap support.

FIG 3

Estimates of the evolutionary rate of RHDV and the timing of its emergence. A set of 174 complete VP60 capsid sequences from globally circulating RHDV strains was used to estimate a maximum-likelihood phylogeny with PhyML. Panel A shows the root-to-tip genetic distances plotted against the sampling time with linear regression performed to examine the clock-like structure of the data. A Bayesian Markov chain Monte Carlo approach was then used to estimate the evolutionary rate (number of substitutions per site per year), shown in panel B, and the date of emergence (i.e., time to the most recent common ancestor [MRCA]), shown in panel C, for the same VP60 data set and a range of different models. This analysis incorporated both strict and relaxed (uncorrelated lognormal [UCLN]) molecular clocks, as well as constant-size (Con), exponential-growth (Exp), and Gaussian Markov random-field skyride (Sky) coalescent models. In panels B and C, the mean values are shown with black circles with the ranges reporting the 95% highest probability densities and the best-fit model, in this case, the relaxed clock (UCLN) with an exponential-growth population model (bold).