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

Ethics statement and legal agreement

This research was approved by the Experimental Animal Ethics Committee of Jinzhou Medical University (No. 20180016) and the Agricultural and Rural Comprehensive Service Center of Jinzhou (No. 20180619).

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 and China2-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 μL of diethyl pyrocarbonate-treated 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™ 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 μL 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), China (China-2017, MF036686; China1-2017, MF770715; China2-2018, MH165180), England (VDV-2013, KC786222), and the Netherlands (VDV-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 (Wang et al., 2015; 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.

Nucleotide sequence analysis of DWV isolates

The nucleotide sequences of China1-2017 and China2-2018 were 10,141 and 10,105 nt, respectively. The whole nucleotide sequences of China1-2017 and China2-2018 are enriched in A/U (China1-2017: A-29.17%, U-32.17%, G-22.76%, C-15.90%; China2-2018: A-29.09%, U-32.00%, C-16.16%, G-22.74%). The isolates contained a single major ORF from the 5′–3′ end, which was composed of 8,682 nt (ORF position in China1-2017: 1,146–9,827 bp; China2-2018: 1,135–9,816 bp), encoding 2,894 amino acids. Multiple sequence comparisons showed that the sequences of China1-2017 and China2-2018 were similar to those of previously reported DWVs/VDV-1 strains. Furthermore, compared to other DWV isolates, the nucleotide sequence identity and deduced amino acid sequence identity of China1-2017 and China2-2018 ranged from 84.3% to 96.7% and 94.8% to 98.6%, and 84.3% to 96.7% and 94.7% to 98.4%, respectively (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). The structural and non-structural proteins were positioned in the genomes as follows: VP3 (China1-2017: 1,800–2,537 nt, China2-2018: 1,789–2,526 nt), VP1 (China1-2017: 2,601–3,848 nt, China2-2018: 2,590–3,837 nt), VP2 (China1-2017: 3,849–4,622 nt, China2-2018: 3,838–4,611 nt), helicase (China1-2017: 5,010–6,428 nt, China2-2018: 4,999–6,417 nt), and 3C-RdRp-protease (China1-2017: 7,686–9,827 nt, China2-2018: 7,675–9,816 nt). Moreover, six conserved domains in the helicase, 3C-pro, and RdRp were identified from China1-2017 and China2-2018. Three conserved helicase regions were found in the deduced amino acid sequences of the China1-2017 and China2-2018 ORF, ranging from 1,472 to 1,575, namely, domain A (¹⁴⁷²GxxGxGKS¹⁴⁷⁹), domain B (¹⁵¹⁸Qx₅DD¹⁵²⁵), and domain C (¹⁵⁶¹KKx₄Px₅NTN¹⁵⁷⁵). The 3C-pro conserved domains included the cysteine protease motif (²³⁰⁵GxCG²³⁰⁸) and the putative substrate-binding motif (²³²²GxHxxG²³²⁷). The highly conserved RdRp region, ²⁴⁹⁵TSxGxP²⁵⁰⁰, 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).