ABSTRACT
Human B19 erythrovirus is a ubiquitous viral pathogen, commonly infecting individuals before adulthood. As with all autonomous parvoviruses, its small single-stranded DNA genome is replicated with host cell machinery. While the mechanism of parvovirus genome replication has been studied in detail, the rate at which B19 virus evolves is unknown. By inferring the phylogenetic history and evolutionary dynamics of temporally sampled B19 sequences, we observed a surprisingly high rate of evolutionary change, at approximately 10−4 nucleotide substitutions per site per year. This rate is more typical of RNA viruses and suggests that high mutation rates are characteristic of the Parvoviridae.
Human B19 erythrovirus was first discovered in the serum of healthy blood donors (3) and has since been detected worldwide. Infection usually occurs in childhood, through respiratory droplets, and by age 15 approximately 50% of children have antibodies to the virus. Most childhood infections are asymptomatic, with erythema infectiosum, probably caused by the formation of immune complexes, the most frequent complication. In adults, however, infection often results in arthropathy and in some cases causes transient aplastic crisis. In immunocompromised individuals a persistent infection usually occurs, resulting in red cell aplasia, while transplacental transmission can lead to miscarriage or hydrops fetalis (22).
Human erythroviruses belong to the family Parvoviridae, whose single-stranded (ss) DNA genomes are among the smallest of all DNA viruses. B19 virus is ∼5.5 kb in length with two major open reading frames flanked by inverted terminal repeats, part of which form hairpin stems for priming replication through a double-stranded (ds) intermediate (2). The first open reading frame encodes the nonstructural (NS1) protein, while the second encodes the capsid proteins VP1 and VP2, which are in frame and colinear with the exception of an additional 227 amino acids at the N terminus of VP1. Amino acid variability is high in the VP1 unique region (19), which is surface exposed and the target site of neutralizing antibodies (8, 13).
Unlike large dsDNA viruses, all autonomous parvoviruses replicate with host cell machinery (10). A common assumption in studies of viral evolution is that DNA viruses have low rates of evolutionary change, near those of their hosts, as observed in the large dsDNA herpesviruses (7, 14) and the small dsDNA human papillomaviruses (1). However, it was recently observed that one group of small ssDNA viruses, the carnivore parvoviruses, have a rate of nucleotide substitution many orders of magnitude higher, at approximately 1 × 10−4 substitutions/site/year, that is within the range seen in RNA viruses (17). It is currently unclear whether this unexpectedly high rate is characteristic of the carnivore parvoviruses alone or whether it typifies the entire Parvoviridae family. To address this issue we estimated the rate of nucleotide substitution in the human virus B19, a distant relative of the carnivore parvoviruses.
B19 sequences (human erythrovirus genotype 1) were compiled and aligned, and isolate sampling dates were obtained through direct communication with the author or from published material. Only one sequence was taken in the case of a patient with multiple sampling times (I1) and from an outbreak with identical sequences (USA8). Two data sets were compiled: (i) 43 VP1 gene sequences (with nucleotides 1429 to 1461 and 2281 to 2346 removed because of incomplete sequencing), and (ii) 27 coding regions, beginning at position 39 of NS1 and running through the VP1 sequence described above (data not shown). Phylogenetic trees of these data were then inferred using the maximum likelihood (ML) method in PAUP* (18), employing the GTR+I+Γ4 model of nucleotide substitution, and rooted with the oldest sequence. Bootstrap support values were estimated using 1,000 replicate neighbor-joining trees under the same substitution model.
Before estimating substitution rates, we first tested for recombination using Sawyer's run test (GENECONV program [16]), which detects gene conversion events between pairs of sequences. Two possible gene conversion events were found among three Japanese isolates found in 1997 and 1998: AN40, AN41, and AN66. Specifically, genome positions 1107 to 3863 of AN40 and AN66 were significantly similar (P = 0.00070), as were positions 3805 to 4902 of AN40 and AN41 (P = 0.04972), suggesting AN40 is a recombinant of AN41 and AN66 with a breakpoint between 3805 and 3863. This was confirmed by phylogenetic analyses of positions 654 to 3804 and 3864 to 4902, which showed AN40 changing topological position (data not shown). We also inferred trees for the VP1 and NS1 components of the coding region data set. These topologies were similar with the exception of AN66 (data not shown). Thus, both recombinants, AN40 and AN66, were removed from subsequent analyses.
To determine if the B19 virus phylogenies exhibited adequate temporal structure for an analysis of substitution rates, we plotted the year of isolation of each sequence against its distance from the root of the tree. Figure 1 shows three possible outliers in the large VP1 data set: CHI1, CHI2, and AN85, the latter of which is also an outlier in the coding region data set. These three sequences, which formed a distinct clade with long branches, were therefore removed from the analysis. Final ML trees were constructed (Fig. 2 and 3), and the isolation year was again plotted against root-to-tip distances, with a highly significant correlation observed in both data sets (P = 5.94 × 10−6 and 5.21 × 10−6, respectively).
To estimate viral substitution rates precisely, we used a Bayesian Markov chain Monte Carlo approach (the BEAST package [http://evolve.zoo.ox.ac.uk ]). This method, which assumes a molecular clock, considers differences in branch lengths among viruses sampled at different times and explores different models whose parameters include tree topology, substitution rate, and substitution model. A Bayesian skyline plot, with 10 population size groups, was used to depict demographic structure from which substitution rates could be deduced (5). Input phylogenies were inferred using the HKY85+Γ substitution model with parameters optimized during multiple runs and with Markov chain lengths of 10 million or 100 million in analyses of 24 and 38 sequences, respectively. Mean values are reported as well as 95% high probability density (HPD) intervals. As uncertainty in the data is reflected in the HPD interval, it is imperative that all values within this interval are considered.
Analyses of both the large VP1 and coding region data sets gave high and comparable mean rates of evolutionary change at 1.14 × 10−4 and 1.83 ×10−4 nucleotide substitutions/site/year (Table 1). That the rates and growth curves of each data set were similar also indicated that the demographic model was robust. However, VP1 showed a slightly higher substitution rate than NS1 when analyzed separately. To determine if this difference is significant when the genes are constrained to share the same phylogenetic history, as must be the case for linked loci, we conducted a second analysis on the coding region, this time estimating separate substitution rates for the NS1 and VP1 portions of the alignment (excluding the 8-bp overlap). While VP1 still showed a higher mean rate than NS1 (Table 1), the difference was not significant, as the HPDs overlapped.
That a nonstructural gene evolves at a similar rate as a surface capsid gene suggests that immune selection is not the primary cause of the high substitution rate in B19 virus. This was further supported by an analysis of the ratio of nonsynonymous (dN) to synonymous (dS) nucleotide changes per site, estimated using an ML method (program CODEML) (20, 21). Mean dN/dS ratios (estimated with the M0 model) were 0.074 for the large VP1 data set and 0.135 and 0.079 for the NS1 and VP1 genes of the coding region data set, respectively, showing that purifying selection is dominant. Moreover, we found no evidence of positive selection acting on any individual amino acid residue (P > 0.5). Finally, to confirm that our molecular clock models adequately approximate the evolutionary history of this virus, we also employed the “relaxed” clock model available in the BEAST package and observed similar substitution rates (Table 1) and population growth dynamics.
That our inferred substitution rate, ∼1 × 10−4 substitutions/site/year, is so similar to those estimated for the carnivore parvoviruses suggests high mutation rates may be characteristic of all autonomous parvoviruses, irrespective of their lifestyle. While B19 virus and the carnivore parvoviruses share a genome structure, they exhibit low levels of sequence similarity, infect different hosts, and have unique routes of infection. For the future it will be important to determine whether this elevated rate is restricted to the Parvoviridae or is characteristic of all ssDNA viruses, as suggested by the high levels of diversity found among these viruses, including members of the Geminiviridae and the Circoviridae, most notably human TT virus (6, 9, 11, 12, 15).
It is clear that DNA viruses do not share a common rate of nucleotide substitution, although this does not preclude the existence of a universal rate of mutation per genome per replication (4). As substitution rates are determined by a combination of forces, including the intrinsic frequency of mutation per round of replication, viral generation time, and the extent of natural selection, it will ultimately be necessary to document the contribution of each to the high rates of evolutionary change seen in some DNA viruses. Although cellular proteins are responsible for the replication of all known ssDNA viruses (and indeed all small DNA viruses), it is possible that the required polymerases and/or proofreading proteins do not replicate or repair these unique genomes as accurately or efficiently as they replicate/repair cellular genomes. Consequently, the nature of the viral genome (ssDNA versus dsDNA) and/or the coding capacity—that is, the presence or absence of genes to supplement the host replication and repair machinery and modulate the host immune response—may be partially responsible for different rates of evolutionary change observed in DNA viruses.
Linear regression plots depicting the relationship between the isolation date (year) of each B19 sequence and its distance from the root of the phylogeny (given as number of substitutions per site) for 41 VP1 sequences (A) and 25 coding region sequences (B). Outliers are labeled.
Phylogeny of the 38 B19 VP1 sequences used in the analyses of substitution rates. The tree is rooted with the oldest sampled isolate, and branch lengths are drawn to scale, with nodes showing >70% bootstrap support marked with an asterisk. Names of sequences are given along with the location and date of isolation. Abbreviations: en, England; us, United States; br, Brazil; ir, Ireland; ge, Germany; fi, Finland; uk, United Kingdom; ja, Japan; ko, Korea; sw, Sweden; ch, China.
Phylogeny of 24 B19 virus coding regions used in the substitution rate analyses. The tree is rooted with the oldest isolate. Branch lengths are drawn to scale, and nodes with >70% bootstrap support are marked with an asterisk. Names of sequences are given along with the location and date of isolation. Abbreviations are as described in the legend for Fig. 2.
Nucleotide substitution rates in human B19 erythrovirus
ACKNOWLEDGMENTS
We thank Andrew Rambaut for assistance with BEAST.
This work was supported by a Howard Hughes Medical Institute Fellowship to L.A.S.
FOOTNOTES
- Received 7 November 2005.
- Accepted 13 January 2006.
- Copyright © 2006 American Society for Microbiology
