The honeybee (Apis mellifera) developmental state shapes the genetic composition of the deformed wing virus-A quasispecies during serial transmission

The main biological threat to the western honeybee (Apis mellifera) is the parasitic mite Varroa destructor, largely because it vectors lethal epidemics of honeybee viruses that, in the absence of this mite, are relatively innocuous. The severe pathology is a direct consequence of excessive virus titres caused by this novel transmission route. However, little is known about how the virus adapts genetically during transmission and whether this influences the pathology. Here, we show that upon injection into honeybee pupae, the deformed wing virus type-A (DWV-A) quasispecies undergoes a rapid, extensive expansion of its sequence space, followed by strong negative selection towards a uniform, common shape by the time the pupae have completed their development, with no difference between symptomatic and asymptomatic adults in either DWV titre or genetic composition. This suggests that the physiological and molecular environment during pupal development has a strong, conservative influence on shaping the DWV-A quasispecies in emerging adults. There was furthermore no evidence of any progressive adaptation of the DWV-A quasispecies to serial intra-abdominal injection, simulating mite transmission, despite the generation of ample variation immediately following each transmission, suggesting that the virus either had already adapted to transmission by injection, or was unaffected by it.

The requirement of viruses for a living host in order to replicate axiomatically sets limits to the amount of damage a virus can do to a host at either individual or population level. Virulence is the quantitative expression of such damage [1][2][3] : essentially the compromise between the host's ability to control the virus or its damage versus the ability of the virus to circumvent this control 1 . It is a highly plastic and adaptive trait, particularly through the virus, subject to selection at multiple levels of host organization (molecular, cell, tissue, individual, population, species) with different possible outcomes at the different organizational levels, including levels where the interests of virus and host may be aligned, for example as a weapon in ecological competition with more susceptible hosts 4 . Consequently, virulence evolution is equally varied, resulting in all types of host-virus relationships ranging from parasitic host exploitation by the virus through mutualism to enslavement of the virus in function of the host's interests, including genetic integration of viral and host genomes 5 . The main mechanism for virulence evolution is through the virus genome and in particular its quasispecies, i.e. the population of virus genome variants generated by the high rates of mutation and recombination associated with RNA replication 6 . Through the quasispecies, a virus can leverage its small physical genome into a highly variable sequence space (a theoretical representation of all possible variants of a sequence 6 ) with which to respond dynamically to adaptive challenges, through a combination of natural selection and complementation 7,8 . Consequently, it is the quasispecies as a

Results
Serial DWV-A transmission produces distinct waves of attenuated symptoms. The two surviving Leksand colonies 36 were placed in quarantine in autumn 2009 and moved to Uppsala in summer 2011 for experimentation. A single DWV symptomatic adult from one of these colonies was used as the original source inoculum. Mite-mediated transmission was simulated by microinjection of controlled amounts of virus into pink-eyed honeybee pupae 40 , which were either left to complete development to adults or sacrificed after 24 hours to provide the inoculum for the next transmission ( Supplementary Fig. S1). While DWV-B replicates inside varroa 22,25 , DWV-A does not 41 , such that the microinjection accurately mimics the mechanical vectoring of DWV-A by varroa 41 . The serial transfers produced mostly DWV-symptomatic adult bees, confirming the direct causal relationship between DWV and its symptoms 25,29,30,34 , but with occasional 'waves' of asymptomatic bees in varying proportions (Fig. 1).
Symptom attenuation in adult bees is unrelated to DWV-A titre, consensus sequence or quasispecies composition. The symptomatic and asymptomatic bees at each transmission had similar DWV titres, which were at least 5 orders of magnitude larger than for the mock-inoculated control bees ( Fig. 1) while sequence analysis showed that the DWV in the mock-inoculated samples was identical to the DWV of the serial inocula ( Supplementary Fig. S2), and thus most likely derived from experimental contamination, rather than from pre-existing local DWV in the pupae. These observations indicate that the elevated DWV titres at each transmission were due to the virus present in the serial inocula. The full DWV genome sequences of the first (T0) and final (T13) 24-hour transmission inocula were 100% identical, indicating that there had been no long-term genetic adaptation by DWV to simulated serial varroa transmission. These DWV sequences were also quite distinct from the original Leksand DWV sequences ( Supplementary Fig. S2), casting a degree of uncertainty over the provenance of the original T0 inoculum. The absence of reliable biogeographic clustering throughout the tree (and in DWV phylogenies in general 28,42,43 ) makes it furthermore impossible to either confirm or deny a possible Leksand origin for the T0 inoculum, and thus its primordial status. The question of whether a naïve, primordial DWV quasispecies would adapt genetically to varroa-mediated transmission therefore remains unanswered. the adult DWV-A quasispecies composition is shaped by strong, conservative selection from initially diverse post-inoculation quasispecies. In the absence of consensus DWV sequence change, we decided to investigate possible changes in the genetic composition of the DWV quasispecies through the first wave of asymptomatic bees, in both the 24-hour transmission inocula and the resulting (a)symptomatic adult bees. The DWV quasispecies in each sample was represented by a matrix of Single Nucleotide Polymorphism (SNP) frequencies along the entire DWV genome, which allows all the major types of variation found in viral quasispecies (point mutations, deletions, major strains and recombinants [6][7][8] ) to be represented. However, the DWV quasispecies analysed in these experiments was grouped around only one of the major DWV strains (DWV-A 44 ), and therefore involves a narrower range of genetic variability and more refined resolution of effects than quasispecies consisting of multiple major strains and their recombinants [20][21][22]31,32,35 . Although the DWV consensus sequence remained constant and near-identical to the T0/T13 consensus sequence throughout this first wave of transmissions and adult progeny ( Supplementary Fig. 2), radical and consistent patterns of change were observed in the DWV quasispecies underlying these consensus sequences, as represented by changes in the frequency distributions of the four nucleotides (or a deletion) along the DWV genome, relative to the original inoculum (Fig. 2). The principal distinction is between the relatively diverse range of quasispecies for the 24-hour serial pupal inocula (blue) and the more uniform range for the (a)symptomatic adult progeny of each serial transfer (green). Particularly striking is the strong uniformity of the adult bee quasispecies, both between symptomatic (light green) and asymptomatic (dark green) bees of each serial transfer and, more interestingly, also between different serial transfers, despite a diverse range of 24-hour quasispecies (blue) serving as the inocula for these serial transfers. The overall nucleotide diversity was similar for the DWV quasispecies of the serial inocula and the adult bees, but distributed differently within the genome, with consistent evidence from multiple indicators ( Table 1) that the adult quasispecies are shaped mostly by negative selection, while there is little evidence of any selection (positive, negative or neutral) within inocula quasispecies. This observation does not preclude the possibility of positive, directional selection during quasispecies development; rather that the evidence for this can be quickly obscured by subsequent quasispecies changes.
DWV-A quasispecies composition is linked to developmental state rather than inoculation history. The most parsimonious explanation is that injection caused a rapid burst of replication that generated extensive variation, seen in the 24-hour inocula, followed by largely negative selection as replication slowed down, resulting in a convergence of quasispecies shapes 45 , as seen in the adult bees. The surprising feature is the strength and consistency of this convergence through multiple serial transmissions, considering the diversity of shapes and variation generated by the 24-hour inocula. These observations are confirmed by more formal analyses (Fig. 3). Despite the much longer evolutionary distance between different adult bees than between any bee and its 24-hour inoculum, the adult bee quasispecies cluster together, with little difference between symptomatic and asymptomatic bees, and well separated from the cluster of 24-hour inocula (Fig. 3A,B). Of special interest here is that the original T0 inoculum, classified in the DAPC as an inoculum but derived from an adult bee, locates both in the plot (Fig. 3A) and by genetic assignment (Fig. 3B) very clearly with the adult quasispecies rather than with its pre-assigned group, the inocula. This is the first indication that the host's status may be highly influential in determining the DWV quasispecies composition 46 . A phylogenetic analysis of the composition of the different quasispecies shows that there is no correspondence between the inferred phylogenetic relationships, based on shared genetic signatures, and their known parentage through serial transmission, particularly for the adult bee quasispecies (Fig. 3C). This attests to the extreme plasticity of the quasispecies and its capacity for rapid and extensive change in composition, thereby obscuring the stable phylogenetic signal required for successful reconstruction of the history of descent. Even though the potential for change exists at every transmission event, shown by the extensive initial changes to the DWV quasispecies 24 hours after injection, the end result in the adult host is a conserved standard shape. A similar conservation of quasispecies shape is also observed through the multiple serial transmissions, mimicking sequential transmission within a colony, since the adult DWV quasispecies after the seventh serial transmission, particularly 7AS, 8 S, 8AS, 9 S and 9AS, were very similar to the original T0 inoculum (Figs. 2, 3C). Both the adult quasispecies and those of the 24-hour inocula show a slight phylogenetic clustering by their serial transmission number, with robust assignment of the individual taxa (Fig. 3C), suggesting that at least some of the genetic signatures in the inocula were transmitted through several consecutive transfers  Graphical representation of DWV quasispecies composition change during transmission. Bubblegraph representation of the difference in frequency (f) of the consensus nucleotide (n) at each nucleotide position along the 10 kb DWV genome (nt) between the experimental sample (f(n i )) and the original T0 inoculum (f(n T0 )), shown for each of the successive 24-hour transmission inocula (blue; T1-T8) and their corresponding 7-day symptomatic (light green) and asymptomatic (dark green) adult transmission progeny. The size of each bubble represents the extent of the change in the sub-consensus nucleotide frequency along the DWV genome, i.e. those quasispecies with the smallest bubbles and the least amplitude are most similar to the original T0 inoculum. to subsequent inocula, while different genetic signatures were retained through successive transfers in the adult quasispecies that develop from these inocula (Fig. 3D, Supplementary Fig. S3). However, the fleeting volatility of such signatures, both between transfers and during the transformation of the quasispecies towards the shapes encountered in adult bees, makes it difficult to capture them by clustering analyses, unless they are persistent and extensive.

Discussion
Despite many decades of co-existence between honeybee, V. destructor and DWV, in millions of colonies and billions of honeybees world-wide, the pathology and virulence of DWV remain stubbornly linked to direct vectored transmission by V. destructor in developing pupae 14,[22][23][24][25]34 , with DWV pathology, virulence and titres all disappearing from honeybee colonies when varroa is removed 23,29 . This suggests that DWV genetic adaptation to transmission by V. destructor is either non-existent, and DWV virulence largely determined by the transmission route and viral titres 25 , or ephemeral; present only during varroa-mediated transmission with the DWV quasispecies rapidly returning to a default genetic shape and composition when varroa transmission is removed. The evidence accumulating thus far strongly favours the second mechanism, i.e. the ability of the DWV quasispecies to adapt very rapidly to different selective challenges through compositional changes to the genetic structure of its quasispecies. Most of this evidence comes from studies on the relative virulence of DWV-A and DWV-B, the two major strains of DWV 14,44 . Both strains are capable of producing symptoms and disease upon injection into pupae, either by varroa or experimental inoculation 25,34,47 , but DWV-B is considerably more virulent than DWV-A at both individual bee 22 and colony level 21 and is more closely associated with varroa 22,48 and winter mortality 21,49 . Recombinants between the two strains are possibly even more virulent than either parent 20,31,32, and have made it possible to map this virulence to the RNA-dependent RNA polymerase region of DWV-B 22 . However, virulence also comes at a cost, both individual 22,50-53 and social [54][55][56] , which may explain why the natural incidence of DWV-A/DWV-B recombinant viruses is low 57 , despite having many selective advantages over their parents in laboratory studies 20,22 , and why natural recombinants often contain the "avirulent" DWV-A polymerase region 20,31,32,47,58,59 . For example, both DWV-A and DWV-B are transmitted by varroa, but only DWV-B also replicates in varroa brain tissues 41,48,60 . While this enhances transmission 20 , it also impairs varroa neurological function 61 , which is essential for host finding and reproduction, thus neutralizing some of DWV-B's selective advantage over DWV-A. A similar, but inverse, cost-benefit relation exists in honeybees. DWV-B is more virulent than DWV-A in honeybee pupae 22 , and thus more likely to be detected and removed by honeybee hygienic behaviour 50,52 and, at a higher level, by elevated colony winter mortality 21,49 . However, DWV-B also impairs to a greater extent than DWV-A the cognitive behaviour of adult bees 22 , an essential prerequisite for effective hygienic behaviour 51 , increasing the chance of DWV-infected pupae escaping hygienic surveillance 50 . In both varroa and honeybees therefore, the cost of DWV-B's higher virulence is greater neurological impairment of its host, which can be either beneficial or detrimental to the virus depending on host, developmental stage and colony context. This is a natural foundation for the appearance of multiple strains adapted to different contexts. Since there is little difference between DWV-A and DWV-B in pathogenicity 47,49 , titres 20,33,47,49 or adult mortality 22 , the relative persistence of DWV-A, DWV-B and/or their recombinants in bee colonies may well be driven largely by pupal mortality and the effectiveness of social hygiene 50 , where the strains differ most.
Like others before us, we also could not distinguish DWV quasispecies from symptomatic and asymptomatic adult bees 24 and found the DWV quasispecies dynamics to be dominated by rapid selective processes 20,22,33 . Results elsewhere suggest a limited role for the cooperative structure of viral quasispecies 8,10,62 in DWV quasispecies dynamics 33,60 . No long-term genetic adaptation to transmission by injection was observed in our experiments, suggesting that our DWV-A quasispecies was suitably adapted to mechanical vectoring 41 .
We identified the host developmental state as an additional, distinct and powerful modulator of DWV quasispecies composition, after serially transmitting DWV-A into successive generations of white-eyed pupae, and following the fate of the DWV-A quasispecies through subsequent development into symptomatic or asymptomatic adult bees. The waves of asymptomatic adult bees in the serial transmission progeny suggest that the appearance of these attenuated symptoms has a systemic origin. This could have been host or virus-derived factors   www.nature.com/scientificreports www.nature.com/scientificreports/ by the difference between the pupal and adult developmental states. Such strong host control over the genetic character of viral quasispecies has been demonstrated repeatedly in other experimental systems, usually through predictable and reversible changes associated with transfer between different hosts 46,69 or host tissues [70][71][72] . Despite minor indications that some of the DWV-A quasispecies genetic signatures are transferred serially between inocula, there is no concrete evidence of any progressive directional change in the 24-hour inocula quasispecies composition that would indicate an adaptation to transmission by injection, either from the indicators of selection (Table 1) or from the genetic composition of the quasispecies (Fig. 3D). Both theoretical 9,62 and experimental studies 22,45,46,73,74 have shown that once produced, favourable genetic features or shapes can become fixed very rapidly in a viral quasispecies, often within a few serial passages. Such rapid genetic adaptation has also been inferred from analysing the fate of different DWV master variants and their recombinants 20,21,34 , with different replication characteristics 21,33 , in different transmission environments [20][21][22]31 or hosts 22,75,76 . The current experiments were therefore of sufficient duration and scope to have identified genetic adaptations specific to the transmission route, if these were either available in the original DWV-A quasispecies or produced during the experiment. These experiments also provide a possible mechanistic explanation for those instances where adaptation to transmission route 17 or host species 75 did not result in divergence of the DWV consensus sequence if, as was the case here, the adaptability of the virus involved primarily the re-distribution of its internal genetic variability, rather than fixed changes to the consensus sequence or competition between major strains and recombinants [20][21][22]33 . The results also imply that any features of the DWV quasispecies present during early wing development that may be relevant for the DWV-symptomatic changes in wing morphology will have largely disappeared from the quasispecies by the time the symptoms become visible in the adult bees, helping explain why no genetic markers for DWV symptoms have yet been identified on the DWV genome 24,34 . Higher resolution mapping of the transformation of the DWV quasispecies during pupal development, tissue tropism or host-range studies would provide greater clarity of the extent of host control over the DWV quasispecies, while the mechanisms of such host-specific control would require more comprehensive analyses of the RNA samples produced from such studies. The rapid adaptability and resilience of the DWV quasispecies also helps explain why DWV functionally reverts to its low-titre, benign status when mite-mediated transmission is removed 23 , if this represents the preferred, ancestral relationship between virus and host. Such consistent reversal to low virulence, or to a common quasispecies shape (as shown here), hints at a more complex host-pathogen relationship than mere pathological exploitation and perhaps other, more mutualistic 77 roles for DWV to bee life at some level of organization (cell 51 , organ, bee, colony, population) or context (health 14,16,17 , sociality 50 , ecology 4,78 ), as one of the predicted possible outcomes of host-pathogen co-evolution 8,27,79 . Since these viruses are part of a shared pathosphere involving other bees, wasps and bumblebees 4,57,78 , the question also has relevance outside the beekeeping community 4,75,78 .

Source of bees and virus.
The two surviving colonies in the Leksand study 36 were transferred in September 2009 from Leksand to an isolated forest area near Forsmark (60.393735, 18.208015), at that time also still beyond the varroa front. When Varroa destructor was detected in Forsmark in 2009, the colonies were relocated in March 2010 to a quarantine apiary in Grötlingbo Nature reserve on Gotland (57.121910, 18.419933), 4.8 km from the nearest registered bee colonies, since they could no longer be moved back into varroa-free areas on mainland Sweden or its archipelagos. One year later, in June 2011, the colonies were brought back to Uppsala for experimentation. The DWV starting inoculum (T0), prepared as described below, was obtained from a single symptomatic adult honeybee from one of these colonies. By the end of 2011 both colonies had acquired varroa. A Varroa destructor-free and DWV-free honeybee colony from Utö (58.915134, 18.211738), an island in the Swedish archipelago that has never been infested with V. destructor, was used to obtain the experimental worker pupae for serial transmission. Throughout the experiments, control worker pupae from this colony tested negative for DWV by RT-qPCR.
Serial DWV transmission. The serial transmission experiments were conducted during the summer of 2011, according to the experimental design shown in Supplementary Fig. S1. Each serial inoculum was prepared by macerating the worker pupa in 500 µL PBS buffer, homogenizing with 100 µL chloroform and strong vortexing, and clarifying the extract by centrifugation at 13000 rpm for 10 minutes, retaining the supernatant as inoculum and for RNA extraction 40 . At each transmission passage, around 15 pink-eyed, white-bodied pupae were each microinjected laterally in the abdomen, in between the 2 nd and 3 rd integuments, with 2 µL DWV inoculum obtained from the previous passage, using a 50 µl microsyringe (Hamilton Microliter ™ Syringes, Nevada, USA) and very thin, 30 gauge disposable needles. The needles were replaced between individual pupae. The pupae were incubated on Whatman filter paper in disposable petri dishes in an incubator at 35 °C and 80% relative humidity 40 . After a 24 hour incubation, a single random pupa was selected as the DWV inoculum for the next serial passage, while the remaining pupae were incubated a further 6 days to complete their development into adult bees. All collected bees were stored at −80 °C until further use. This procedure was repeated for 14 serial transmissions per trial, and for two independent trials. For technical controls, ten pink-eyed pupae were microinjected with just 2 μL PBS at each transmission. From the serial transmission assays detailed above, all the serial 24-hour inocula pupae and at least one of the corresponding adult transmission progeny were selected for DWV quasispecies analysis through IonTorrent sequencing. www.nature.com/scientificreports www.nature.com/scientificreports/ at −80 °C. RNA was extracted from these 24-hour and 7-day sample homogenates using NucleoSpin ® RNA Kit (Macherey-Nagel, Düren, Germany) and was eluted into 50 μL nuclease-free water. The approximate RNA concentration was determined by NanoDrop. Reverse transcription was performed in 20 µL reactions containing approximately 2.5 µg RNA, 10 mM dNTPs, 100 µM random hexamer primers and 400 units SuperScript ™ III RT (Invitrogen ™ , Life Technologies Corporation, New York, USA). The resulting cDNA was purified using the High Pure PCR purification Kit (Roche Life Sciences, Penzberg, Germany), eluted into 100 μL elution buffer (10 mM Tris-HCl, pH 8.5) and stored at −80 °C until further use.

RT-qPCR quantification of honeybee viruses.
The levels of several common honeybee viruses (Acute bee paralysis virus -ABPV; Israeli acute paralysis virus -IAPV; Kashmir bee virus -KBV; Black queen cell virus -BQCV; Chronic bee paralysis virus -CBPV; Deformed wing virus strain A (DWV-A) and strain B (DWV-B); Sacbrood virus -SBV; Slow bee paralysis virus -SBPV; Bee Macula-like virus -BeeMLV) were quantified by qPCR of the cDNA using previously described protocols and primers 23 , as was the Apis mellifera β-actin mRNA as an internal reference gene for quality control and data Normalization 40,80 (Supplementary Table S1). The reactions were run in triplicate, using KAPA SYBR ® FAST Universal qPCR kit (KAPA Biosystems), 0.3 μL cDNA and 0.2 μM each of forward and reverse primer in a 12 μL reaction volumes. The qPCR was processed in an ECO ™ Real-Time PCR machine (Illumina, SD, USA) with the following qPCR cycling profile: 3 min at 95 °C followed by 40 cycles of [3 sec at 95 °C -30 sec at 57 °C -data collection]. The amplification was followed by a melting curve analysis by reading the fluorescence at 0.5 °C increments from 55 °C to 95 °C, to verify the specificity of the qPCR products. Apart from DWV, only low levels of BQCV and SBV were detected in these samples, which is consistent with current knowledge of the viruses found naturally in Swedish honeybees 23,81 . DWV whole genome amplification. The whole DWV genome was amplified from one symptomatic adult bee and (where possible) one asymptomatic adult bee from each of the first 9 serial transmissions, as well from the corresponding 24-hour serial transmission inocula, each of which was also derived from a single bee. The DWV genome was amplified in four overlapping sections covered by primer-pairs F27/B1806, F1725/B3095, F3018/ B7295 and F7243/B10101 (Supplementary Table S1). In a few instances fragment F3018/B7295 was replaced by three smaller fragments covered by primer-pairs F3018/B4329, F4220/B5625 and F5625/B7295; and fragment F7243/B10101 by smaller fragments F7243/B8794 and F78688/B10101 (Supplementary Table S1). These primer pairs were designed conservatively to amplify all known DWV variants 44  Sequencing. First, the DWV whole-genome PCR fragments of the T0 and T13 24-hour inocula ( Fig. 1;   Supplementary Fig. S1) were sequenced using Sanger technology and primer-walking, to identify any consensus changes over the entire adaptation experiment. The 694 bp DWV-Lp gene fragment of the Control sample was also sequenced with Sanger technology. When the T0 and T13 DWV full genome sequences were found to be 100% identical by Sanger sequencing, a NGS strategy was designed to assess the variability underlying this consensus sequence uniformity. NGS sequencing was performed using the Ion PGM system (ThermoFisher) at the National Genomics Infrastructure in Uppsala. 0.1-1 µg DNA from each sample was fragmented using the S2 system (Covaris). The fragment ends were repaired and ligated to unique bar-coded adaptors using the Applied Biosystems Library builder (ThermoFisher). The adaptor-linked fragments were amplified using the Ion Xpress ™ Plus gDNA Fragment Library Preparation protocol and selected for a 470 bp target size range (Blue Pippin TM , Sage Science). Library size and concentration were assessed by a Bioanalyzer High Sensitivity Chip (Agilent Technologies) and the Fragment Analyzer system (Advanced Analytical). The DNA samples were pooled together in sets of nine, followed by template preparation using the Ion PGM ™ Template OT2 400 Kit on the Ion OneTouch ™ 2 system (ThermoFisher). Samples were then sequenced on the Ion PGM ™ System with Ion PGM ™ Sequencing 400 Kit on Ion 316v2 chips (ThermoFisher) aimed at a 400 bp read length. Data analysis. The IonTorrent sequence reads were first filtered by bar-code, to identify the sample origin, and then mapped to reference sequences AY292384 (DWV-A), AY251269 (DWV-B) and CEND01000001 (DWV-C) using Tmap included in TorrentSuite 3.6.2 with LifeTechnologies recommended parameters. Variability data and consensus sequences were created using mpileup from SAMtools (v.0.1.8) 82 and an in house script. This resulted in a numerical file for each sample containing the frequencies of each nucleotide or a deletion at each position of the DWV genome between nucleotides 10 and 10,120 of the mapping reference sequence for all samples, as well as the sequencing coverage across the DWV genome (average/max/min per sample = 8952×/14292×/2772×), thus capturing the full genetic complexity of the DWV quasispecies 83 . These files form the basis for the data analyses. Statistical analyses were conducted in R unless stated otherwise.
Diversity and selection analyses. The DWV consensus sequences of the 23 NGS sequenced samples were, with the exception of 3 individual SNPs, identical among each other (Supplementary Fig. S2). Their combined consensus sequence was also identical to the T0 and T13 Sanger sequences, and this combined consensus Scientific RepoRtS | (2020) 10:5956 | https://doi.org/10.1038/s41598-020-62673-w www.nature.com/scientificreports www.nature.com/scientificreports/ sequence has been submitted to GenBank under accession number MN746311. The genetic analyses therefore concern primarily the distribution of the genetic variation underlying these consensus sequences, and how they are related to each other. A sub-consensus sequence database was created for each individual sample by calling the second most frequent nucleotide at all polymorphic nucleotide sites after filtering out indels (insertions/ deletions) and those nucleotides with frequencies lower than 0.2% or higher than 50%. The consensus and sub-consensus sequences thus obtained were aligned using MEGA 84 and saved in NEXUS format. Molecular evolutionary analyses were conducted using DnaSP (v5) 85 and MEGA. Genetic variability from the coding region (polyprotein) within inocula and 7-day sample populations was estimated as nucleotide diversity (π) 86 , consisting of the average number of nucleotide differences per site between all pairs of analyzed sequences, using DnaSP. Similarly, Tajima's D 87 and Fu and Li's D and F 88 were used to test for neutrality and to infer positive (beneficial) or negative (purifying) selection as implemented by DnaSP.
The Codon-based Z test of Selection was performed calculating the number of synonymous substitutions per synonymous site (d S ) and the number of non-synonymous substitutions per non-synonymous site (d N ) between inocula and 7-day sample populations using the Nei and Gojobori's method 89 with MEGA in order to reject the null hypothesis of strict neutrality (H o : d S = d N ) and infer, as well, the balance between neutral mutations, negative selection and positive mutations acting on the DWV quasispecies. clustering analyses. The clustering analysis of the genetic variation between the sequenced samples was performed using Discriminant Analysis of Principal Components (DAPC), a model-free method implemented in the adegenet package for R 90,91 . This method relies on reducing a large SNP dataset to a number of uncorrelated variables by Principal Component Analysis (PCA), which are subsequently used by discriminant analysis (DA) to assign the samples to the different genetic clusters by maximizing between-group variance and minimizing within-group variance. The input for the DAPC consisted of frequencies of polymorphic SNPs (columns) for each of the 23 virus quasispecies populations (rows). The data for sample 8AS was based on a reduced genome, since the sequence coverage for the last 1,302 nucleotides at the 3′ end of the genome was below the ≥100x cut-off value for nucleotide frequency calculations (see Fig. 2). The criteria for SNP inclusion was a MAF > 0.02, which is well above the technical error rates for IonTorrent sequencing 92 . Two SNP datasets were constructed. The first dataset included the 23 DWV quasispecies populations, for which there was a total of 1,212 polymorphic SNPs (Fig. 3B). The second dataset included the eight 24-hour transmission inocula, for which there was a total of 773 polymorphic SNPs (Fig. 3D). The sequential K-means clustering algorithm was run in the second dataset for K = 1 to 8 using 1 × 10 6 iterations per run. The different clustering solutions were compared using the Bayesian Information Criterion (BIC) to identify the optimal number of genetic clusters that describe the data. The optimal K is associated with the lowest BIC value 91 . phylogenetic analysis. The phylogenetic analyses of the SNP frequency data were conducted with various programmes implemented by PHYLYP (v.3.65c) 93 . Nei's genetic distance 86 between the 23 DWV quasispecies samples was calculated from the SNP frequencies of the 1,212 characters in the DPC data matrix using the programme 'Gendist' . This distance matrix was used to construct a phylogenetic topology using the Neighbour-Joining algorithm in the programme 'Neighbour' and the consensus tree was built using the programme 'Consense' . Statistical support for the branches in the consensus topology was provided by 1,000 replicate bootstraps of the data, using the programme 'Seqboot' . The phylogenetic analysis describing the relationship between the DWV samples from the current experiment and other biogeographic isolates from Sweden and Gotland 14,36,93 relative to a selection of DWV-A 44,94 isolates representing Europe and the Middle East 14,20,32,43,58,94,95 , North and South America 33,94,96 and East Asia 14,97-100 was inferred using the Maximum Likelihood method and the Tamura-Nei model of evolution 101 , as implemented by MEGA-X 85 , using the DWV-B 44,48 and DWV-C 44 reference strains as outgroup sequences. The statistical confidence of each of the nodes was determined by bootstrap analysis involving 500 replicates. Initial tree(s) for the heuristic search were obtained automatically by applying Neighbor-Join and BioNJ algorithms to a matrix of pairwise distances estimated using the Maximum Composite Likelihood (MCL) approach, and then selecting the topology with superior log likelihood value. All positions containing gaps and missing data were eliminated resulting in a total of 694 positions in the final dataset.

Data availability
The datasets generated and/or analysed during the current study are available from the corresponding author on reasonable request.