Phylogenetic and recombination analyses of two deformed wing virus strains from different honeybee species in China

Background Deformed wing virus (DWV) is one of many viruses that infect honeybees and has been extensively studied because of its close association with honeybee colony collapse that is induced by Varroa destructor. However, virus genotypes, sequence characteristics, and genetic variations of DWV remain unknown in China. Methods Two DWV strains were isolated from Jinzhou and Qinhuangdao cities in China, and were named China1-2017 (accession number: MF770715) and China2-2018 (accession number: MH165180), respectively, and their complete genome sequences were analyzed. To investigate the phylogenetic relationships of the DWV isolates, a phylogenetic tree of the complete open reading frame (ORF), structural protein VP1, and non-structural protein 3C+RdRp of the DWV sequences was constructed using the MEGA 5.0 software program. Then, the similarity and recombinant events of the DWV isolated strains were analyzed using recombination detection program (RDP4) software and genetic algorithm for recombination detection (GARD). Results The complete genomic analysis showed that the genomes of the China1-2017 and China2-2018 DWV strains consisted of 10,141 base pairs (bp) and 10,105 bp, respectively, and contained a single, large ORF (China1-2017: 1,146–9,827 bp; China2-2018: 1,351–9,816 bp) that encoded 2,894 amino acids. The sequences were compared with 20 previously reported DWV sequences from different countries and with sequences of two closely related viruses, Kakugo virus (KV) and V. destructor virus-1. Multiple sequence comparisons revealed a nucleotide identity of 84.3–96.7%, and identity of 94.7–98.6% in amino acids between the two isolate strains and 20 reference strains. The two novel isolates showed 96.7% nucleotide identity and 98.1% amino acid identity. The phylogenetic analyses showed that the two isolates belonged to DWV Type A and were closely related to the KV-2001 strain from Japan. Based on the RDP4 and GARD analyses, the recombination of the China2-2018 strain was located at the 4,266–7,507 nt region, with Korea I-2012 as an infer unknown parent and China-2017 as a minor parent, which spanned the entire helicase ORF. To the best of our knowledge, this is the first study to the complete sequence of DWV isolated from Apis cerana and the possible DWV recombination events in China. Our findings are important for further research of the phylogenetic relationship of DWVs in China with DWV strains from other countries and also contribute to the understanding of virological properties of these complex DWV recombinants.

virus (KV) (Mcmahon et al., 2016). V. destructor virus-1 (VDV-1) was classified as DWV-B, with an overall RNA genome of 84% to those of the classic DWV-A. In addition, the first 1,455 nt of the ORF encoding the lower molecular mass structural proteins shows the greatest diversion from those of the classic DWV-A, with an RNA identity of 79%, and translates to a polypeptide of 485 aa with an identity of 90% (Ongus et al., 2004;Ryabov et al., 2014). Recent studies have shown that recombination events have frequently occurred between these two main variants (Mordecai et al., 2015;Zioni, Soroker & Chejanovsky, 2011). Thus far, the complete genome sequences of approximately 20 strains have been sequenced, with isolates mostly obtained from A. mellifera and Vespa crabro; however, there are no complete genome sequences of DWV isolates from A. cerana (Reddy et al., 2013;Lamp et al., 2016;Forzan et al., 2017). The epidemiology of DWV in A. cerana and A. mellifera has been investigated in China; however, virus genotypes, sequence characteristics and genetic variations of DWV remain unknown (Zheng et al., 2015;Chao et al., 2017).
In the present study, we isolated two DWV strains (China1-2017 and China2-2018) from A. mellifera and A. cerana, respectively, and produced the complete nucleotide sequences of DWV from A. cerana. The sequences of the China1-2017 and China2-2018 isolates were analyzed and compared to the reference nucleotide sequences of DWV genotypes from other countries. Furthermore, we analyzed molecular biological characteristics and phylogenetic relationship of DWV structural and non-structural polyprotein regions of the China1-2017 and China2-2018 strains and reference strains. Moreover, we performed a recombination analysis of DWV from Chinese isolates using recombination detection program (RDP4) software and a genetic algorithm for recombination detection (GARD).

Sample collection
Samples were collected from two different regions in China. A total of 257 A. mellifera samples originating form Jinzhou of Liaoning Province (41 14′34″N, 121 8′15″E) were collected in March 2016, and 463 A. cerana worker bee samples originating from Qinhuangdao of Hebei Province (40 3′26″N, 119 33′4″E) were collected in April 2017. All sampled individuals exhibited typical deformed wing symptoms. Upon collection, the bees were stored on ice and immediately transported to the laboratory, where they were kept frozen at -80 C until analysis.

Virus isolation
Virus isolation was performed according to previously published methods (Ying et al., 2016;Mingxiao et al., 2011;Jakubowska et al., 2016) Briefly, 50-60 adult bees were completely homogenized with a mortar and pestle in 10 mL of phosphate-buffered saline (pH 7.4), containing 0.5% nonylphenol ethoxylate, then the homogenate was incubated for 30 min at 20 C. After this, the homogenate was centrifuged at 5,000Âg for 20 min at 4 C. Large debris was removed, and the supernatant was again centrifuged at 8,000Âg for 30 min at 4 C. After discarding the precipitate, the supernatant was placed in an ultracentrifuge tube and centrifuged at 82,000Âg for 1 h at 4 C. The resulting pellet was resuspended in two mL STE buffer (10 mM Tris-HCl, 0.1M NaCl, and one mM EDTA; pH 7.3). The viral strains of A. mellifera and A. cerana were named China1-2017 andChina2-2018, respectively. Viral RNA isolation and DWV screening by PCR Total RNA was extracted from the purified virus of China1-2017 and China2-2018 strains using a TIANamp Virus DNA/RNA Kit (Tiangen, Beijing, China), according to the manufacturer's instructions. Total RNA was eluted in 30 mL of diethyl pyrocarbonatetreated water and stored at -80 C until further analyses. To determine the presence of DWV in A. mellifera and A. cerana, a RT-PCR assay was performed using a PrimeScript TM RT-PCR Kit (TaKaRa, Dalian, China), according to the manufacturer's instructions.
The primers DWV-F (5′-TTTGCAAGATGCTGTATGTGG-3′) and DWV-R (5′-GTCGTGCAGCTCGATAGGAT-3′) were used to amplify a 395 base pairs (bp) fragment of the DWV RdRp gene (accession number: AY292384). The PCR amplification was carried out under the following conditions: 94 C for 2 min, followed by 25 cycles of 94 C for 30 s, 55 C for 30 s, 72 C for 30 s, and a final extension of 72 C for 5 min. A sample of six mL of PCR product was loaded on a 1.2% agarose gel containing Gelstain and analyzed using a Tanon 2500 Digital Gel Image Analysis System. The PCR products were sequenced commercially (Sangon Biotech, Shanghai, China).

RT-PCR amplification and genome sequencing
Viral genomic RNA of China1-2017 and China2-2018 were reverse-transcribed to cDNA with an oligo (dT) primer (Sambrook & Russell, 2001) (TransGen Biotech, Beijing, China), according to the manufacturer's recommendations. A total of 17 primer pairs were designed to amplify the complete genome sequence of China1-2017 and China2-2018, using the complete DWV sequence of the USA and KOR strains (accession number: AY292384 and JX878304; Table 1) (Lanzi et al., 2006;Reddy et al., 2013). PCR amplifications were performed in Eppendorf tubes, and the cycling protocol for RT-PCR amplification was as follows: 45 min at 42 C (reverse transcription), followed by 30 cycles at 95 C for 60 s, 50-55 C for 30 s, and 72 C for 60 s. The 3′ and 5′ termini of the China1-2017 and China2-2018 strains were obtained by employing a rapid amplification of cDNA ends (RACE) technique, using a SMARTer RACE 5′/3′ kit (Clontech, Beijing, China). PCR products were electrophoresed and purified, and subsequently sequenced commercially (Sangon Biotech, Shanghai, China). The nucleotide sequences of all fragments were assembled to compile the genomic sequences of the China1-2017 and China2-2018 strains using Lasergene software (DNASTAR), based on published complete DWV sequences (accession number: AY292384, JX878304, and JX878305).

Sequence and phylogenetic analyses
To research the DWV genome sequence characteristics of the two novel strains, 20 complete genome sequences of DWV were obtained online from the National Center for Biotechnology Information (http://www.ncbi.nlm.nih.gov/BLAST/) ( Table 2). Multiple alignments of nucleotides and amino acid analyses were performed using the ClustalW program in the MegAlign software program (DNASTAR, Madison, WI, USA). At the same time, the conserved domains of the nonstructural proteins in the C-terminal portion of the polyprotein were analyzed with reference to the functional domains of mammalian picornaviruses. To investigate the genetic variation and evolutionary characteristics of the isolates, DWV complete genome sequences of the classic Type A and Type B strains from 2001 to 2018 and two closely related viruses, KV (Fujiyuki et al., 2004) and VDV-1 (Ongus et al., 2004) in the GenBank database were selected for further analysis, including four isolates in Asia, 13 in Europe, three in America, and one in Oceania. Then, phylogenetic trees of ORF, VP1 and 3C+RdRp were reconstructed by the maximum likelihood method using MEGA5. The best-fitting nucleotide substitution model (GTR+I+G) was used for the all alignment datasets, which was determined by the lowest BIC score in MEGA 5.0 (Tamura et al., 2011). A total of 1,000 replicates of bootstrap resampling were used to ensure the reliability of individual nodes in each phylogenetic tree.

Recombination analysis
Recombination analysis was performed using eight full-length genomes of DWV from Japan (KV-2001, AB070959), Korea (Korea I-2012, JX878304; Korea II-2012, JX878305),  -2004, AY251269), and the detection of inter-strain recombination, identification of closest parental sequences and localization of possible recombination break points were assessed using RDP4 software, which comprises RDP, GeneConv, Bootscan, MaxChi, Chimaera, SiScan, and 3SEQ algorithms. The standard settings of these algorithms were used with the default values of RDP4. The likelihood of recombination events was significant in at least four algorithms at P < 1.0E-6 or recombination consensus scores (RCS) above 0.6, based on the RDP4 analysis. Recombination events were considered possible when the P-value of at least three algorithms was below 0.05, and the RCS was between 0.4 and 0.6, and the likelihood of recombination events was considered insignificant when the RCS was under 0.4 with P < 0.05 Lee et al., 2017;Gao et al., 2018). Moreover, the recombination events were further verified by GARD implemented in the Datamonkey web interface (Delport et al., 2010) and the credibility of the recombination breakpoints was assessed by the KH test.

Nucleotide sequence accession number
The DWV nucleotide sequences of the strains China1-2017 and China2-2018 are accessible on GenBank (accession number: MF770715 and MH165180).

RT-PCR detection of DWV samples
After virus isolation, specific primer pairs were used to detect DWV in samples from Jinzhou and Qinhuangdao by RT-PCR, and fragments of approximately 395 bp were produced (Fig. 1). After sequencing and alignment, the nucleotide sequence identity exceeded 97% by BLAST; thus, we successfully isolated two novel DWV strains from A. mellifera and A. cerana in China.  (Table S1).

Amino acid sequence analysis of DWV isolations
As previously reported for mammalian picornaviruses, the ORF of China1-2017 and China2-2018 strains encoded a 2,894-amino-acid polyprotein, with structural proteins at the N-terminus and non-structural proteins at the C-terminus (Lanzi et al., 2006;Organtini et al., 2016). TSxGxP 2500 , was recognized between the deduced amino acid positions 2,495 and 2,500.

Phylogenetic relationships of the DWV isolates
To assess the genetic relationships of the DWVs, three phylogenetic trees were constructed based on the VP1, 3C+RdRp segments, and ORF gene sequences. The results showed that the ORF and VP1 groups better explained the geographical distribution of DWV, and the 3C+RdRp-coding region better explained the genotype and diversity of DWVs.
In the ORF gene phylogenetic tree, the 24 DWV isolates (including China1-2017 and China2-2018) were divided into two groups (lineage A and B; Fig. 2A). The first group contained two master lineages (lineage A1 and A2), one of which included eight isolates from America and Europe; the second lineage included six isolates from Asia. China1-2017 and China2-2018 belonged to lineage A2. Lineage B contained six isolates from Europe. We found that the phylogenetic tree of the ORF among DWV isolates correlated with the geographical distribution.
The phylogenetic tree based on the VP1 segment produced two groups (lineages C and D; Fig. 2B) and was similar to the ORF tree ( Fig. 2A). The first group (lineage C) was further divided into two sub-groups (lineages C1 and C2): ten isolates from America, Europe, and Oceania formed lineage C1, whereas eight isolates formed lineage C2. The second group (lineage D) contained six isolates from China, Korea, Japan, and Europe.
The phylogenetic tree based on the 3C+RdRp segments produced two distinct groups (lineages E and F; Fig. 2C). The 22 isolates in lineage E belonged to DWV Type A, which included variants from America, Europe, Oceania, and Asia. Lineage F belonged to the classic DWV Type B, which only contained two isolates from Belgium and the Netherlands.

Recombination analysis of the DWV isolates
To explore potential recombination signals in the DWV isolates from Asia, recombination signals were assessed using the RDP4 software. Using seven algorithms, nine recombination events were detected in Asian strains ( Fig. 3B; Table 3). In all potential recombination events, three recombination events (events 2, 3, and 5) had a high degree of certainty based on the RDP4 software standard (Table 3). However, GARD analyses indicated that only one isolate (event 3), China2-2018, was identified as a recombinant at the breakpoint in the positions 4,266 and 7,507 nt with a high level of confidence (LHS, RHS P-values < 0.01). Based on the above analysis, event 3 was identified as the real recombination event. In event 3, the recombination of strain China2-2018 was located at the 4,266-7,507 nt region, with Korea I-2012 as an infer unknown parent and China-2017 as a minor parent (Figs. 4A-4C), which spanned the entire helicase ORF (Fig. 3).

DISCUSSION
Deformed wing virus is one of the most prevalent, pathogenic honeybee viruses in the world, and has been directly linked to colony collapse disorder (Organtini et al., 2016;Kevill et al., 2017). Despite the importance of DWV as a honeybee virus, only a limited number of complete genome sequences for DWV are available. In the present study, we determined the entire genome sequences for two DWV isolates collected from A. mellifera and A. cerana in China. Generally, DWV is only prevalent in A. mellifera and V. destructor and is not common in A. cerana (Tentcheva et al., 2004;Xie, Huang & Zeng, 2016;Zhang & Han, 2018). However, we obtained the complete genome sequence from the China2-2018 strain A. cerana for the first time in the present study, which helps to investigate the host range of DWV. The comparison of China1-2017 and China2-2018 showed 96.7% identity of nucleotide sequences and 98.1% identity of amino acids. Although the identity of the two isolates was relatively high, they did not belong to the same strain. The isolated strain sequences of similarities with 20 reference strains ranged from 84.3% to 96.4%, and the highest sequence identity was assigned to strain UK-2009 (GU109335), whereas Belgium-2016 (KX783225) and VDV-2004 had the lowest sequence identities. Therefore, the two novel isolates were significantly different from DWV Type B. The amino acid identity ranged from 94.8% to 98.6%, and the highest similarity belonged to the  Table S1, virus strains from the same continents or from the same countries showed higher levels of identity of the nucleotides and amino acids, including the novel isolates; therefore, these viruses have been present in the honeybee populations for a long time and the viruses have evolved more or less independently. Moreover, six highly conserved motifs were identified from the two novel isolates, featuring the typical characteristics of iflaviruses. Based on the phylogenetic analyses of ORF, VP1, and 3C+RdRp segments, the two novel isolates from China clustered within the same clade as other DWV strains form Asia; therefore, the novel isolates had a closer relationship with other Asian strains, and the DWV strains from Asia might have originated from the ancestor KV-2001 (AB070959). Furthermore, the phylogenetic analysis of DWV based on ORF and VP1 revealed numerous different geographically determined clades and the phylogenetic tree of VP1 also presented a clearer pattern regarding the geographical distribution; therefore, the strains with the closer relationship had different evolutionary rates in different environments and hosts. Moreover, the phylogenetic tree based on the VP1 sequences confirmed that China1-2017 and China2-2018 were more closely related to the isolates from Korea than to those from Japan. In RNA viruses, 3C-pro and RdRp genes tend to be highly conserved; therefore, 3C-pro and RdRp genes have usually been used to distinguish subtype classification of RNA viruses (Baker & Schroeder, 2008;Geng et al., 2014;Kevill et al., 2017). In the present study, the phylogenetic analysis of 3C+RdRp showed that the two novel isolates all belonged to DWV Type A. With a difference in the ORF and VP1 tree, the China2-2018 (MH165180) strain was more closely related to the Korea-2012 strain (JX878304) than to the China1-2017 strain (MH770715). We suspect that the VP1 and 3C+RdRp genes had generated at a different evolutionary rate in different environments, which led to the diversity of the phylogenetic trees, because the structural proteins genes (VP1, VP2, and VP3) are more likely to mutate in the picornavirus (Amin et al., 2014;Hu et al., 2016).
Natural recombination is an important strategy for viruses to adapt to new environmental conditions and hosts. Recombination events have been observed in DWVs (Seo et al., 2009;Lian et al., 2013;Dalmon et al., 2017); however, the recombination events of Chinese DWVs have never been described. Based on geographical location, we used all DWV strains from Asia, the classic DWV Type B and two novel isolates from China in the present study to analyze the potential recombination events using RDP4 software and GARD, and the recombination event of the China2-2018 strain was confirmed (Table 3; Fig. 4). In this recombination event, the strains from China and Korea were mainly involved, including the novel isolates; therefore, the DWV recombination widely exists in honeybee colonies of East Asia. In general, irregular and complicated Figure 4 Recombination analysis of full-length DWV genome sequences of China-2017 and China2-2018 isolations by RDP4. An alignment of the eight DWV genomes from Asia and Europe was analyzed using the RDP4 software. Two potential recombination events were identified for Chinese strains (illustrated by A-F, respectively). The left part illustrates the results of RDP analyses (A and D). The right part presents phylogenetic analyses based on the full-length genome (excluding one of the terminal direct repeats) excluding the region of recombination (B and E) or based on the recombination region only (C and F) using UPGMA in MEGA5.0 with 1,000 replicates. Values on internal branches refer to the percentage of bootstrap replicates in which the branch was found; only values greater than 50% are shown. The scales illustrate the number of substitutions per nucleotide. The color code used is described at the top.
Full-size  DOI: 10.7717/peerj.7214/ fig-4 recombination patterns indicate that the recombination events are usually random, although a detailed understanding of the mechanism involved in such recombination phenomenon must be clarified. When compared with the DWV recombinant isolated by Dalmon et al., we found that all recombinant sites located in the region encoded non-structural proteins. The apiculture characteristic of China is probably the greatest factor that led to the DWV epidemic in A. cerana colonies. In China, there are a large number of A. mellifera and A. cerana colonies found in the same region during the nectar collecting season, which could promote the transmission of DWV from A. mellifera to A. cerana by pollen. In this interaction among the virus, host, and infectious vector, the recombination strains continuously appear. Future studies are required to discover more recombination phenomenon and the occurrence of recombination events may contribute to the high levels of genetic diversity and viral adaptability to host in DWVs, which may increase the potential of this virus to threaten successful beekeeping.

CONCLUSIONS
In summary, we found two novel DWV isolates from China and reported the first complete genome sequence of DWV from A. cerana. Based on phylogenetic trees, the novel DWV isolates from China were confirmed to be the closest related to the strain from Korea. Furthermore, the recombinant phenomenon was discovered in the novel isolates of China2-2018, which is the first description of DWV recombination in China. Our study has not only revealed the presence of novel DWV recombinants in China but also provides information that may be useful for further research on the phylogenetic origins of Chinese DWV strains.

ADDITIONAL INFORMATION AND DECLARATIONS
Funding