# Genotype B of deformed wing virus and related recombinant viruses become dominant in European honey bee colonies

Fabrice Sircoulomb, Eric Dubois, Frank Schurr, Pierrick Lucas, Marina Meixner, Alicia C. Bertolotti, Yannick Blanchard, Richard Thiéry

► To cite this version:

Fabrice Sircoulomb, Eric Dubois, Frank Schurr, Pierrick Lucas, Marina Meixner, et al.. Genotype B of deformed wing virus and related recombinant viruses become dominant in European honey bee colonies. Scientific Reports, 2025, 15, 10.1038/s41598-025-86937-5. anses-04969098

HAL Id: anses-04969098

https://anses.hal.science/anses-04969098v1

Submitted on 27 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Check for updates

OPEN

# Genotype B of deformed wing virus and related recombinant viruses become dominant in European honey bee colonies

Fabrice Sircoulomb<fnref n="1" />, Eric Dubois<fnref n="1" />✉, Frank Schurr<fnref n="1" />, Pierrick Lucas², Marina Meixner³, Alicia Bertolotti<fnref n="1" />, Yannick Blanchard²,⁴ & Richard Thiéry<fnref n="1" />,⁴

The *Varroa destructor* mite's transmission of deformed wing virus (DWV) to honey bees is responsible for most winter mortalities of colonies worldwide. Four DWV genotypes (A, B, C and D) and numerous recombinants have been described. The most recent studies have reported the greater prevalence of DWV-B over DWV-A in several countries, including European ones, while C and D genotypes appear rare or extinct. However, no global evaluation of DWV-A and DWV-B distribution was available at the European level to date. In this study, we quantified both DWV genotypes by real-time PCR from pools or individual honey bees and from *V. destructor* mites sampled in 15 European countries between 2010 and 2017. These data and the sequencing of the viral RNA provide a first insight into DWV diversity, with a clear dominance of DWV-B and recombinants (A/B) in Europe. Chimeric sequencing reads were used to locate the recombinant junctions along the DWV genome. These were not randomly distributed, but mainly clustered in three genomic areas: the 5'UTR, leader peptide and helicase coding sequences. In our study, the DWV recombinant genomes shared at least the VP1-VP3 coding sequences with the DWV-B. Further studies are needed to explore the apicultural context explaining these differences in DWV genotype dominance.

The health of honey bee colonies is threatened by numerous parasites and pathogens, including viruses<fnref n="1" />. Deformed wing virus (DWV, species: *Iflavirus aladeformis*), originally regarded as a low-impact stressor to honey bees, is now considered as one of the main viral threats responsible for winter losses in honey bee colonies around the world²,³. The clinical impact of DWV on honey bee health began following the introduction of the *Varroa destructor* mite in the second part of the 20th century into Europe and America³. It is unclear whether DWV infected the western honey bee or not before the *V. destructor* mite switched hosts from *Apis cerana* (eastern cavity-nesting honey bee) to *Apis mellifera* (western honey bee)⁴. It is known that DWV can be transmitted vertically or horizontally by oral acquisition and persist as a covert infection even in the mite's absence⁵⁻⁹. By directly inoculating virus particles into the haemolymph while feeding on larvae or adult honey bees, the *V. destructor* mite acts as a biological vector and significantly enhances the pathogenicity of DWV¹⁰⁻¹². The horizontal transmission of DWV by the *V. destructor* mite has been shown to be associated with morphological symptoms such as deformed wings, a shortened abdomen, behavioural disorders such as foraging at an earlier age¹³, reduced distance and duration of flights¹⁴, a decrease in learning performance¹⁵,¹⁶ and a reduced lifespan¹⁷. It has also been suggested that the synergistic suppression of honey bee immunity by both the mite and the virus could be an additional factor leading to this change in DWV pathogenicity¹⁸,¹⁹. Moreover, the mite's transmission of DWV has been associated with a reduction in DWV genetic diversity¹⁰,¹² and a higher viral load observed in colonies and single honey bee heads²⁰, suggesting the selection of better-adapted viral variants to this new context of host-pathogen-vector interactions²¹⁻²⁴. The *V. destructor* mite has notably been associated with the spread of an emerging DWV-related virus, Varroa destructor virus-1 (VDV1)⁹. This emerging virus was genotyped as DWV-B, distinguishing it from the original DWV variants (genotyped as DWV-A)²⁵. However, in the case of co-infection by DWV-A and DWV-B, DWV recombinants are produced²³, establishing a phylogenetic bridge between the two main genotypes²⁴. Genotypes DWV-A and DWV-B are most commonly detected²⁴⁻²⁷ whereas the most recently identified DWV genotypes, C and D, are considered rare or

<footnote n="1">¹ANSES Sophia Antipolis Laboratory, Unit of Honey Bee Pathology, Sophia Antipolis 06902, France. ²ANSES Ploufragan - Plouzané - Niort Laboratory, Unit of Viral Genetics and Biosecurity, Ploufragan 22440, France. ³LLH Bee Institute Kirchhain, 35274 Kirchhain, Germany. ⁴Yannick Blanchard and Richard Thiéry are jointly supervised this work. ✉email: eric.dubois@anses.fr</footnote>
extinct, respectively²⁸. Although DWV-C (or related recombinants), has been found in some countries, there is no clear data about its prevalence or impact on honey bee health²⁷⁻³⁰. Virulence analyses of the DWV-A and B genotypes have suggested that DWV-B might be more virulent in adult bees²²,²³,³¹, even though no significant differences in the prevalence of clinical symptoms or death were observed after inoculating DWV-A or DWV-B into honey bee pupae³²,³³. These results support previous studies suggesting that the overall prevalence of virulent viruses in honey bee brood could be reduced over the period of *V. destructor* mite infestation¹²,³⁴,³⁵. These field and experimental studies have shown that the introduction of the *V. destructor* mite changes the viral landscape, resulting in reduced prevalence of the most virulent viruses such as the sacbrood virus (SBV) and acute bee paralysis virus (ABPV), with its related viruses, the Kashmir bee virus (KBV) and Israeli acute paralysis virus (IAPV). The authors also suggested that a negative selection on *V. destructor* promoting the most virulent forms of DWV could impair pupae survival and, therefore, mite reproduction. Nevertheless, the main difference between DWV genotypes A and B is their ability to multiply in the *V. destructor* mite. Posada-Florez et al. (2019) observed that DWV-A did not replicate in the mite and that DWV-A spread by the mite was correlated with highly-efficient viral transmission by inoculation during mite feeding³⁶. In contrast, two studies have confirmed that DWV-B does replicate in the mite³⁷,³⁸. Gusachenko et al. (2020) investigated a recombinant genome (5' non-coding sequences and structural protein-coding sequences from DWV-B and non-structural protein-coding sequences from DWV-A) and suggested that the genomic organisation of several DWV recombinants could be a key factor for replication in both the mite and the honey bee³⁸. A recent study of the antiviral RNA interference pathway of *V. destructor* suggested that both DWV-A and DWV-B replicate in the mite but DWV-A infections may be more efficiently suppressed by mite antiviral defences³⁹. Moreover, this finding supports the enhanced capability of DWV-B and several other recombinant DWVs to be vectorised by the mite, amplifying the virus in salivary glands for inoculation above all into bee larvae³⁷. Since its first detection in Europe or the USA, the prevalence of DWV-B has increased over time⁴⁰. Ryabov et al. (2017) estimated that DWV-B spread within the United States between 2010 and 2016 ⁴¹. Nevertheless, DWV-A was still the most prevalent genotype in samples collected in 2016. In Europe, where DWV-B spread earlier, the dominance of DWV-B over DWV-A was reported in Germany between 2016 and 2017 ⁴², in the United Kingdom in 2017, and in Italy in 2019 ⁴⁰, and is considered to have taken place in 2014 in France⁴³. Therefore, the prevalence of DWV genotypes may differ among countries, resulting in an initial bottleneck by selection of the *Varroa*-adapted DWV variants followed by viral diversification, thus evading antiviral RNA interference⁴⁴. Natsopoulou et al. (2017) regretted the lack of knowledge regarding the European distribution of DWV genotypes⁴⁵.

In this study, we first evaluated the variations in viral load of DWV genotypes A and B at the colony level using pooled bee samples, and at the individual level using single bee heads. All the samples were collected between 2010 and 2017 from 15 different European countries (Fig. 1). We then compared the viral loads observed in single bee heads from honey bees either with or without deformed wings, and in *V. destructor* mites. We selected a subset of samples and performed RNA sequencing to study the genomic structure of the DWV genotypes. We characterised recombinant and parental complete consensus genomic sequences. We mapped the recombinant junctions and performed a phylogenetic study to evaluate the genetic relationship between the DWV genotypes obtained from these different European locations. Our study highlights the increasing importance of DWV-B and related recombinant viruses in the honey bee viral landscape.

## Results

### Viral prevalence and loads of DWV-A and DWV-B in the RNA collection (honey bee samples collected between 2010 and 2011)

Deformed wing virus (DWV) detection and quantification at the honey bee colony level were estimated by two RT-qPCRs, specific for either the A- or B-genotype VP3 coding sequences, performed by analysing the RNA collection prepared from pooled bee heads (n = 44; Supplementary Figure S1A, Supplementary Table S1). No DWV-B was detected in samples from Austria. DWV-A was detected more frequently or at higher mean loads than DWV-B in samples from North Macedonia and Bulgaria. On the contrary, DWV-B was detected more frequently and at higher loads than DWV-A in samples from Germany and Croatia. Among the samples from Germany, 14 were from colonies treated against *V. destructor* mites (oxalic acid). These 14 samples had significantly lower mite infestation levels (Mann-Whitney U test, *p* = 3.9 × 10⁻⁷) and either no detection or lower loads of DWV-A or DWV-B (max 4.8 log₁₀ copies/bee) than untreated colonies (max 8.7 log₁₀ copies/bee). Finally, we found a positive correlation (Spearman's test, *p* = 9.4 × 10⁻⁷) between the mite infestation level and the maximum viral load per bee head for DWV-A or DWV-B (Supplementary Figure S1B).

### Viral prevalence and loads of DWV-A and DWV-B in pooled samples of whole honey bees collected in Europe

The viral loads quantified in the pools of whole honey bee were analysed separately, as they were processed using different protocols (Fig. 2A and Supplementary Table S1).

DWV-B was not detected in any of the ten samples collected from the *V. destructor*-free European islands: Ushant (France; six colonies), Colonsay (Scotland; one colony) and Åland (Finland; three colonies). Only one sample from Ushant Island tested positive for DWV-A at a low viral load of 2.9 log₁₀ copies/bee.

The observed prevalence of DWV-A was 70% (n = 43) and that of DWV-B was 82% (n = 50) after testing 61 pooled samples collected in *V. destructor*-infested areas. Four distinct patterns were observed: samples containing both DWV-A and DWV-B (56%, n = 34); samples with only DWV-A (15%, n = 9); samples with only DWV-B (26%, n = 16); and samples (3%, n = 2) with no DWV genotypes detected. The positive samples (n = 59) were quantified at viral load levels ranging from 2.9 to 10.5 log₁₀ copies/bee for DWV-A, and from 5.4 to 11.9 log₁₀ copies/bee for DWV-B. The mean viral loads in the infected pools were 7.2 and 8.7 log₁₀ copies/bee for DWV-A and DWV-B, respectively. The DWV-B load was significantly (Mann-Whitney U test, *p* = 3.9 × 10⁻³) higher than
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Geographical location of the 116 colonies included in this study. Each dot represents a colony positioned according to the city closest to the apiary. Colonies from *V. destructor*-free and -infested areas are represented in blue and red dots, respectively.</figcaption></figure>

the DWV-A load. Additionally, we observed that the infected whole honey bees could be subdivided into low and high viral loads with thresholds set at 6.0 and 7.3 log₁₀ copies/bee for DWV-A and DWV-B, respectively (Supplementary Figure S2). Using these thresholds, we observed exclusively high DWV-A load levels in 16 samples (27%) and high DWV-B load levels in 22 samples (37%). Both genotypes were observed at high viral load levels in 14 samples (24%). Finally, seven samples (12%) showed low levels of both DWV genotypes.
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Colony-level viral loads of DWV-A and DWV-B estimated from honey bee pools. **A:** Variability of the observed viral loads. Each dot represents a pool of foragers from a colony according to DWV-A (x-axis) and DWV-B (y-axis) copies per bee. Colonies from mite-free or mite-infested areas are in blue and black dots, respectively. RT-qPCR quantifications were expressed in equivalents of viral genome copies per bee (copies per bee). **B:** Viral loads of DWV-A and DWV-B according to the country of origin. The red and blue bars represent the DWV-A and DWV-B loads, respectively. Viral loads were expressed in equivalents of viral genome copies per bee (copies per bee). Country codes are on the x-axis and initials are as follows: FR = France, UK = United Kingdom, ES = Spain, GR = Greece, HR = Croatia, MK = North Macedonia, RS = Republic of Serbia, RO = Romania, NO = Norway, IT = Italy, MD = Moldova.</figcaption></figure>
We investigated the relationship between DWV-A and DWV-B loads and the location of apiaries (Fig. 2B). DWV-B was detected in all the European countries included in this study. However, this genotype was detected at higher loads than DWV-A in samples from Greece, Spain, the United Kingdom and finally France, where no DWV-A was detected. In contrast, DWV-A was detected at higher loads than DWV-B in samples from Moldova, Italy, Norway and Romania. Samples collected in 2017 from Croatia, Macedonia or Serbia had similar DWV loads for both A and B genotypes.

## Viral load quantification in single honey bee heads and V. destructor mites

Both DWV-A and DWV-B viral loads were quantified by RT-qPCR on RNA extracted from the single heads of honey bees with either deformed or normal wings (Fig. 3A, Supplementary Figure S3, and Supplementary Table S1). The relative abundances of both viruses were consistent between honey bee heads and the viral load quantified in the pooled sample from the same colony. The viral loads found in single heads were distributed over a wider range than those in the pooled honey bees, specifically for the lower viral loads. Viral loads ranged from 1.8 to 10.7 and from 1.8 to 11.0 log₁₀ copies/head for DWV-A and DWV-B, respectively. The DWV-B load (mean = 7.4 log₁₀ copies/head) was significantly higher (Mann-Whitney U test, *p* = 1.5 × 10⁻⁸) than the DWV-A load (mean = 5.7 log₁₀ copies/head). We examined the relationship between the country of origin and the distribution of DWV genotype viral loads in honey bee heads (Supplementary Table S1). DWV-B viral loads

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 3.** DWV-A and DWV-B loads quantified in single honey bee heads and V. destructor. **A:** Global distribution of viral loads measured in single heads from honey bees with deformed (red dots) and normal wings (blue dots). Each dot represents one head of a single honey bee according to DWV-A (x-axis) and DWV-B (y-axis) load level. RT-qPCR was performed on RNA extracted from single honey bee heads. RT-qPCR quantification levels were expressed in equivalents of viral genome copies per bee head (copies per bee head). **B:** DWV-A or DWV-B load (in log₁₀) and maximum load from DWV-A and DWV-B quantified in single heads of honey bees with either normal or deformed wings. The box plots show the distribution of populations, with the first quartile (25%) and third quartile (75%) (box), the median (50%, white dash), the mean (cross), the minimal and maximal values (whiskers), and outliers (circle). ***: Mann-Whitney U test, *p* < 0.001.</figcaption></figure>
were higher than DWV-A loads in the heads of honey bees originating from France (four colonies), the United Kingdom (two colonies), Spain, Germany and Macedonia. In contrast, DWV-A was detected at higher load levels than DWV-B in samples from Moldova, Italy and Romania.

Taking into account the observed wing morphology of the tested honey bees, DWV-A or DWV-B loads were significantly higher (Mann-Whitney U test, *p* = 4.9 × 10⁻⁵ or *p* = 2.9 × 10⁻⁷, respectively) in the head of honey bees with deformed wings (*n* = 109; DWV-A mean load = 5.1 log₁₀ copies/head; DWV-B mean load = 7.2 log₁₀ copies/head) than in the head of honey bees with normal wings (*n* = 74; DWV-A mean load = 3.0 log₁₀ copies/head; DWV-B mean load = 4.4 log₁₀ copies/head) (Fig. 3B and Supplementary Table S1). The greater virus load in honey bees with deformed wings (mean = 9.3 log₁₀ copies/head) than in honey bees with normal wings (mean = 5.2 log₁₀ copies/head) was more striking when using the maximum virus load from DWV-A or DWV-B (Mann-Whitney U test, *p* = 0). Finally, the results showed that DWV-B viral loads were significantly higher than DWV-A viral loads in the infected heads of honey bees with either normal wings (Mann-Whitney U test, *p* = 3.5 × 10⁻³) or deformed wings (Mann-Whitney U test, *p* = 1.9 × 10⁻⁷).

Viral loads in bee heads and *V. destructor* mites were quantified for 11 colonies (Supplementary Figure S3 and supplementary Table S1). We observed that the relative proportion of DWV-A and DWV-B, as well as the viral loads, were concordant between *V. destructor* mites and single bee heads for five colonies: ES-Arbeiza21, FR-Caumont8, MK-MalinoSvetiNikole15, RO-TarguMures31 and UK-SherburnInElmet35.

## Honey bee virus sequencing

In order to obtain deeper insights into the variability of DWV genomes, we selected and sequenced 93 isolates. The reads were aligned against already known honey bee viruses (Fig. 4). DWV-B was the most prevalent virus (94%, *n* = 87), while DWV-A was identified in 74% (*n* = 69) of the samples. Moreover, by analysing the sequencing data, we were able to identify other honey bee viruses within our samples. We found ABPV (23%, *n* = 21), apis rhabdovirus type 1 (ARV-1; 26%, *n* = 24) and type 2 (ARV-2; 9%, *n* = 8), black queen cell virus (BQCV; 3%, *n* = 3), lake Sinai virus type 2 (LSV-2; 1%, *n* = 1) and SBV (11%, *n* = 10).

## Analysis of DWV genotypes

Next-generation sequencing (NGS) provided clear evidence of recombinant genomes containing fragments from both DWV-A and DWV-B, with a coverage depth around 4 log₁₀ (Fig. 5A and Supplementary Figure S4). DWV recombinants junctions were defined using LUMPY, further confirmed by ViReMa analysis, and the variation in both DWV-A and DWV-B reference genome coverage depth profiles (Fig. 5B, Supplementary Table S4 and Supplementary Table S2). LUMPY identified 85 recombinant junctions from 26 samples. Sixty-seven of these were also observed in ViReMa results. The average difference between the position junctions identified by LUMPY and by ViReMa was 79 and 60 nucleotides on DWV-A and DWV-B genomes, respectively. Among the identified junctions, 49 were used to guide the reconstruction of the recombinant genomes. Thirty one were located within a single contig aligning partially to DWV-A and DWV-B (Supplementary Table S2). Six main types of DWV recombinant profiles were identified (Rec-I to Rec-VI, Fig. 5A) based on the DWV-A sequence observed in the recombinant genomes. The DWV-A genomic fragments involved in the recombinant genomes often overlapped the 5'UTR, L-protein domain, non-structural domain and more rarely the 3'UTR domain. In contrast, DWV-B genomic fragments from the L-protein domain to the helicase were more often observed in recombinant genomes. Consequently, recombinant junctions were not randomly distributed and clustered in three genomic areas along DWV-A and DWV-B genomes (Fig. 5B). The first region is located between nucleotides 820–2201 of the DWV-A genome sequence. It overlaps with the 5'UTR, L-protein and VP2 domains. The second region was located at positions 4590–6003, mainly overlapping with the helicase domain. The third region was located at positions 9628–9891, overlapping with the end of the 3 C-pro/RdRp and 3'UTR domain. A few additional breakpoints were also found outside these regions (positions 2645, 7282 and 8547).

After sequencing RNA from honey bee pools from varroa-free areas, it was not possible to identify the recombinant junctions because the average DWV-A or B coverage was too low. Nevertheless, there was an increase in the coverage depth of the non-structural portion of the DWV-A genome in the pooled sample of honey bees from the Åland Islands, suggesting the presence of a DWV recombinant genome (Supplementary Figure S4).

We then investigated the relationship between the different DWV genotypes assembled from our samples (Supplementary Tables S3 and S4). We found a significant association between the absence of DWV-A and the detection of DWV recombinants (chi-square test, *p* = 7.9 × 10⁻⁶) or DWV-B (chi-square test, *p*-value = 1.6 × 10⁻³). An analysis of the pooled samples revealed that DWV-A was the main genotype in eight colonies from Italy, Moldova, the Republic of Serbia and Romania. DWV-B was the main genotype in nine colonies from Croatia, France, Greece, Macedonia, Spain and the United Kingdom. Both genotypes were also observed with equal coverage depth in two colonies, one from Croatia and the other from Italy. One DWV recombinant was observed to be dominant in three colonies from the same apiary in France. DWV recombinants and DWV-B genotypes were found to have the same coverage depth in colonies from France, Greece, Macedonia, Spain and the United Kingdom. Lastly, we also observed one colony from Italy where DWV-A, DWV-B and DWV recombinants were detected with equal coverage depth.

Results from recombinant junctions and variation in coverage depth profiles from single bee head sequencing data were similar to those found at the colony level (by testing pooled honey bee samples) (Fig. 5, Supplementary Figure S4). Several coverage depth profiles (e.g. ES-Arbeiza22_2017_H62) suggested the presence of multiple DWV recombinants co-existing in the same honey bee head. Genomes isolated from heads of honey bees with either deformed or normal wings were found in the entire spectrum of DWV genotypes (DWV-A, DWV-B and DWV recombinants) without any specific association.
<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Detection of honey bee virus reads by next-generation sequencing. The reads obtained from bee pools or single bee heads were aligned against the genomes of known honey bee viruses. The sequences of eight viral species were covered to above 80%. The darkness of the red indicates the coverage depth (scale from 0 to 15). The results are given by country: UK = United Kingdom, RS = Republic of Serbia, RO = Romania, MK = North Macedonia, MD = Moldova, IT = Italy, HR = Croatia, GR = Greece, FR = France, FI = Finland, ES = Spain, DE = Germany.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Genome profiles and recombinant junctions. **A:** Example of coverage depth and LUMPY profiles observed in six single bee heads. DWV-A and DWV-B coverage depth in log₁₀ scale calculated in 150 nucleotide windows along the genome are represented by red and blue graphs, respectively. Dots represent the allele balance of the break end identified by LUMPY and SVtyper for DWV-A (red) and DWV-B (blue). The main genomes deduced are symbolised below each graph. **B:** Positions of the recombination junctions along the DWV-A genome. Each red dot (located on the viral genome by the black line) represents a recombination junction identified in DWV recombinant. The genomic domains⁴⁶ of DWV-A are colour coded and ordered in the legend below the graph according to their position from the 5'-end to the 3'-end.</figcaption></figure>

## Phylogenetic analysis of DWV genomes

To investigate the phylogeographic relationships of the DWV genotypes, a neighbour-net splits graph (Fig. 6) and a maximum likelihood tree were constructed (Supplementary Figure S5). The neighbour-net splits graph showed a main split separating the DWV-A genomes from the reticulated topologies containing the recombinants and the DWV-B genomes. The branches within the DWV-A cluster were longer than in the DWV-B cluster, suggesting more diversity in the former. The DWV recombinant cloud was subdivided into two clusters: one clustering Rec-I, II and III genomes close to DWV-B, and the second clustering Rec-IV, V and VI in a median position between DWV-A and DWV-B. DWV recombinants clustering in similar groups were observed in samples from different countries. For example, recombinant genomes of type Rec-I were found in samples collected in North Macedonia and Spain, and recombinant genomes of type Rec-III were found in samples collected in France, Spain and the United Kingdom.
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Neighbour-net splits graph of DWV genotype genomes. The dashed line represents the main split. The scale bar represents the unit for the expected number of substitutions per site. Red, blue, orange and pink correspond to DWV-A, DWV-B, DWV recombinant and Kakugo virus genomes, respectively. Oval and rectangular shapes represent the reference and reconstructed genomes, respectively. The long, vertical branch connects to the DWV-C reference genome and was split for clarification. The different types of recombinant genomes are symbolised together with their position in the graph.</figcaption></figure>

In the maximum likelihood tree, we frequently found clusters of DWV genotypes from the same colony or the same apiary. Less frequently, we observed the clustering of genomes isolated from samples collected in neighbouring countries in the DWV-A group, such as Romania and Moldova (Supplementary Figure S5C). These genomes cluster with genomes isolated in Asian countries (China, South Korea and Japan). We also observed DWV genotype genomes from distant countries within the same cluster, such as DWV-B genomes isolated from samples collected in Germany and Croatia. As observed in the splits graph (Fig. 6), the branches from the common ancestor were shorter for DWV-B than DWV-A. The clustering of DWV variants from the same colony or the same apiary was strengthened by the phylogenetic analysis of partial sequences (Supplementary Figure S6). The analysis of VP1-VP3 coding sequences separated with more accuracy the DWV-A from the DWV-B or DWV recombinants. By analysing the 3 C-protease and RNA dependant RNA polymerase (3Cpro-RdRp) coding sequences, the DWV recombinants were located in both phylogenetic groups (DWV-A and DWV-B). In addition, these partial genome sequence analyses established the phylogenetic relationship between the parental and the recombinant DWV variants detected in the same colony or in the same honey bee head.

## Discussion

In this study, we compared the distribution of the deformed wing virus (DWV) genotypes A and B, and the viral load levels quantified by specific real-time PCR by analysing samples (n = 308) collected between 2010 and 2017 from 116 colonies in 15 European countries. The latest genotype to emerge, DWV-B, was reported for the first time in V. destructor mites collected in 2001 from the Netherlands⁹, but has already become widespread in Europe some 20 years later (74% of the tested honey bee pools). Except in the samples from Austria (collected from five colonies of the same apiary in 2010), DWV-B was detected in all the sampled countries, suggesting its rapid spread in Europe. Our results are supported by previous studies that have reported a higher prevalence of DWV-B in the United Kingdom³¹,⁴⁷, Germany²⁷,³¹,⁴⁰, France²⁴,⁴⁸ and Lithuania⁴⁹. Nevertheless, differences between countries were observed when comparing the viral loads of both DWV genotypes. DWV-B was quantified at higher loads than DWV-A in samples from France, the United Kingdom, Spain, Germany and Croatia. DWV-A was detected in 64% of the pooled samples. Based on viral quantification, we could not exclude that a better-adapted DWV-A still infects honey bees from Italy, Romania, Moldova, Bulgaria, North Macedonia and Austria. This could also suggest that DWV-A was dominant in these countries, because DWV-B might have been introduced only shortly before our study was conducted⁵⁰. For example, DWV-A was dominant in 2010–2011 in North Macedonia whereas DWV-B became dominant in 2017, although the apiaries used in this comparison were different (Fig. 2B, Supplementary Figure S1 and Supplementary Table S1). Recently, Paxton et al. (2022) suggested a switch in the prevalence of DWV-B over DWV-A in Italian honey bees between 2011 and 2019⁴⁰. Considering the DWV recombinants previously described in Europe and in the US, which contain in their genome the DWV-B VP3 coding sequence²²,²⁴,⁴⁴,⁵⁰⁻⁵², we expected that most of the DWV recombinants would be detected by our DWV-B specific qPCR. This inclusivity of qPCR for DWV-B or DWV recombinants was supported by our sequence analysis of the DWV genotypes. Even when new DWV recombinants were characterised, all of them contained the DWV-B VP1-VP3 coding sequence (Fig. 5 and Supplementary Figure S4). Therefore, the spatial distribution of the DWV genotypes observed in our study was the result of progressive replacement of DWV-A by DWV-B and/or DWV recombinants. This improved efficiency of DWV-B replication
in bees was recently demonstrated in a study using infectious plasmids of both DWV genotypes⁴⁴, which also revealed that co-infection of honey bees by both DWV genotypes leads to the production of DWV recombinants.

Ten pools of honey bees from colonies established in three different *V. destructor*-free areas were sequenced. DWV-A was detected at an extremely low level ($2.94 \log_{10}$ copies per bee) in one pool from a colony on Ushant Island. Previous studies report the finding of DWV genotypes prior to the arrival of *V. destructor* at low virus loads and with low colony-level prevalence in Hawaii, New Zealand and Newfoundland¹⁰,¹²,⁵³. However, a recent study on Australian honey bees reported the absence of DWV in the absence of the mite⁵⁴. One explanation for these differences might be the importation of DWV-infected queens and the horizontal or vertical transmission of the virus without any involvement of *V. destructor* mites⁶,⁵⁵. Unfortunately, the recent introduction of the *V. destructor* mite on Ushant Island⁵⁶ and in Australia may have changed the DWV loads and diversity.

Analysis of DWV genotype distribution in 105 pools (of whole bees or bee heads) from varroa-infested areas revealed that only DWV-A was detected in 13% of samples (n=14) and only DWV-B was detected in 24% of samples (n=25). However, we did not find any statistically significant differences in either proportion (chi-square test, p=0.05). Ryabov et al. (2017) previously found that colonies from the US with detectable DWV-B and undetectable DWV-A were significantly under-represented⁴¹. This discrepancy could be due to the target sequence of the primers used by Ryabov, et al.⁴¹ (non-structural protein coding sequence: 8,645 to 8,780 base on DWV genome, AY292384) being different from that of the primers in our study (VP3 coding sequence: 4,245 to 4,296 base on DWV genome, AY292384). Since the genomic structure of the previously described viral recombinants contained DWV-A non-structural protein domains and DWV-B structural protein domains²²,²⁴, we considered that recombinants were detected by the DWV-A qPCR in the Ryabov, et al.⁴¹ study and by the DWV-B qPCR in our study. Therefore, DWV recombinants may modify the distribution of DWV-A and DWV-B viral loads according to the targeted genome sequence. Considering that the 5'NC and structural protein-coding sequences from DWV-B could be associated with the ability to multiply in the *V. destructor* mite³⁸, it would make sense to assume that these recombinants originate from DWV-B rather than DWV-A. The assembly of sequencing reads obtained from a single bee head is much simpler than that of reads from pools of whole bees for reconstructing DWV genomes. The pools contained DWV sequence mixtures that were difficult to identify using our pipeline (Supplementary Figure S8). Nevertheless, we detected co-infections by different DWV genotypes and recombinant DWV even in bee heads. However, we detected the presence of a recombinant virus without association with a DWV-A or a DWV-B by qPCR and sequencing from two bee heads (FR-Caumont6_2017_H78 and H111; Supplementary Table S1 and Supplementary Figure S4). In addition to previous observations made under experimental conditions showing that recombinant DWVs could overcome DWV-A⁵⁷, our results from natural infections suggested the ability of some DWV recombinants to spread from bee to bee independently of DWV-A or DWV-B.

Our RT-qPCR data showed that even if one DWV genotype was not detected in the pool, it was detected at a low level in the head of honey bees from the same colony and even in *V. destructor* mites (Fig. 3 and Supplementary Table S1). This suggests that the sensitivity of quantification for individual samples (honey bee head or *V. destructor* mite) was higher than for the pool made from whole honey bees. This conclusion is also supported by the viral loads quantified in several heads, that were below the limit of detection we had previously reported for DWV-A and DWV-B in pools (3.4 and 5.0 $\log_{10}$ copies/honey bee, respectively⁴³). Indeed, the matrix complexity and the abundance of bee contaminants are higher in the pool of whole honey bees than in single honey bee heads and may reduce the yield of viral RNA extracted or reduce the sensitivity of RNA detection. Nevertheless, the representativeness of the bee sample must also be considered, and at least 60 bees should be collected to detect an intra-colony infection with a prevalence limit of 5% (with 95% of confidence). Thus, testing a pool of honey bee heads may be more useful than testing a pool of whole honey bees for DWV genotype detection in colonies. Therefore, both DWV-A and DWV-B were detected in the head of symptomatic and asymptomatic honey bees as reported previously²⁰. The obvious difference between symptomatic and asymptomatic honey bees was the viral load level. In symptomatic honey bees, the viral load was calculated to be above $6.5 \log_{10}$ and $7.7 \log_{10}$ virus/head for DWV-A and DWV-B, respectively. This is consistent with previously reported thresholds separating covert and overt infections, and suggests that a high quantity of viral genome copies is associated with wing deformities, regardless of the DWV genotype³³,⁴³,⁵⁸. Therefore, an RT-qPCR method quantifying both DWV genotypes, such as the one recently developed⁵⁹, would not be sufficient to establish which DWV genotype is associated with overt infection. We observed that the relative proportions of DWV-A and DWV-B as well as the viral loads detected in *V. destructor* mites and in honey bee heads from the same colony were correlated (Supplementary Figure S3). Nevertheless, the virus detected in individual mites could have been viruses either accumulated by feeding³⁶ or by replicating in mite salivary glands³⁷. A deeper analysis of the DWV genotypes present and replicating in the mite would clarify the dynamics of DWV genotype transmission by *V. destructor*.

In conclusion, it was suggested that the *V. destructor* mite was modifying the prevalence of DWV genotypes by vectorising better-adapted viruses (genotype B and DWV recombinants). Our study provides new sequences of DWV-A, DWV-B and DWV recombinants detected in honey bees from 15 European countries exposed to the mite. New hot points for genomic recombination were identified and support the plasticity of the viral genome. All DWV genotypes could be detected in the heads of honey bees, whether they had deformed or normal wings. The main characteristic distinguishing DWV-A from DWV-B and DWV recombinants was the viral loads. The highest levels of DWV-B or DWV recombinants detected in honey bees and *V. destructor* are a robust characteristic reported by several studies. Our study did not allow us to conclude on the date of emergence of DWV-B and associated DWV recombinants. However, it suggests that the emerging DWV-B genotype and related DWV recombinants are continuing to spread in Europe and might replace the historical DWV-A genotype.
# Methods

## Sample collection

In 2017, we collected samples from 116 colonies across 15 European countries (Fig. 1). Complementary samples from viral collections were provided as purified RNA from pools of 10 honey bee heads (44 colonies) or as whole adult honey bees (72 colonies). The detection of DWV in total head RNA was suggested to be a molecular marker for overt DWV infection<sup>20,60,61</sup>. Samples were received frozen (in dry ice) in our laboratory and were kept at -80 °C before analysis. Details of the samples used in this study are listed in Supplementary Table S6.

## Viral genomic RNA preparation

To analyse pooled samples, seven to 20 whole honey bees were crushed in 0.01 mM phosphate buffer (PB, 1 mL per bee) with 5- to 8-mm tungsten beads using a Mixer Mill MM 400 (RETSCH) for 3 min at 30 Hz. Bee homogenates were then cleared by two successive centrifugations at 8,000 x g for 10 min at 4 °C. Two millilitres of supernatant were used to perform polyethylene glycol precipitation of macromolecules (10% PEG 8,000 and 300 mM NaCl, final concentrations) after incubation for 2 h at 4 °C. Viruses were recovered by centrifugation at 10,000 x g for 30 min at 4 °C. Each pellet was suspended in 400 µL of phosphate buffer (PB). To analyse single bee heads, each head was crushed in PB (500 µL per bee head) with three 3-mm tungsten beads using a Mixer Mill MM 400 (RETSCH) for 1 min at 30 Hz. The homogenate was cleared by two successive centrifugations at 8,000 x g for 10 min at 4 °C.

Unprotected nucleic acids co-purified with viruses from the pools of whole bees or bee heads were cleaved in 1X DNase reaction buffer (Invitrogen, AM2238) by the Turbo nuclease (25 U/ml, Sigma-Aldrich T4330) at 37 °C for 2 h.

Individual *V. destructors* were crushed in PB (500 µL per mite) with 1 mL Ten Broeck homogeniser. The homogenates were incubated on ice for 10 min before viral RNA extraction.

RNA was extracted from 140 µL of viral suspension with the QIAamp Viral RNA Mini Kit (Qiagen), without RNA carrier, according to the manufacturer's instructions. Purified RNA from bee pool or bee head was recovered in 60 µL of elution buffer (80 µL for purified RNA from *V. destructor* mite).

## RT-qPCR analysis

Retro-transcription into cDNA was performed using random primers and the SuperScript II reverse transcriptase (SSRT) kit (Invitrogen) from 12.5 µL of purified RNA according to a previous report (final volume 20 µl)<sup>62</sup>. DWV-A and DWV-B were quantified (using 5 µL of complementary DNA) by two genotype-specific real-time PCRs (qPCR) amplifying the VP3-coding sequences (final volume 25 µl; Supplementary Table S6)<sup>43</sup>. Before log₁₀-transformation, the RT-qPCR results were converted into copies per bee, bee head or *V. destructor* mite taking into account the volume of homogenate used for RNA purification, the volume of the elution buffer used for RNA recovery from the spin column, the volume of purified RNA used for cDNA synthesis and the volume of cDNA tested by qPCR. The final quantities were calculated by multiplying the instrument's measurement with the dilution factor: 27.4 for pool load, 68.5 for bee head load or 91.4 for mite load. The quantification method's accuracy was also analysed (Supplementary Figure S7<sup>63</sup>). Both RT-qPCR methods quantified the recombinant plasmids within the acceptability limits, i.e. ±1.0 log₁₀ copies per whole bee or per bee head.

## Statistical analysis

The number of DWV-A or DWV-B log₁₀ copies estimated from the heads of honey bees either with or without deformed wings was compared using the Mann-Whitney U test (used for comparing not normally distributed viral loads of equal or different sample sizes). The relationship between the viral loads in bees and of *V. destructor* infestation level, were tested with Spearman's test for small numbers of data points and non-normal distributions. The DWV-A and DWV-B prevalences were compared by a chi-square test. The p-values were calculated using R package lawstat (v3.2). Differences were considered significant at *p* < 0.05.

## RNA sequencing and data pre-processing

The NEBNext rRNA depletion kit (New England Biolabs, E6310) and NEBNext Ultra II directional RNA kit (New England Biolabs, E7760) were used to prepare the sequencing libraries using 12 µL of RNA eluate. The libraries were subject to quality control using a Tapestation HS D5000 (Agilent). Qubit (Invitrogen) was used for quantification. Libraries were pooled and sequenced with a Nextseq 500 (Illumina) using the 150-cycle high output kit at the sequencing core facility of the Brain & Spine Institute, Paris, France.

Sequencing reads were cleaned with Trimmomatic (Bolger et al., 2014, v0.36) applying parameter settings (ILLUMINACLIP: oligos.fasta: 2:30:5:1: true; LEADING: 3; TRAILING: 3; MAXINFO: 40:0.2; MINLEN: 36).

## Detection of honey bee viruses

Cleaned reads were aligned with 13 viral honey bee reference genomes (Supplementary Table S7) using Megablast<sup>64</sup>. Reads aligning more than once were removed from the output files. When more than 90% of a given viral genome was covered in a given sample, the reads per kilobase million (RPKM) was calculated using viral genome length and total number of cleaned reads produced as normalisation values. The results were visualised using plothmap in R (v1.0.12).

## DWV genotype genomic structure

For recombinant junction identification, cleaned reads were aligned simultaneously with the most relevant DWV-A and DWV-B parental genomes using Speedseq (v0.1.2)<sup>65</sup> (Supplementary Table S8). Then LUMPY (v0.2.13)-smoove (v0.2.5) was used to identify the recombination junctions from discordant and split reads<sup>66</sup>. The number of sampled reads used to estimate mean and standard deviation of insert size was adjusted to
optimise the analysis run time. Outputs from LUMPY were processed with SVtyper (v0.7.0) to obtain additional structural genotype statistics⁶⁵. Potential recombinant junctions were defined from break end events with over ten supporting reads (SU), over 1,000 X coverage depth (DP) and over 0.05 allele balance (AB; fraction of observations from alternate allele). A graphical representation of genome coverage depth and the position of complete LUMPY-defined break end events for each sample was generated with the ggplot2 (v3.1) package in R. Graphical representations of all LUMPY-defined junctions were generated using the trackViewer (v1.18) package in R.

## De novo assembly and genome reconstruction of DWV genotypes

Honey bee associated reads were removed after mapping to the Amel 4.5 reference genome (GCA_000002195.1). Reads were mapped to the most homologous DWV-A or DWV-B genotype reference genome (Supplementary Table S8) to estimate coverage depth. Average coverage depth was used to calculate the number of reads to sample in order to get an average coverage depth of 80X. Subsample reads were *de novo* assembled using SPAdes (v3.10.0) into contigs and aligned to the reference genome to generate a consensus sequence⁶⁷. This subsampling strategy was used to optimise the run time of *de novo* assembly and facilitate assembly of the main genotype. The alignment of consensus and contigs against the reference genome were visualised with Integrative Genomics Viewer (IGV, v2.3.91)⁶⁸. The assembly of a full-length sequence was achieved by using overlapping junctions. The graphical representations containing the coverage depth and LUMPY- and ViReMa⁶⁹-defined junctions were used to guide genome reconstruction specifically for recombinant genomes (Supplementary Figure S8). In the case of a mix of DWV-A or DWV-B and DWV recombinant, assembly of the recombinant genome was achieved by joining the contigs of the main background genome (DWV-A or DWV-B) with remaining contigs. The reconstructed genome were validated by alignment of the reads with BWA-MEME (v0.7.17) on the proposed genome⁷⁰. The putative junctions were particularly scrutinized and must display a perfect alignment of the reads on the reconstructed genome to be ascertained. Finally, the open reading frames of each reconstructed genome were validated by comparing them with the known start and end positions of the reference genomes.

## Phylogenetic analysis

Whole-genome nucleic acid sequences of reconstructed and reference genomes (Supplementary Table S9) were aligned using MAFFT (v7.221.3) in automatic mode with Galaxy. The alignment block was edited with MEGA (v10.0.5) by removing nucleotides at the start (110 nucleotides) and end (65 nucleotides). The alignment length was 10,066 nucleotides. The splits graph network was generated with Splitstree 4 using multiple sequence alignments. The network was constructed using uncorrected *p* distances transformed using the neighbour net method. For the phylogenetic tree, best-fit models of nucleotide substitutions were statistically selected using JmodelTest (v2.1.10) with default settings (ML optimised base tree and Best base tree search). GTR+G+I was selected as the best-fit substitution model according to Akaike criteria. Maximum likelihood trees were constructed using raxmlGUI v1.5 (RAxML v8.1.2) without outgroups using the ML+rapid bootstrap analysis type, with the number of bootstrap replicates set to autoMRE. The bootstrap convergence test defines 600 replicates as sufficient to obtain stable support values. The RaxML partition tree was edited with Dendroscope (v3.5.10).

Partial coding sequences (VP1-VP3 and 3Cpro-RdRp coding sequences, nucleotide position 2594 to 4616 and 7679 to 9821 on the DWV genome Y292384, respectively) were extracted from each of the 107 assembled consensus and 26 reference DWV genome sequences using BLAST. The Multiple Sequence Alignment Software MAFF was used in auto mode with default parameters to align the VP1-VP3 or 3 C-RdRp sequences⁷¹. Each alignment was trimmed with TrimAl using nogap mode⁷². The trimmed alignment Fasta file was loaded in the PhyML 3.1/3.0 aLRT from the Phylogeny.fr platform using aLRT statistical test and GTR substitution model⁷³. The phylogenetic tree was then edited with TreeDyn 198.3⁷³.

## Data availability

All the NGS data files used in this paper are publicly available in the NCBI's short read archive (SRA), which is accessible under BioProject ID: PRJNA1055031.

Received: 19 March 2024; Accepted: 15 January 2025

Published online: 08 February 2025

## References

1. Cornman, R. S. et al. Pathogen webs in collapsing honey bee colonies. *PLoS One*. 7, e43562. https://doi.org/10.1371/journal.pone.0043562 (2012).

2. de Miranda, J. R. & Genersch, E. Deformed wing virus. *J. Invertebr Pathol.* **103** (Suppl 1), 48–61. https://doi.org/10.1016/j.jip.2009.06.012 (2010).

3. Wilfert, L. et al. Deformed wing virus is a recent global epidemic in honeybees driven by Varroa mites. *Science* **351**, 594–597. https://doi.org/10.1126/science.aac9976 (2016).

4. Hasegawa, N. et al. Evolutionarily diverse origins of deformed wing viruses in western honey bees. *Proc. Natl. Acad. Sci. U S A.* **120**, e2301258120. https://doi.org/10.1073/pnas.2301258120 (2023).

5. de Miranda, J. R. & Fries, I. Venereal and vertical transmission of deformed wing virus in honeybees (*Apis mellifera* L). *J. Invertebr Pathol.* **98**, 184–189. https://doi.org/10.1016/j.jip.2008.02.004 (2008).

6. Yue, C., Schroder, M., Gisder, S. & Genersch, E. Vertical-transmission routes for deformed wing virus of honeybees (*Apis mellifera*). *J. Gen. Virol.* **88**, 2329–2336. https://doi.org/10.1099/vir.0.83101-0 (2007).

7. Mazzei, M. et al. Infectivity of DWV associated to flower pollen: experimental evidence of a horizontal transmission route. *PLoS One*. **9**, e113448. https://doi.org/10.1371/journal.pone.0113448 (2014).

8. Chen, Y. P., Pettis, J. S., Collins, A. & Feldlauber, M. F. Prevalence and transmission of honeybee viruses. *Appl. Environ. Microbiol.* **72**, 606–611. https://doi.org/10.1128/AEM.72.1.606-611.2006 (2006).
9. Ongus, J. R. et al. Complete sequence of a picorna-like virus of the genus Iflavirus replicating in the mite Varroa destructor. *J. Gen. Virol.* **85**, 3747–3755. https://doi.org/10.1099/vir.0.80470-0 (2004).

10. Martin, S. J. et al. Global honey bee viral landscape altered by a parasitic mite. *Science* **336**, 1304–1306. https://doi.org/10.1126/science.1220941 (2012).

11. Bowen-Walker, P. L., Martin, S. J. & Gunn, A. The transmission of deformed wing virus between honeybees (*Apis mellifera* L.) by the ectoparasitic mite varroa jacobsoni oud. *J. Invertebr Pathol.* **73**, 101–106. https://doi.org/10.1006/jipa.1998.4807 (1999).

12. Mondet, F., de Miranda, J. R., Kretzschmar, A., Le Conte, Y. & Mercer, A. R. On the front line: quantitative virus dynamics in honeybee (*Apis mellifera* L.) colonies along a new expansion front of the parasite Varroa destructor. *PLoS Pathog.* **10**, e1004323. https://doi.org/10.1371/journal.ppat.1004323 (2014).

13. Benaets, K. et al. Covert deformed wing virus infections have long-term deleterious effects on honeybee foraging and survival. *Proc. Biol. Sci.* **284** https://doi.org/10.1098/rspb.2016.2149 (2017).

14. Wells, T. et al. Flight performance of actively foraging honey bees is reduced by a common pathogen. *Environ. Microbiol. Rep.* **8**, 728–737. https://doi.org/10.1111/1758-2229.12434 (2016).

15. Iqbal, J. & Mueller, U. Virus infection causes specific learning deficits in honeybee foragers. *Proc. Biol. Sci.* **274**, 1517–1521. https://doi.org/10.1098/rspb.2007.0022 (2007).

16. Gisder, S. & Genersch, E. Viruses of commercialized insect pollinators. *J. Invertebr Pathol.* **147**, 51–59. https://doi.org/10.1016/j.jip.2016.07.010 (2017).

17. Dainat, B., Evans, J. D., Chen, Y. P., Gauthier, L. & Neumann, P. Dead or alive: deformed wing virus and Varroa destructor reduce the life span of winter honeybees. *Appl. Environ. Microbiol.* **78**, 981–987. https://doi.org/10.1128/AEM.06537-11 (2012).

18. Annoscia, D. et al. Haemolymph removal by Varroa mite destabilizes the dynamical interaction between immune effectors and virus in bees, as predicted by Volterra's model. *Proc. Biol. Sci.* **286**, 20190331. https://doi.org/10.1098/rspb.2019.0331 (2019).

19. Nazzi, F. et al. Synergistic parasite-pathogen interactions mediated by host immunity can drive the collapse of honeybee colonies. *PLoS Pathog.* **8**, e1002735. https://doi.org/10.1371/journal.ppat.1002735 (2012).

20. Zioni, N., Soroker, V. & Chejanovsky, N. Replication of Varroa destructor virus 1 (VDV-1) and a Varroa destructor virus 1-deformed wing virus recombinant (VDV-1-DWV) in the head of the honey bee. *Virology* **417**, 106–112. https://doi.org/10.1016/j.virol.2011.05.009 (2011).

21. Di Prisco, G. D. et al. Dynamics of persistent and acute deformed wing virus infections in honey bees, *Apis mellifera*. *Viruses* **3**, 2425–2441. https://doi.org/10.3390/v3122425 (2011).

22. Moore, J. et al. Recombinants between deformed wing virus and Varroa destructor virus-1 may prevail in Varroa destructor-infested honeybee colonies. *J. Gen. Virol.* **92**, 156–161. https://doi.org/10.1099/vir.0.025965-0 (2011).

23. Wood, G. R., Burroughs, N. J., Evans, D. J. & Ryabov, E. V. Error correction and diversity analysis of population mixtures determined by NGS. *PeerJ* **2**, e645. https://doi.org/10.7717/peerj.645 (2014).

24. Dalmon, A. et al. Evidence for positive selection and recombination hotspots in deformed wing virus (DWV). *Sci. Rep.* **7**, 41045. https://doi.org/10.1038/srep41045 (2017).

25. Mordecai, G. J., Wilfert, L., Martin, S. J., Jones, I. M. & Schroeder, D. C. Diversity in a honey bee pathogen: first report of a third master variant of the deformed Wing Virus quasispecies. *ISME J.* **10**, 1264–1273. https://doi.org/10.1038/ismej.2015.178 (2016).

26. Brettell, L. E. et al. A comparison of deformed Wing Virus in Deformed and Asymptomatic Honey bees. *Insects* **8** https://doi.org/10.3390/insects8010028 (2017).

27. Kevill, J. L., Highfield, A., Mordecai, G. J., Martin, S. J. & Schroeder, D. C. ABC Assay: Method Development and Application to quantify the role of three DWV Master variants in Overwinter colony losses of European Honey bees. *Viruses* **9** https://doi.org/10.3390/v9110314 (2017).

28. de Miranda, J. R. et al. Cold case: the disappearance of Egypt bee virus, a fourth distinct master strain of deformed wing virus linked to honeybee mortality in 1970's Egypt. *Virol. J.* **19**, 12. https://doi.org/10.1186/s12985-022-01740-2 (2022).

29. Brettell, L. E., Schroeder, D. C. & Martin, S. J. RNAseq analysis reveals virus diversity within Hawaiian Apiary insect communities. *Viruses* **11**, (2019). https://doi.org/10.3390/v11050397

30. Lopes, A. R., Low, M., Martin-Hernandez, R., Pinto, M. A. & De Miranda, J. R. Origins, diversity, and adaptive evolution of DWV in the honey bees of the azores: the impact of the invasive mite Varroa destructor. *Virus Evol.* **10**, veae053. https://doi.org/10.1093/ve/veae053 (2024).

31. McMahon, D. P. et al. Elevated virulence of an emerging viral genotype as a driver of honeybee loss. *Proc. Biol. Sci.* **283** https://doi.org/10.1098/rspb.2016.0811 (2016).

32. Tehel, A. et al. The two prevalent genotypes of an emerging infectious disease, deformed wing virus, cause equally low pupal mortality and equally high wing deformities in host honey bees. *Viruses* **11**, 114. https://doi.org/10.3390/v11020114 (2019).

33. Dubois, E. et al. Outcomes of honeybee pupae inoculated with deformed wing virus genotypes a and B. *Apidologie* **51**, 18–34. https://doi.org/10.1007/s13592-019-00701-z (2020).

34. Remnant, E. J., Mather, N., Gillard, T. L., Yagound, B. & Beekman, M. Direct transmission by injection affects competition among RNA viruses in honeybees. *Proc. R. Soc. B* 286, 20182452. (2019). https://doi.org/10.6084/m9.fishshare.c.4365839

35. Schroeder, D. C. & Martin, S. J. Deformed wing virus: the main suspect in unexplained honeybee deaths worldwide. *Virulence* **3**, 589–591. https://doi.org/10.4161/viru.22219 (2012).

36. Posada-Florez, F. et al. Deformed wing virus type A, a major honey bee pathogen, is vectored by the mite Varroa destructor in a non-propagative manner. *Sci. Rep.* **9**, 12445. https://doi.org/10.1038/s41598-017-47447-3 (2019).

37. Gisder, S. & Genersch, E. Direct evidence for infection of Varroa destructor mites with the Bee-Pathogenic deformed wing virus variant B - but not variant A - via fluorescence-in situ-hybridization analysis. *J. Virol.* **95** https://doi.org/10.1128/JVI.01786-20 (2021).

38. Gusachenko, O. N. et al. Green bees: reverse genetic analysis of deformed wing virus transmission, replication, and Tropism. *Viruses* **12** https://doi.org/10.3390/v12050532 (2020).

39. Damayo, J. E. et al. Virus replication in the honey bee parasite, Varroa destructor. *J. Virol.* **97**, e0114923. https://doi.org/10.1128/jvi.01149-23 (2023).

40. Paxton, R. J. et al. Epidemiology of a major honey bee pathogen, deformed wing virus: potential worldwide replacement of genotype A by genotype B. *Int. J. Parasitol. Parasites Wild* **18**, 157–171. https://doi.org/10.1016/j.ijppaw.2022.04.013 (2022).

41. Ryabov, E. V. et al. Recent spread of Varroa destructor virus-1, a honey bee pathogen, in the United States. *Sci. Rep.* **7**, 17447. https://doi.org/10.1038/s41598-017-17802-3 (2017).

42. D'Alvise, P., Seeburger, V., Gihring, K., Kieboom, M. & Hasselmann, M. Seasonal dynamics and co-occurrence patterns of honey bee pathogens revealed by high-throughput RT-qPCR analysis. *Ecol. Evol.* **9**, 10241–10252. https://doi.org/10.1002/ece3.5544 (2019).

43. Schurr, F. et al. Validation of quantitative real-time RT-PCR assays for the detection of six honeybee viruses. *J. Virol. Methods.* **270**, 70–78. https://doi.org/10.1016/j.jviromet.2019.04.020 (2019).

44. Ryabov, E. V. et al. Dynamic evolution in the key honey bee pathogen deformed wing virus: novel insights into virulence and competition using reverse genetics. *PLoS Biol.* **17**, e3000502. https://doi.org/10.1371/journal.pbio.3000502 (2019).

45. Natsopoulou, M. E. et al. The virulent, emerging genotype B of deformed wing virus is closely linked to overwinter honeybee worker loss. *Sci. Rep.* **7**, 5242. https://doi.org/10.1038/s41598-017-05596-3 (2017).

46. Lamp, B. et al. Construction and rescue of a molecular clone of deformed Wing Virus (DWV). *PLoS One.* **11**, e0164639. https://doi.org/10.1371/journal.pone.0164639 (2016).
47. Kevill, J. L. et al. DWV-A Lethal to Honey Bees (*Apis mellifera*): A Colony Level Survey of DWV Variants (A, B, and C) in England, Wales, and 32 States across the US. *Viruses* 11, (2019). https://doi.org/10.3390/v11050426

48. Manley, R. et al. Knock-on community impacts of a novel vector: spillover of emerging DWV-B from Varroa-infested honeybees to wild bumblebees. *Ecol. Lett.* **22**, 1306–1315. https://doi.org/10.1111/ele.13323 (2019).

49. Amsjektj-Graziani, P. et al. Molecular characterization and phylogenetic analysis of honeybee (*Apis mellifera*) Mite-Borne Pathogen DWV-A and DWV-B isolated from Lithuania. *Microorganisms* **12** https://doi.org/10.3390/microorganisms12091884 (2024).

50. Norton, A. M. et al. Adaptation to vector-based transmission in a honeybee virus. *J. Anim. Ecol.* **90**, 2254–2267. https://doi.org/10.1111/j.1365-2656.13493 (2021).

51. Wang, H. et al. Sequence recombination and conservation of Varroa destructor virus-1 and deformed wing virus in field collected honey bees (*Apis mellifera*). *PLoS One.* **8**, e74508. https://doi.org/10.1371/journal.pone.0074508 (2013).

52. Barroso-Arévalo, S., Vicente-Rubiano, M., Molero, F., Puerta, F. & Sanchez-Vizcaino, J. M. Nucleotide sequence variations may be associated with virulence of deformed wing virus. *Apidologie* **50**, 482–496. https://doi.org/10.1007/s13592-019-00660-5 (2019).

53. Shutler, D. et al. Honey bee *Apis mellifera* parasites in the absence of Nosema ceranae fungi and Varroa destructor mites. *PLoS One.* **9**, e98599. https://doi.org/10.1371/journal.pone.0098599 (2014).

54. Roberts, J. M. K., Anderson, D. L. & Durr, P. A. Absence of deformed wing virus and Varroa destructor in Australia provides unique perspectives on honeybee viral landscapes and colony losses. *Sci. Rep.* **7**, 6925. https://doi.org/10.1038/s41598-017-0729-0-w (2017).

55. Chen, Y., Evans, J. & Feldlaufer, M. Horizontal and vertical transmission of viruses in the honey bee, *Apis mellifera*. *J. Invertebr. Pathol.* **92**, 152–159 (2006).

56. Ménage, A. & L'Hostis, M. Ouessant colonisée par *Varroa destructor* ! *La. Santé de l'Abeille* **306**, 565–571 (2021).

57. Ryabov, E. V. et al. A virulent strain of deformed wing virus (DWV) of honeybees (*Apis mellifera*) prevails after Varroa destructor-mediated, or in vitro, transmission. *PLoS Pathog.* **10**, e1004230. https://doi.org/10.1371/journal.ppat.1004230 (2014).

58. Tentcheva, D. et al. Comparative analysis of deformed wing virus (DWV) RNA in *Apis mellifera* and *Varroa destructor*. *Apidologie* **37**, 41–50. https://doi.org/10.1051/apido:2005057 (2006).

59. Bradford, E. L., Christie, C. R., Campbell, E. M. & Bowman, A. S. A real-time PCR method for quantification of the total and major variant strains of the deformed wing virus. *PLoS One.* **12**, e0190017. https://doi.org/10.1371/journal.pone.0190017 (2017).

60. Mockel, N., Gisder, S. & Genersch, E. Horizontal transmission of deformed wing virus: pathological consequences in adult bees (*Apis mellifera*) depend on the transmission route. *J. Gen. Virol.* **92**, 370–377. https://doi.org/10.1099/vir.0.025940-0 (2011).

61. Yue, C. & Genersch, E. RT-PCR analysis of deformed wing virus in honeybees (*Apis mellifera*) and mites (*Varroa destructor*). *J. Gen. Virol.* **86**, 3419–3424. https://doi.org/10.1099/vir.0.81401-0 (2005).

62. Schurr, F. et al. Trueness and precision of the real-time RT-PCR method for quantifying the chronic bee paralysis virus genome in bee homogenates evaluated by a comparative inter-laboratory study. *J. Virol. Methods.* **248**, 217–225. https://doi.org/10.1016/j.jvir. omet.2017.07.013 (2017).

63. Piou, V., Schurr, F., Dubois, E. & Vetillard, A. Transmission of deformed wing virus between *Varroa destructor* foundresses, mite offspring and infested honey bees. *Parasit. Vectors.* **15**, 333. https://doi.org/10.1186/s13071-022-05463-9 (2022).

64. Morgulis, A. et al. Database indexing for production MegaBLAST searches. *Bioinformatics* **24**, 1757–1764. https://doi.org/10.1093 /bioinformatics/btn322 (2008).

65. Chiang, C. et al. SpeedSeq: ultra-fast personal genome analysis and interpretation. *Nat. Methods.* **12**, 966–968. https://doi.org/10.1 038/nmeth.3505 (2015).

66. Layer, R. M., Chiang, C., Quinlan, A. R. & Hall, I. M. LUMPY: a probabilistic framework for structural variant discovery. *Genome Biol* **15**, R84, (2014). https://doi.org/10.1186/gb/2014-15-6-r84

67. Bankevich, A. et al. SPAdes: a new genome assembly algorithm and its applications to single-cell sequencing. *J. Comput. Biol.* **19**, 455–477. https://doi.org/10.1089/cmb.2012.0021 (2012).

68. Thorvaldsdottir, H., Robinson, J. T. & Mesirov, J. P. Integrative Genomics Viewer (IGV): high-performance genomics data visualization and exploration. *Brief. Bioinform.* **14**, 178–192. https://doi.org/10.1093/bib/bbs017 (2013).

69. Sotcheff, S. et al. ViReMa: a virus recombination map of next-generation sequencing data characterizes diverse recombinant viral nucleic acids. *Gigascience* **12** https://doi.org/10.1093/gigascience/giad009 (2023).

70. Li, H. & Durbin, R. Fast and accurate short read alignment with Burrows-Wheeler transform. *Bioinformatics* **25**, 1754–1760. https://doi.org/10.1093/bioinformatics/btp324 (2009).

71. Katoh, K. & Standley, D. M. MAFFT multiple sequence alignment software version 7: improvements in performance and usability. *Mol. Biol. Evol.* **30**, 772–780. https://doi.org/10.1093/molbev/mst010 (2013).

72. Capella-Gutierrez, S., Silla-Martinez, J. M. & Gabaldon, T. trimAl: a tool for automated alignment trimming in large-scale phylogenetic analyses. *Bioinformatics* **25**, 1972–1973. https://doi.org/10.1093/bioinformatics/btp348 (2009).

73. Dereeper, A. et al. Phylogeny.fr: robust phylogenetic analysis for the non-specialist. *Nucleic Acids Res.* **36**, W465–469. https://doi.o rg/10.1093/nar/gkn180 (2008).

Acknowledgements

We are grateful to the Smartbees partners who supported this work by collecting honey bees and mites in their apiary: Elke Genersch and Sebastian Gisder (Länderinstitut für Bienenkunde e.V., Germany), Andone Estonba and Egoitz Galarza (University of the Basque Country, Spain), Eva Forsgreen (Swedish University of Agricultural Sciences, Sweden), Benjamin Basso (Institut technique et scientifique de l'apiculture et de la pollinisation, France), Olivier Lambert and Monique L'Hostis (Ecole Nationale Vétérinaire, Agroalimentaire et de l'Alimentation de Nantes-Atlantique, France), Nicolas Cougoule (French Aagency for Food, Environmental and Occupational Health & Safety, France), Fani Hatjina (Agricultural University of Athens, Greece), Janja Filipi (University of Zadar, Croatia), Francesco Pennacchio and Gennaro Di Prisco (Università degli Studi di Napoli "Federico II", Italy), Raffaele Dall'Olio (Università degli Studi di Udine, Italy), Valentina Cebotari (Academy of Sciences of Moldova, Maldova), Borche Pavlov (MacBee, North Macedonia), Bjørn Dahle and Elisabeth Furuseth Hansen (Norwegian Beekeepers Association/Norwegian School of Veterinary Science, Norway), Eliza Cauia (Beekeeping Research and Development Institute, Romania), Nebosjsa Nedic and Kazimir Matovic (University of Belgrade, Serbia), Ewan Campbell and Alan Bowman (University of Aberdeen, United Kingdom), and finally Maureen Wakefield, J Wilford and Karyn Burke (Food and Environment Research Agency, United Kingdom). This study was supported by the European commission through the 7th framework collaborative project SMARTBEES (Sustainable Management of Resilient Bee Populations, grant agreement KBBE-7-61393). Finally, the authors wish to thank Delphine Libby-Claybrough, professional translator and native English speaker, for her editorial assistance in revising the manuscript.
**Author contributions**

Conception and design of the study: E.D., F.Si., R.T. Coordination of sample collection: F.Si., M.M. Sample treatment and PCR analysis: F.Sc., M.M., F.Si. Sequencing and data interpretation: F.Si., P.L., F.Sc., Y.B., E.D., A.B. Writing original draft: F.Si., E.D., A.B. All the authors reviewed and approved the final version of the manuscript.

**Declarations**

**Competing interests**

The authors declare no competing interests.

**Additional information**

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1038/s41598-025-86937-5.

**Correspondence** and requests for materials should be addressed to E.D.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2025