Phylogenetic Analysis of Deformed Wing Virus Genotypes from Diverse Geographic Origins Indicates Recent Global Distribution of the Virus

Thirty-seven honeybee specimens, some exhibiting clinical symptoms characteristic of DWV infection, were collected in 2003 and sent to us by collaborating colleagues. The specimens originated from five European countries, Austria (nine samples), Poland (five samples), Germany (three samples), Hungary (five samples), and Slovenia (four samples), as well as from Nepal (two samples), Sri Lanka (two samples), the United Arab Emirates (two samples), Canada (two samples), and New Zealand (three samples). The samples were stored at −20°C until they were investigated.

The bees were homogenized in sterile ceramic mortars with sterile sand; thereafter, diethyl pyrocarbonate-treated water was added. The homogenates were centrifuged at 20,000 × g for 1 min, and 140 μl of the supernatant was used for RNA extraction by employing the QIAamp viral RNA mini kit according to the instructions of the manufacturer (QIAGEN, Hilden, Germany).

Eight primer pairs were selected based on the complete DWV genomes deposited in the GenBank database under accession numbers NC_004830 and AY292384 (Fig. 1B) by employing a primer designer program (Scientific and Educational Software, version 4.1). At first the primers were tested with one Austrian DWV strain by amplifying different genomic regions of the virus. The amplification products were sequenced and compared with the previously mentioned complete genome records from GenBank. For a more comprehensive analysis, four primer pairs that were expected to produce overlapping amplification products corresponding to the partial VP2, putative VP4, and VP1 structural-protein genes and the partial RNA helicase enzyme coding genome region were selected (Fig. 1B). These primer pairs were applied in reverse transcription-PCR (RT-PCR) assays to the remaining 36 DWV samples. Since in the case of the first two primer pairs, the overlapping area was only 19 bases long, another primer pair corresponding to this region was designed in order to create a longer overlapping genome stretch. The primer sequences, orientations, and locations and the expected sizes of the amplified products are shown in Table 1.

DWV RNA was reverse transcribed, and selected regions of the genome were amplified with a continuous-RT-PCR method by using a One Step RT-PCR kit according to the recommendations of the manufacturer (QIAGEN). The amplifications were performed in GeneAmp PCR System 2400 thermocyclers (Applied Biosystems). The temperature profile for the RT-PCR was as follows: 30 min at 50°C (reverse transcription) and 15 min at 95°C (denaturation and HotStarTaq activation), followed by 40 cycles of amplification with 30 s at 94°C, 50 s at 55°C, and 1 min at 72°C. Reactions were completed with a final elongation step for 7 min at 72°C. RNA extracts from bees without DWV infection were used as negative controls. The PCR products were electrophoresed on a 1.2% Tris-acetate-EDTA-agarose gel and stained with ethidium bromide. Bands were photographed with a Kodak DS electrophoresis documentation and analysis system using the Kodak Digital Science 1D software program. Fragment sizes were determined with reference to a 100-bp ladder (Promega).

The PCR products were excised from the gel and extracted by employing the QIAquick gel extraction kit according to the instructions of the manufacturer (QIAGEN). Fluorescence-based direct sequencing of the amplicons in both directions was performed by using the ABI PRISM BigDye Terminator cycle-sequencing ready-reaction kit (21). The nucleotide sequences were identified by the basic local alignment search tool BLAST (2) and aligned by employing the Align Plus program (Scientific and Educational Software, version 4.1). Multiple alignments for phylogenetic analyses were created with the help of the ClustalX program (33). Phylogenetic analyses were conducted by using the PAUP*4.0 beta 10 version (29) with maximum-parsimony, distance (neighbor-joining), and maximum-likelihood criteria. Bootstrap analyses of 1,000 replicates of the parsimony and distance trees were performed. The same data sets were also analyzed by using the PHYLIP package 3.6 beta version (15) with maximum-parsimony, distance (neighbor-joining, fitch, and kitsch), and maximum-likelihood algorithms. Trees were drawn with the help of the TreeView 1.6.6. software.

The DWV sequences described in this paper were submitted to the GenBank database under accession numbers DQ224278 to DQ224311.

At first, eight DWV-specific primer pairs were applied to one DWV-positive Austrian honeybee sample (the infecting strain was designated Austria-1) (Fig. 1B). The PCR products were observed after electrophoresis as clear and distinct bands corresponding to the expected molecular sizes. Amplification products never occurred in the negative controls. The amplification products were sequenced in both directions, and the sequences were aligned with the reference sequences with accession numbers NC_004830 and AY292384 (23). In total, 3,985 nt of the DWV strain Austria-1, or 39% of the complete genome, were determined (Fig. 1C). Sequence alignments revealed a low level of sequence divergence among the DWV genomes investigated and a rather regular distribution of the few nucleotide differences. The Italian strain (NC_004830) exhibited 27 nucleotide substitutions (0.7%) compared to the strain from the United States (AY292384), while strain Austria-1 differed at 28 nt (0.7%) and 40 nt (1%) from the Italian and U.S. strains, respectively. It is also important that the complete genome record of DWV deposited under the accession number NC_004830 (derived from the sequence deposited under accession number AJ489744) contains 69 not-exactly-identified nucleotides (indicated by K, R, W, and Y), of which 23 are located in the genome regions investigated in this study; they were excluded from the previous calculations. Within these loci, the U.S. and Austrian sequences always confirmed the choices, but at eight positions these two strains differed from each other as well; therefore, it is not possible to decide which of the two sequences is more similar to the Italian sequence.

All partial nucleotide sequences of DWV strains found in GenBank were also aligned with the complete genome records and, when applicable, with the sequences obtained in this study (Fig. 1C). A 4,700-nt-long sequence covering the complete nonstructural-protein coding region of a DWV strain from France (France-1) was previously deposited in GenBank under accession number AY224602. The sequence is 98% similar to those of the corresponding regions in the complete genome records of the Italy-1 and USA-1 strains (23) and 98 to 99% similar to the sequences obtained in this study and to other sequences deposited in GenBank. The 5′ untranslated regions of the Italian, U.S., and Austria-1 sequences were aligned with a sequence from a DWV strain detected in Costa Rica (accession number DQ139981). These sequences exhibited 98 to 99% identity to one another. Another genome stretch of the same Costa Rican genotype in the putative VPg protein coding region is available under accession number DQ139982. This region is 98% similar to the corresponding sequences of the Italy-1, USA-1, and France-1 strains. Sequence information on the partial RdRp coding regions of 15 DWV strains was found in GenBank: besides those from the Italy-1, USA-1, and France-1 strains, sequences of the RdRp coding regions from two additional strains from France (accession numbers DQ139981 and DQ139982) (31), one strain from Uruguay (accession number DQ364631) (3), and nine strains from Spain (accession numbers DQ385499 to DQ385507) were available. The RdRp coding regions in these genotypes showed again 98 to 99% similarity to the corresponding regions of the reference sequences.

Because the aforementioned comparisons did not reveal regions of higher diversity than others in the investigated genome stretches of the Italian, U.S., and Austrian DWV strains, two regions were selected for further analyses on the basis of their putative functions. In viruses of vertebrates, frequently the structural proteins (antigens) are subjected to immunoselection, and therefore the corresponding genes often contain hypervariable regions. In contrast, nonstructural proteins (enzymes) are usually more conserved, because mutations may disturb their efficacy. Such tendencies in viruses of invertebrates have not been recorded so far, and therefore, we decided to compare the levels of variability of corresponding partial structural-protein coding regions and partial nonstructural-protein coding regions from DWV samples collected from all over the world: a genome stretch between nucleotide positions 1967 and 2731, located in the VP2, putative VP4, and VP1 structural-protein coding region, and a genome stretch between nucleotide positions 6053 and 6614, located in the putative RNA helicase coding region. RT-PCR assays of the remaining 36 honeybee samples were performed with primer pairs which result in overlapping PCR amplification products. The selected regions were amplified in 33 samples. Interestingly, all three samples originating from New Zealand proved to be negative for these regions in all RT-PCR assays employed. The amplicons from the positive samples were directly sequenced in both directions. The sequences were compiled and aligned by using the complete DWV sequences deposited in GenBank as references. After the completion of these investigations, continuous sequence information was available from every investigated DWV sample for a 764-nt fragment within the putative VP2-VP1 capsid protein gene region and a 562-nt fragment within the helicase gene region. Together, the two fragments cover 13% of the complete DWV genome. The sequences were 98 to 99% similar to those of the reference strains; nucleotide substitutions at 50 positions within the structural-protein coding fragment and at 49 positions within the helicase coding fragment were observed. A more detailed investigation of these positions revealed that at each position, the nucleotides in the studied DWV genotypes always varied between two alternatives. The distribution of the transitions did not show any tendencies: for instance, at position 6430 in 14 genotypes, cytosine was found, and in 23 genotypes, uracil was present; at position 6565, 17 genotypes contained adenine and 20 contained guanine; and at position 6604, 20 sequences exhibited uracil and 17 exhibited cytosine, regardless of the geographic origins of the strains. The putative amino acid sequences corresponding to the two investigated DWV genome stretches were also determined, and these amino acid sequences were aligned and compared. Amino acid substitutions at four positions within the capsid protein region and at seven positions in the helicase protein region were found (Table 2).

Phylogenetic trees were constructed based on the aforementioned nucleotide sequences to illustrate the probable genetic relationships among the genotypes. The corresponding region of the DWV strain France-1 (accession number AY224602) (17) was also included in the alignment and phylogenetic tree of the helicase coding regions. Two closely related viruses, KGV (accession number NC_005876) (16) and Varroa destructor virus 1 (VDV-1; NC_006494) (27), were included in the analyses, and VDV-1 was used as the outgroup. The phylogenetic trees based on the structural-protein coding regions and the helicase coding regions are shown in Fig. 2 and 3, respectively. As was already indicated by the high levels of genetic similarity, the DWV genotypes exhibit short genetic distances between each other. Certain sequences proved to be essentially identical to one another, and in general very close genetic relationships are indicated by the trees, independent of the geographic origins of the samples. Different statistical algorithms (maximum-parsimony, distance matrix, and maximum-likelihood criteria) were applied; nonetheless, essentially the same tree structures were obtained by using PAUP* (29) and PHYLIP (15) software packages. Therefore, only the phylogenies resulting from the PAUP analysis are shown in Fig. 2 and 3. Due to the low level of sequence divergence, in many cases the bootstrap analysis revealed low values of confidence, indicating that the present clustering is statistically poorly supported. Since the bootstrap values are of no significance in this case, they have been omitted in Fig. 2 and 3. KGV and VDV-1 form distinct monophyletic branches, and their separation from the DWV genotypes is supported by bootstrap analyses. The DWV sequences derived from bee samples collected in remote geographic regions do not show phylogenetic segregation: in some cases, viruses from Europe cluster together with viruses from Asia. The North American sequences are relatively separated; however, in the tree based on sequences within the helicase coding region, they show a close relationship to the sequence of a DWV strain collected in the United Arab Emirates. In comparing the two trees, very few similarities are found in the structures of clustering and in the members of the various subbranches. To summarize, the DWV phylogeny based on the genetic analysis and comparison of viruses collected in various parts of the world many thousands of kilometers apart indicates a high degree of genome conservation in DWV.