Prevalence and Molecular Characterization of Deformed Wing, Acute Bee Paralysis and Black Queen Cell Viruses Infecting Honey Bees and Varroa Mites

Surveys were conducted during 2018 and 2019 to forty different apiaries in the Malatya (25 apiaries) and Elazığ (15 apiaries) provinces and honey bees and parasite Varroa mite feeding on the honey bees were collected. In total, 147 honey bee and thirty Varroa samples were obtained from Malatya (Battalgazi, Doğanşehir and Yeşilyurt districts) and Elazığ (Arıcak, Kovancılar, Sivrice and Palu districts). Eighty-four honey bee and fifteen Varroa mite samples were obtained from twenty-five different apiaries in Malatya province, while sixty-three honey bee and fifteen Varroa mite samples were taken from fifteen different apiaries in the Elazığ province. There were at least ten hives in each surveyed apiary. Honey bee and parasite Varroa mite samples were collected by sampling a representative number of bees from several hives in each apiary. In apiaries with ten to forty at least two samples were collected for every ten hives depending on size.

Honey bee samples showing paralysis, sudden death, change in colour (becoming dark brown to black), trembling, enlarged abdomens with wings dislocating, becoming hairless, dark to shiny black in colour were collected as well as the symptomless samples. Collected live honey bee and Varroa mite samples were taken into sterile tubes covered with cheesecloth in cold chain and were taken to the laboratory. They were numbered, photographed and classified according to the symptom type. Vector Varroa mites were collected via pliers under the microscope and taken into sterile micro tubes for RNA extraction.

Total RNA was extracted from live honey bee and Varroa mite samples using a commercial genomic RNA purification kit (GeneJET Plant RNA Purification Kit). Samples were prepared and extracted individually. Virus-specific primers were used to perform reverse transcription. DWV-specific primers used in the study (F-5′-TTG GTA TGC TCC GTT GAC TG-3′, R-5′-ATT CCT CAG AAG TTG GTT TCG-3′) were encoding an RdRP gene region with an expected band size of 488. BQCV-specific primers (F-5′-GAC AGC GTG CCA AAG AGA G-3′, R-5′-GCG AAC CCG TCC AAT ACT TA-3′) and ABPV-specific primers (F-5-GTA TGG AAG TGG GCT GAG GA-3′, R-5-CGC GGT ACT AAA AAG CTA CGA-3′) were encoding a part of the coat protein gene with an expected band size of 567 bp and 476 bp, respectively (Rüstemoğlu & Sipahioğlu, 2019). To run RT-PCR, a total of 25 μl of reverse transcription mixture was used, which included 5 μl of 5x RT buffer, 1 μl of MMLV reverse transcriptase enzyme, 1 μl 100 nmol of each primer, 5 μl of dNTP (each 2.5 mM), 0.5 μl of RNase inhibitor, 7 μl RNase free water, and 2 μl total RNA.

PCR tests were performed according to Rüstemoğlu & Sipahioğlu (2019) with slight modifications. A 25 μL of PCR reaction volume containing 14.5 µl RNase free water, 2.5 µl 10x GoTaq Green Buffer (Promega, Madison, WI, USA), 0.5 µl Taq DNA Polymerase, 1 µl dNTP (20 mM each), 1.5 µl 25 mM MgCl, 1 μl of each forward and reverse primer (100 mM each), and 3 µl cDNA for DWV, BQCV and ABPV were prepared. Thermo Scientific Arktik Thermal Cycler (Waltham, MA, USA) was used to perform PCR with the following cycling parameters: initial denaturation at 94°C for 2 min followed by thirty-five cycles of denaturation at 94°C for 1 min, annealing at 55°C for 30 seconds for DWV and BQCV and 57°C for 30 seconds for ABPV, extension at 72°C for 1 min, and final extension at 72°C for 10 min. The electrophoresis of fifteen μl of each amplicon was performed on 2% agarose gel in 1x TAE buffer (40 mM Tris-acetate, 1 mM EDTA, pH 8.0) and stained with Pronasafe nucleic acid staining solution (CondaLab). As a marker, a DNA ladder was employed (100bp DNA ladder, Thermo Scientific).

Amplified PCR products were sequenced bidirectionally by a commercial firm (BM Labosis, Ankara/Türkiye) using the dideoxy chain termination reaction. The nucleotide sequences were compared to other isolates from the NCBI database online, revealing their phylogenetic relationships with viruses from throughout the world.

Randomly chosen viral isolates that were positive in RT-PCR tests (eight DWV, four ABPV and two BQCV isolates) were sequenced and compared to other viral isolates in the NCBI GeneBank, with differences and similarities investigated. Sequences of DWV, BQCV and ABPV generated in this study, as well as sequences of closely related species downloaded after a BLAST search from the GenBank database (http://www.ncbi.nlm.nih.gov/), were included in separate datasets in the Geneious Prime (GP) software (Version 2023.0.1). On the Geneious alignment, each dataset was individually aligned, and discrepancies were manually rectified. For phylogenetic inference, Maximum Likelihood analysis (ML) was used with 1000 bootstrap repetitions in the same software for each dataset. Phylogenetic analyses were done with Geneious tree builder.

According to the RT-PCR tests performed with specific primers, thirty-five honey bee samples from the survey region were found to be infected with DWV, nineteen honey bee samples were found to be infected with ABPV, and fifteen honey bee samples were found to be infected with BQCV. Viruses were found in almost all of the symptomatic samples as well as the symptomless ones. Two of the thirty-five DWV positive samples, four of the nineteen ABPV positive samples, and three of the fifteen BQCV positive samples were asymptomatic. Almost all of the samples with enlarged abdomens with wings dislocating, trembling and paralysis were found to be infected with DWV. Honey bees found to be infected with BQCV showed change in colour (becoming dark brown to black) and ABPV positive honey bees had hairless abdomens and showing trembling, paralysis and death symptoms.

The infection rates were found to be 23.81%, 12.93% and 10.20% for DWV, ABPV and BQCV, respectively (Tab. 1). Based on the results of the viral-specific PCRs, out of 147 samples, four (2.72%) were found to be coinfected with DWV and BQCV, whereas one (0.68%) was found to be coinfected with DWV+BQCV+ABPV in Malatya province. No mixed infections were detected in Elazığ province. Totally 30 Varroa mite samples were also collected and tested against the mentioned viruses in both Malatya and Elazığ provinces in this study. Five out of thirty Varroa samples (two collected from Malatya and three collected from the Elazığ province) were found to be infected with DWV. The total infection rate of Varroa mites was found to be 16.66% for DWV. None of the samples tested were found to be positive for ABPV and BQCV (Tab. 1).

Randomly selected DWV positive isolates were sequenced and compared with other viral isolates in the NCBI GeneBank and differences and similarities with them were investigated. The resulting DWV sequences (OP805878, OP805879, OP805880, OP805887, OP805888, OP805889, OP805890, OP805891) were named as Malatya 1, Malatya 2, Malatya 3, Elazig 1, Elazig 2 (mite), Elazig 3, Elazig 4 and Elazig 5. Based on nucleic acid sequence similarity, the phylogenetic tree was divided into five major clusters. The isolates in this study were in the first group, together with DWV isolates from A. mellifera honey bees from the United Kingdom, Lebanon, Türkiye (Hakkari), France, Germany, and Saudi Arabia (Fig. 2).

Fig. 2.

Multiple alignments showed that DWV isolates revealed from this study has 96.34–98.06% nucleotide identity with the other global isolates retrieved from GenBank (Tab. 2). It was interesting that Malatya and Elazığ isolates were found to be similar with a Japan isolate accession no. B070959, which is a Kakugo virus (KV) (Iflavirus) with the similarity rate of 96.61% (Tab. 2).

When DWV isolates revealed in this study were compared with ten different isolates with the world and one isolate from Türkiye, it was observed that Elazığ 1 isolate had an insertion mutation in 107^th nucleotide, Malatya 1 and Malatya 2 isolates had five deletions in 62^nd–63^rd nucleotides and six frameshift mutations, 32 deletions and 29 frameshift mutations in Malatya 3, Elazığ 2, Elazığ 3 and Elazığ 4 isolates in 128^th and 129^th nucleotides. There was no mutation in Elazığ 5 isolate. There was one insert mutation in Elazığ 1 isolate. In total thirty-seven deletions in Malatya 1 (two mutations), Malatya 2 (three mutations), Malatya 3 (eight mutations), Elazığ 2 (eight mutations), Elazığ 3 (eight mutations) and Elazığ 4 (eight mutations) mutations were observed. There was also thity-five frameshift mutations in Malatya 1, (four mutations), Malatya 2 (two mutations), Malatya 3 (seven mutations), Elazığ 2 (seven mutations), Elazığ 3 (eight mutations) and Elazığ 4 (seven mutations) isolates (data not shown).

ABPV sequences (OP805881, OP805882, OP805883 and OP805884) were named as Malatya 1, Malatya 2, Malatya 3 and Malatya 4 and their multiple alignments showed that these isolates revealed 96.76% and 96.70% nucleotide identity with Ankara and Tekirdağ (Türkiye) isolates and 89.62–96.46% with the other world isolates retrieved from GenBank (Tab. 3). Based on nucleic acid sequence similarity, the phylogenetic tree was divided into two major clusters. Malatya 1 and Malatya 2 isolates were in the first group forming another branch different than the other isolates, whereas Malatya 3 and Malatya 4 were in the second group distinct from all the other isolates (Fig. 3). The similarities and differences were determined by comparing the base sequence of coat protein genes with other ABPV isolates deposited in GeneBank and a shift mutation in the 96^th nucleotide of ABPV Malatya 1 isolate, and an insertion mutation and a shift mutation in 98^th nucleotide of ABPV Malatya 2 isolate were observed, whereas forty-four deletion mutations and twenty shift mutations were observed in the 164^th and 167^th nucleotides of ABPV Malatya 3 and ABPV Malatya 4 isolates (data not shown).

Fig. 3.

BQCV sequences (OP805885 and OP805886) were named as Malatya 1 and Malatya 2 in the GenBank/EMBL database and their multiple alignments showed that BQCV isolates revealed 96.78% nucleotide identity with another Turkish isolate from Hakkari province and 87.40–90.21% with the other world isolates retrieved from GenBank (Tab. 4). Based on nucleic acid sequence similarity, the phylogenetic tree was divided into six major clusters. The isolates revealed in this study were in the second group together with the Türkiye (Hakkari) isolate of BQCV from A. mellifera honey bee (Fig. 4).

Fig. 4.

Also, the BQCV isolates obtained in this study were compared via CP gene sequence with the other isolates deposited in the GeneBank and the similarities and differentiations were determined. According to the results, seven shift mutations were observed in the 525^th nucleotide of Malatya 1 isolate, and fifteen shift mutations and one deletion mutation were observed in the 527^th nucleotide in BQCV Malatya 2 and BQCV Elazığ 1 isolates (data not shown).