rna viruses have rapid replication and high error rates , leading to immense diversity within each viral replication cycle ( domingo and holland , 1997 ) . as a result , many rna viruses are highly genetically heterogeneous and exist within infected population structures known as quasispecies . it has been proposed that this gives these viral pathogens an increased ability to shift to a new environmental niche , such as a new host , as a suitable mutant is more likely to already exist if this opportunity arises . however , if the diversity within a quasispecies becomes too high , deleterious mutations can accumulate leading to loss of overall fitness ( clarke et al . , 1993 ) . the adaptability and host range of a virus are a function of the level of diversity found within a quasispecies ( termed the quasispecies swarm size ) ( schneider and roossinck , 2001 ) . as well as a randomly generated swarm of mutants around one variant , quasispecies can also exist as a number of master variants , each with their own swarm of random mutations ( palacios et al . , 2008 ; lauring and andino , 2010 ) . determining the extent of genetic heterogeneity in virus populations thus has important implications for predicting and preventing emerging viral diseases . the global decline in honey bee populations over the past few decades has been attributed to the ectoparasitic mite varroa destructor and its affiliation with deformed wing virus ( dwv ) type a variant ( dainat et al . , 2012 ; martin et al . , 2012 the role of honey bees as pollinators is vital to the environment and economy , as bees are a key pollinator species for agriculture . the economic value of honey bees is estimated to be more than 225 billion us dollars worldwide ( gallai et al . , 2009 ) . dwv , a member of the single - stranded positive - sense rna genus iflavirus ( lanzi et al , 2006 ) , exists as a group of closely related viruses , often considered as variants of the same species complex ( ryabov et al . , 2014 ) . the ictv ( the international committee on taxonomy of viruses ) database categorises dwv type a as two variants , dwv ( lanzi et al , 2006 ) and kakugo virus ( kv ) ( fujiyuki et al . , 2004 ) . a second closely related virus varroa destructor virus-1 ( vdv-1 ) is also part of the genus iflavirus , which is now designated dwv type b ( martin et al . , 2012 ; mordecai et al . , 2015 ) . dwv type b was designated a separate species based on a nucleotide identity to dwv type a of 84% that , according to demarcation criteria , is sufficient dissimilarity to warrant the creation of a new master variant ( fauquet et al . , 2005 ) . vdv-1 was originally isolated from v. destructor but has since been reported to replicate in honey bees ( ongus et al . , 2004 ; zioni et al . , 2011 ) where it has been shown to cause wing deformities in bees ( zioni et al . , 2011 ) . dwv type b was also recently found to dominate the dwv population in honey bees from an isolated apiary in swindon , uk ( mordecai et al . , 2015 ) ; to the apparent exclusion of any other master variants . in honey bee populations that have never been exposed to varroa mites , dwv exists in a very large variant swarm with numerous master variants ( martin et al . , 2012 ) . however , transmission of dwv by varroa reduces variant diversity to one master variant ( martin et al . , 2012 ) . therefore , dwv exists as an endlessly mutating swarm of variants with these master variants constituting part of this quasispecies and sharing a recent common ancestor ( baker and schroeder , 2008 ; martin et al . , 2012 ) . the dominance of one master variant over another will lead to ultimately different life histories for the colony , that is , death if dwv type a dominates ( martin et al . , 2012 ) or health if dwv type b dominates ( mordecai et al . , 2015 ) . by existing as a diverse swarm of variants , certain biological traits may allow a virus to infect one cell type over another , known as cell tropism ( koyanagi et al . , 1987 ) . whether or not a virus is able to infect a susceptible cell depends firstly on recognition of a cellular receptor on the cell surface and secondly on intracellular host factors that dictate whether the host cell is permissive to virus replication . therefore , amino acid substitutions caused by nucleotide mutations in the structural or nonstructural region of the virus genome can affect both the host range and cell tropism of a virus . when categorising viruses based on a phylogenetic relationship it is important to note that a single amino acid change can have a substantial effect on the phenotypic traits of a virus . therefore , when a virus exists as a collection of variants or quasispecies , although the phylogeny and ancestry of the viruses may be similar , the host range , tropism , pathogenicity and epidemiology of the variants may differ greatly ( domingo et al . in addition , recombination between these variants is a source of further variation ( moore et al . , 2011 ) . there are numerous biological implications of quasispecies occupying large amounts of sequence space that challenge the treatment and control of both established and emerging infectious diseases ( gomez et al . , 1999 ) . for example , quasispecies theory has been used to describe how viruses such as hepatitis c virus and hiv are able to escape host immune responses ( pavio and lai , 2003 ; woo and reifman , 2012 ) . in addition , effective vaccines are yet to be developed for both these viruses because of the rapid emergence of resistant mutants under vaccine selection pressure ( gaschen et al . , 2002 ; law et al , 2013 ) . for similar reasons , the highly divergent nature of rna viruses has been implicated in the emergence of antiviral drug - resistant infections in aids ( metzner et al . , 2009 ) , hepatitis c ( halfon and locarnini , 2011 ) , hepatitis b ( nishijima et al . , 2012 ) and influenza ( boivin et al . , 2002 ) . conventional methods to define and analyse variance within rna viruses include rna  fingerprinting ' ( domingo et al . , 1978 ) and reverse transcriptase - pcr amplification using specifically designed primers ( see , for example , highfield et al . , clone libraries and sanger sequencing were used to identify the dwv type a as being associated with varroa infestation and colony collapse ( martin et al . , 2012 ) . although these techniques are valid for identifying known variants , primer - based methods are prone to missing unidentified variants and are biased towards overrepresented sequences ( gomez et al . , 1999 ) . pcr - based methods are less appropriate to determine the extent of variation in a quasispecies that is not normally distributed that is , where multiple variants exist , each with their own spectrum of mutants ( gomez et al . , 1999 ) . in these instances deep sequencing methods such as illumina platforms are more suited to discovering new variants as well as to diversity analysis ( wood et al . , 2014 ) . illumina sequencing is advantageous for samples with large amounts of genetic variation because of its depth of sequencing , although significant raw read analysis is required . reference - based assembly methods can overlook biological variants because of inaccurate read alignments and loss of data ( archer et al . , 2010 ; iqbal et al . , 2012 ; yang et al . , 2012 ) , although former studies using high - throughput next - generation sequencing of dwv have used these methods ( moore et al . , 2011 ; ryabov et al . , a previous study ( highfield et al . , 2009 ) found that despite controlling varroa populations , high dwv loads were associated with overwintering colony losses ( ocl ) . historical losses due to ocl of 10% were normal ; these have now risen to 20% since the establishment of varroa . this suggests that a  new ' non - varroa - transmitted dwv master variant may be circulating in some colonies during the winter causing ocl . to investigate this hypothesis we used a bespoke de novo assembly pipeline ( mordecai et al . , 2015 ) . the vicuna ( broadinstitute.org/scientific-community/science/projects/viral-genomics/vicuna ) de novo assembler was used as it is designed to assemble highly heterogeneous viral populations and is well suited to the computational challenge that the dwv quasispecies present ( yang et al . , 2012 ) . as well as assembling the dwv type a master variant , a third dwv master variant ( type c ) was assembled and confirmed to be distinct through phylogenetic inference . illumina ( san diego , ca , usa ) hi - seq ( 2  100 ) pair - end sequencing was carried out by tgac ( the genome analysis centre ) and at the university of exeter on four samples from highfield et al . ( 2009 ) ; late june , early october and late october from gd1 , and early october from gd2 ( supplementary table s1 ) . the samples were originally collected from devon in the southwest of england ( highfield et al . , 2009 ) . a bioinformatics pipeline was developed to accommodate the large amount of variation found within the dwv species complex . first , the quality of the raw reads was verified using fastqc ( babraham bioinformatics , cambridge , uk ) . samples were then converted from fastq to fasta using the fastq_to_fasta script that is part of the fastx - toolkit ( hannon lab ) ( http://hannonlab.cshl.edu/fastx_toolkit/ ) . to isolate the dwv complex sequence reads from the host and other contaminating sequences , the reads were searched against a custom blast database containing the dwv , vdv-1 and kv genomes , with an e - value of 10e5 . blast was carried out against read 1 of the illumina data . the ncbi - blast - parser perl script ( http://www.bioinformatics-made-simple.com/2012/07/ncbi-blast-parser-extract-query-and.html ) sed ' and  awk ' scripts were used to delete empty lines and the reads that contained  nohits ' . the corresponding blast hits were extracted from the read 2 raw reads using qiime ( caporaso et al . , 2010 ) . the paired reads were balanced using a custom script written in r version 3.2.0 ( r core team , 2015 ) . finally , the balanced dwv family reads were assembled using the vicuna assembler that was developed to generate consensus assemblies from genetically heterogeneous populations , specifically rna viruses . vicuna contigs > 200  bp in length were imported into geneious ( version 7.04 , created by biomatters , auckland , new zealand ) and the  map to reference tool ' was used to align the contigs with the dwv and vdv-1 reference genomes . for several of the samples the vicuna assembly yielded full - length contigs that covered the whole genome , whereas for others consensus scaffolds were created from two or more contigs . a second novel variant was also assembled by vicuna . in order to create a consensus sequence of the novel variant , any contigs containing sequence other than the novel variant were removed . if a contig contained the novel variant as well as sequence belonging to the type a dwv genome because of recombination or in silico recombination , the type a dwv regions were trimmed and deleted . a consensus genome of the novel variant was created from these sorted and trimmed contigs . accordingly , the same steps were carried out on the type a sequence in order to assemble the type a genome consensus from devon . the dwv scaffolds were then aligned with full - length genome sequences from the ncbi ( national center for biotechnology information ) database using the muscle alignment tool ( edgar , 2004 ) within geneious and the full - length genome was obtained . full genome comparisons were visualised in mvista ( frazer et al . , 2004 ) . to quantify the number of reads attributed to each dwv variant using blast ( supplementary table s1 ) , the novel type c variant genome sequence was added to the custom dwv family database and the blast for each sample was carried out a second time . this allowed the type c variants to be categorised correctly by blast rather than attributed to the closets hit to the three other reference genomes . to validate the quantification of reads attributed to each variant , the balanced read 1 and 2 blast hits were aligned competitively to three reference genomes using the  map to reference tool ' within geneious . consensus genomes of type a and type c were translated using the  live annotation ' tool within geneious . because of repetitive regions in the genome , assembly error took place in several sites , leading to an incomplete open reading frame . these assembly errors were corrected by using the vicuna analysis tool to pull out individual reads covering these regions , and the correct sequence was determined . finally , the rdp4 programme ( martin et al . , 2010 ) was used to determine whether any recombination took place between the sequences . assembled genomes are available from the european nucleotide archive under the accession numbers ers657948 ( type a ) and ers657949 ( type c ) . to establish the phylogenetic relationship of dwv subtypes and closely related iflaviruses , we reconstructed the phylogeny of the conserved rdrp amino acid sequences for seven dwv subtypes , spanning all three types that were either sequenced and assembled from the devon hive or available from genbank ( type a : nc_005876.1 , nc_004830.2 ; type b : kc_786222.1 , nc_006494.1 , jq_413340 ) , as well as formica exsecta virus 1 ( nc_023022.1 ) and sacbrood virus ( nc_002066.1 ) . we used a bayesian approach using mrbayes ( v. 3.1.2 ) ( huelsenbeck and ronquist , 2001 ) . we assumed a fixed rate model of protein evolution and reconstructed the phylogeny using a model jumping method . this method allows for different models of amino acid substitution to be used in the markov chain monte carlo ( mcmc ) procedure , with all models contributing to the final result weighted according to their respective posterior probability . we ran two runs of four chains for 4  000  000 mcmc generations , sampling trees every 1000 generations . all trees were drawn using figtree v.1.4.2 ( http://tree.bio.ed.ac.uk/software/figtree/ ) . to estimate the evolutionary rate of dwv and its subtypes , we collated three independent data sets with temporal information for at least one population : ( a ) partial lp - gene ; 10 populations , n=78 , 328  nt ( 13231650  bp ) ( genbank aj489744 , ay292384 , gu109335 , hm162355 , jf346615-jf346620 , jf346624-jf346629 , jf346633-jf346639 , kf164292 , kf164293 , kj437447 , kp734726 , kp734738 , kp734747 , kp734765-kp734770 , kp734774-kp734787 , kp734817-kp734825 , kp734827-kp734846 ) ; ( b ) partial capsid - gene ; 1 population , n=167 , 1215  nt ( 26343848  bp ) ( genbank ay292384 , hq655502-hq655561 , kf314827- kf314932 ) ; ( c ) partial rdrp - gene ; 1 b. terrestris population , n=145 , 508  nt ( 80168522  bp ) ( genbank kp734326kp734470 ) . to test whether these fragments contain a molecular clock signal , we estimated the root - to - tip divergence using path - o - gen v.1.4 ( http://tree.bio.ed.ac.uk/software/pathogen/ ) to estimate how much genetic variation we constructed maximum likelihood phylogenetic trees without the assumption of a molecular clock in phylip v. 3.695 ( evolution.genetics.washington.edu ; felsenstein , 1989 ) and tested whether the regression between root - to - tip distance in these maximum likelihood trees and the age of the samples indicated a clock - like signal . in addition , we tested for a significant temporal signal by randomising the temporal information across each data set 100 times and compared the resulting random evolutionary rates from beast with the real data set . a temporal signal is supported if there is a significant difference between the real data set and the randomised data sets ( ramsden et al . , 2009 ; alizon and fraser , 2013 ) . both analyses support a temporal signal in these data sets ( supplementary figure s3 and supplementary table s4 ) . to determine the appropriate molecular clock models for each data set , we used the path sampling maximum likelihood estimator implemented in beast 1.8 ( baele et al . , 2012 , 2013 ) . as this method is very computationally intensive , we first used a range of simpler tests to limit the number of models to be compared by this method as suggested by alizon and fraser ( 2013 ) , and drummond and bouckaert ( 2014 ) . for the three large temporal data sets , we used jmodeltest v.2.1.1 ( https://github.com/ddarriba/jmodeltest2 ; posada , 2008 ) to compare substitution models based on the bayesian information criterion ( alizon and fraser 2013 ) ; the resulting substitution models were the tamura - nei ( tnr ) model ( lp - fragment ) and the hasegawa , kashino and yano ( hky ) model ( vp3- and rdrp - fragments ) , both with gamma variation . for the dwv - subtype data set , we compared the general time - reversible model and the srd06 model ( shapiro et al . , 2006 ) using path sampling . we partitioned substitution rates between the first and second and third codon positions as , for all data sets , the third codon position had a significantly higher rate . we tested whether a strict clock rate can be excluded by running models with a lognormal relaxed clock ; if the relaxed clock 's coefficient of variation statistic abuts the zero boundary , a strict clock can not be excluded ( gray et al . , 2011 ) . we ran models with exponential population growth for all models ; if the exponential growth rate was significantly higher than zero , a constant population size can be excluded . we then used path sampling to distinguish between clock models ( exponential and lognormal relaxed clock ) and demographic models ( constant population size , exponential population growth or a gaussian markov random field ( gmrf skyride ) ; drummond et al . , 2002 ; minin et al . , based on these analyses , we chose an exponential relaxed clock and exponential growth prior for the lp- and capsid - fragments and a lognormal relaxed clock and constant growth prior for the rdrp - fragment as well as for the dwv - subtype analyses . to generate a genome - wide estimate for the evolutionary rate in dwv , we calculated the mean of the relaxed clock means for the individual fragments ( lp - fragment : 9.097  10 ( 95% highest posterior density : 4.412  101.394  10 ) , capsid - fragment : 1.845  10 ( 1.159  102.569  10 ) , rdrp - fragment : 1.278  10 ( 4.131  102.513  10 ) , resulting in a mean evolutionary rate of 1.346  10 ( 5.41  102.627  10 ) substitutions / site / year . we implemented this evolutionary rate as an uncorrelated lognormal clock prior for the dwv - subtype analysis with a lognormal distribution with a mean of 1.35  10 in real space and a log s.d . we ran models with 2 runs each of 50 million mcmc generations , sampling every 5000 generations with a burn - in of 5 million generations to obtain effective sample sizes > 200 . we examined traces for convergence using tracer v.1.6 ( http://www.tree.bio.ed.ac.uk/software/tracer/ ) and used treeannotator v.1.8 ( beast.bio.ed.ac.uk/downloads ) to produce a maximum clade credibility tree for the dwv - subtype analysis . ( 2003 ) , as implemented in dambe , to confirm that the alignment had not reached substitution saturation . illumina ( san diego , ca , usa ) hi - seq ( 2  100 ) pair - end sequencing was carried out by tgac ( the genome analysis centre ) and at the university of exeter on four samples from highfield et al . ( 2009 ) ; late june , early october and late october from gd1 , and early october from gd2 ( supplementary table s1 ) . the samples were originally collected from devon in the southwest of england ( highfield et al . , 2009 ) . a bioinformatics pipeline was developed to accommodate the large amount of variation found within the dwv species complex . first , the quality of the raw reads was verified using fastqc ( babraham bioinformatics , cambridge , uk ) . samples were then converted from fastq to fasta using the fastq_to_fasta script that is part of the fastx - toolkit ( hannon lab ) ( http://hannonlab.cshl.edu/fastx_toolkit/ ) . to isolate the dwv complex sequence reads from the host and other contaminating sequences , the reads were searched against a custom blast database containing the dwv , vdv-1 and kv genomes , with an e - value of 10e5 . blast was carried out against read 1 of the illumina data . the ncbi - blast - parser perl script ( http://www.bioinformatics-made-simple.com/2012/07/ncbi-blast-parser-extract-query-and.html ) sed ' and  awk ' scripts were used to delete empty lines and the reads that contained  nohits ' . the corresponding blast hits were extracted from the read 2 raw reads using qiime ( caporaso et al . , 2010 ) . the paired reads were balanced using a custom script written in r version 3.2.0 ( r core team , 2015 ) . finally , the balanced dwv family reads were assembled using the vicuna assembler that was developed to generate consensus assemblies from genetically heterogeneous populations , specifically rna viruses . vicuna contigs > 200  bp in length were imported into geneious ( version 7.04 , created by biomatters , auckland , new zealand ) and the  map to reference tool ' was used to align the contigs with the dwv and vdv-1 reference genomes . for several of the samples the vicuna assembly yielded full - length contigs that covered the whole genome , whereas for others consensus scaffolds were created from two or more contigs . a second novel variant was also assembled by vicuna . in order to create a consensus sequence of the novel variant , any contigs containing sequence other than the novel variant were removed . if a contig contained the novel variant as well as sequence belonging to the type a dwv genome because of recombination or in silico recombination , the type a dwv regions were trimmed and deleted . a consensus genome of the novel variant was created from these sorted and trimmed contigs . accordingly , the same steps were carried out on the type a sequence in order to assemble the type a genome consensus from devon . the dwv scaffolds were then aligned with full - length genome sequences from the ncbi ( national center for biotechnology information ) database using the muscle alignment tool ( edgar , 2004 ) within geneious and the full - length genome was obtained . full genome comparisons were visualised in mvista ( frazer et al . , 2004 ) . to quantify the number of reads attributed to each dwv variant using blast ( supplementary table s1 ) , the novel type c variant genome sequence was added to the custom dwv family database and the blast for each sample was carried out a second time . this allowed the type c variants to be categorised correctly by blast rather than attributed to the closets hit to the three other reference genomes . to validate the quantification of reads attributed to each variant , the balanced read 1 and 2 blast hits were aligned competitively to three reference genomes using the  map to reference tool ' within geneious . consensus genomes of type a and type c were translated using the  live annotation ' tool within geneious . because of repetitive regions in the genome , assembly error took place in several sites , leading to an incomplete open reading frame . these assembly errors were corrected by using the vicuna analysis tool to pull out individual reads covering these regions , and the correct sequence was determined . finally , the rdp4 programme ( martin et al . , 2010 ) was used to determine whether any recombination took place between the sequences . assembled genomes are available from the european nucleotide archive under the accession numbers ers657948 ( type a ) and ers657949 ( type c ) . to establish the phylogenetic relationship of dwv subtypes and closely related iflaviruses , we reconstructed the phylogeny of the conserved rdrp amino acid sequences for seven dwv subtypes , spanning all three types that were either sequenced and assembled from the devon hive or available from genbank ( type a : nc_005876.1 , nc_004830.2 ; type b : kc_786222.1 , nc_006494.1 , jq_413340 ) , as well as formica exsecta virus 1 ( nc_023022.1 ) and sacbrood virus ( nc_002066.1 ) . we used a bayesian approach using mrbayes ( v. 3.1.2 ) ( huelsenbeck and ronquist , 2001 ) . we assumed a fixed rate model of protein evolution and reconstructed the phylogeny using a model jumping method . this method allows for different models of amino acid substitution to be used in the markov chain monte carlo ( mcmc ) procedure , with all models contributing to the final result weighted according to their respective posterior probability . we ran two runs of four chains for 4  000  000 mcmc generations , sampling trees every 1000 generations . to estimate the evolutionary rate of dwv and its subtypes , we collated three independent data sets with temporal information for at least one population : ( a ) partial lp - gene ; 10 populations , n=78 , 328  nt ( 13231650  bp ) ( genbank aj489744 , ay292384 , gu109335 , hm162355 , jf346615-jf346620 , jf346624-jf346629 , jf346633-jf346639 , kf164292 , kf164293 , kj437447 , kp734726 , kp734738 , kp734747 , kp734765-kp734770 , kp734774-kp734787 , kp734817-kp734825 , kp734827-kp734846 ) ; ( b ) partial capsid - gene ; 1 population , n=167 , 1215  nt ( 26343848  bp ) ( genbank ay292384 , hq655502-hq655561 , kf314827- kf314932 ) ; ( c ) partial rdrp - gene ; 1 b. terrestris population , n=145 , 508  nt ( 80168522  bp ) ( genbank kp734326kp734470 ) . to test whether these fragments contain a molecular clock signal , we estimated the root - to - tip divergence using path - o - gen v.1.4 ( http://tree.bio.ed.ac.uk/software/pathogen/ ) to estimate how much genetic variation can be explained by the sampling date . we constructed maximum likelihood phylogenetic trees without the assumption of a molecular clock in phylip v. 3.695 ( evolution.genetics.washington.edu ; felsenstein , 1989 ) and tested whether the regression between root - to - tip distance in these maximum likelihood trees and the age of the samples indicated a clock - like signal . in addition , we tested for a significant temporal signal by randomising the temporal information across each data set 100 times and compared the resulting random evolutionary rates from beast with the real data set . a temporal signal is supported if there is a significant difference between the real data set and the randomised data sets ( ramsden et al . , 2009 ; alizon and fraser , 2013 ) . both analyses support a temporal signal in these data sets ( supplementary figure s3 and supplementary table s4 ) . to determine the appropriate molecular clock models for each data set , we used the path sampling maximum likelihood estimator implemented in beast 1.8 ( baele et al . , 2012 , 2013 ) . as this method is very computationally intensive , we first used a range of simpler tests to limit the number of models to be compared by this method as suggested by alizon and fraser ( 2013 ) , and drummond and bouckaert ( 2014 ) . for the three large temporal data sets , we used jmodeltest v.2.1.1 ( https://github.com/ddarriba/jmodeltest2 ; posada , 2008 ) to compare substitution models based on the bayesian information criterion ( alizon and fraser 2013 ) ; the resulting substitution models were the tamura - nei ( tnr ) model ( lp - fragment ) and the hasegawa , kashino and yano ( hky ) model ( vp3- and rdrp - fragments ) , both with gamma variation . for the dwv - subtype data set , we compared the general time - reversible model and the srd06 model ( shapiro et al . , 2006 ) using path sampling . we partitioned substitution rates between the first and second and third codon positions as , for all data sets , the third codon position had a significantly higher rate . we tested whether a strict clock rate can be excluded by running models with a lognormal relaxed clock ; if the relaxed clock 's coefficient of variation statistic abuts the zero boundary , a strict clock can not be excluded ( gray et al . , 2011 ) . we ran models with exponential population growth for all models ; if the exponential growth rate was significantly higher than zero , a constant population size can be excluded . we then used path sampling to distinguish between clock models ( exponential and lognormal relaxed clock ) and demographic models ( constant population size , exponential population growth or a gaussian markov random field ( gmrf skyride ) ; drummond et al . based on these analyses , we chose an exponential relaxed clock and exponential growth prior for the lp- and capsid - fragments and a lognormal relaxed clock and constant growth prior for the rdrp - fragment as well as for the dwv - subtype analyses . to generate a genome - wide estimate for the evolutionary rate in dwv , we calculated the mean of the relaxed clock means for the individual fragments ( lp - fragment : 9.097  10 ( 95% highest posterior density : 4.412  101.394  10 ) , capsid - fragment : 1.845  10 ( 1.159  102.569  10 ) , rdrp - fragment : 1.278  10 ( 4.131  102.513  10 ) , resulting in a mean evolutionary rate of 1.346  10 ( 5.41  102.627  10 ) substitutions / site / year . we implemented this evolutionary rate as an uncorrelated lognormal clock prior for the dwv - subtype analysis with a lognormal distribution with a mean of 1.35  10 in real space and a log s.d . we ran models with 2 runs each of 50 million mcmc generations , sampling every 5000 generations with a burn - in of 5 million generations to obtain effective sample sizes > 200 . we examined traces for convergence using tracer v.1.6 ( http://www.tree.bio.ed.ac.uk/software/tracer/ ) and used treeannotator v.1.8 ( beast.bio.ed.ac.uk/downloads ) to produce a maximum clade credibility tree for the dwv - subtype analysis . ( 2003 ) , as implemented in dambe , to confirm that the alignment had not reached substitution saturation . illumina hi - seq ( 2  100 ) pair - end sequencing of asymptomatic honey bees from devon was carried out on colonies that either survived ( gd2 ) or collapsed ( gd1 ) because of ocl ( highfield et al . , 2009 ) . to assess dwv diversity within the samples the raw reads were searched against a custom dwv family database ( including kv and vdv-1 genomes ) and reads matching to any of the dwv genomes were extracted . the coverage for the number of reads that matched to a custom dwv blast database was estimated using the lander / waterman equation ( that is , the depth of sequencing ) . genome coverage depth ranged from 457 to 165  927  ( average coverage was 86  838  ) ( supplementary table s1 ) . de novo assembly of the illumina reads yielded a complete genome of a type a variant as well as of a novel dwv variant that we named type c ( supplementary table s2 ) . competitive alignment to the now three master variants revealed that out of the ~30 million dwv blast positive reads , ~27 million were assembled to one of three genomes ( supplementary table s3 and supplementary figure s1 ) . around 3 million assembled to type a , ~3.5 thousand to type b and ~24 million to type c. as expected , more reads align to the 3 region of the genome than the 5 an artefact of the reverse transcription 3 bias ( brooks et al . , 1995 ) . although reads aligned to the whole of the type c genome , the depth of coverage was lower at the 5 region ( supplementary figure s1 ) . in addition , coverage of the 5 region of the type a variant was unusually high , indicating that this disproportionate coverage in the 5 regions of type a and c could be the evidence of recombination between the two variants . because of the low read depth and coverage for type b in our samples , a full genome could not be assembled ( supplementary table s3 and supplementary figure s1 ) . nonetheless , closer examination of 3 region where the rna - dependent rna polymerase ( rdrp ) gene is located ( baker and schroeder , 2008 ) revealed that the devon type b variant shared 100% identity to vdv-1 type b genome in this region ( supplementary figure s2 ) . the sequence identity of our newly assembled type a and c variants was compared with other members of the dwv complex ( table 1 and figure 1 ) . both type b and c differ from the type a nucleotide sequence in similar regions of the genome . however , type b and c share only 79% nucleotide identity in the polyprotein encoding region of the genome and 89% identity in the amino acid sequence . phylogenetic analysis showed the relationship of the novel type a and c variants from devon to previously sequenced viruses ( figures 2 and 3 ) . the type a genome clustered with other type a variants , whereas type c formed a distinct and separate branch , thereby confirming the originality and thus new master variant assignment of dwv type c. given that type a and c can recombine with each other , we screened for any genomes with evidence of recombination and excluded them from the analysis in order to comply with the assumptions underlying phylogenetic reconstruction . the type c genome codes for a full - length polyprotein and translation of the genome permitted the amino acid sequence to be compared with more distant members of the genus iflavirus ( figure 2 ) . the relationship of the dwv complex with f. exsecta virus 2 and sacbrood virus was analysed using a bayesian inference of phylogeny of a conserved region of the rdrp amino acid sequence ( figure 2 ) . the dwv sequences were attributed to three main groups ( type a , b and c ) . the nucleotide consensus sequences created by vicuna for type a was found to be 98.2% identical to the type a reference dwv genome . the analysis suggested that the new type c genome is clearly related to type a and b , and forms a distinct clade to other members of the family iflaviridae ( for example , sacbrood virus ) that are all only distantly related to a dicistrovirus outgroup , f. exsecta virus . in addition , a bayesian analysis of the polyprotein encoding region of the sequences was carried out using an mcmc model , permitting a molecular clock model to be run within beast v1.8.1 ( figure 3 ) ( drummond et al . , 2012 ) . divergence times were calculated based on a tip - dated coalescent model , with an evolutionary rate prior based on three independent tip - dated fragments of dwv type a. the samples in the tree span 11 years ( 20002011 ) , and the samples used for estimating the evolutionary rate span 11 , 13 and 22 years for the rdrp- , capsid- and lp- fragments , respectively . the bayesian tree of the nucleotide sequences had a similar structure to the amino acid tree , showing that the type c viral variant is distinct from type a and b. the molecular clock estimation predicts that type c diverged from the other dwv variants 319 years ago ( 571010 95% highest posterior density ) , and type a and b disassociated from each other 181 years ago ( 38497 95% highest posterior density ) . this estimate is unlikely to be biased by substitution saturation , as there is no evidence for saturation in our data set ( xia et al . , 2003 ) . the phylogenetic analysis concludes that the dwv quasispecies is made up of three distinct master variants . however , the number of variants is not exhaustive and it is conceivable that more variants will be added in the future as the use of ngs becomes increasingly prolific . the result further suggests that type c has not recently emerged , but rather is an established dwv variant . moreover , using reverse transcriptase - pcr to amplify a region of the dwv rdrp gene , martin et al . ( 2012 ) attributed a novel variant to the type b swarm of variants , although it is now clear that this novel variant was in fact type c. the presence of type c in hawaii implies that this variant is widespread and not specific to the united kingdom . it also confirms that it has not emerged recently as the hawaiian sample originates from kauai that has yet to be colonised by varroa ( martin et al . , 2012 ) , suggesting type c is part of the wild - type dwv quasispecies . the high levels of heterogeneity within viral rna populations mean that viruses are able to occupy large areas of sequence space and consequently are able to exist in multiple hosts ( domingo and holland , 1997 ) . dwv is a  generalist ' known to infect bumblebees ( genersch et al , 2006 ; frst et al . , 2014 ) , v. destructor ( ongus et al . , 2004 ) and other insects ( reviewed in manley et al . , 2015 ) . further work is required to ascertain whether type c causes acute wing deformities in honey bees or its presence is because of viral  spill over ' from another host . three of the four samples in this study were taken from a hive ( gd1 ) that went on to collapse from ocl . however , as type a and the type a c 5 recombinant were also present , it remains unclear which variant was responsible for ocl . it is worth noting however that between 85% and 98% of the reads in hive gd1 were type c that make up the type c or the a the second hive ( gd2 ) that survived ocl contained roughly equal levels of type a and c reads ( which includes reads that make up the a c recombinant ) ; yet , the coverage indicates that the viral load was much lower in this hive . moreover , the low abundance of dwv type b in these devon colonies appears to confirm the observation made in the colonies from swindon ( mordecai et al . , 2015 ) ; that is , dwv type b is not present at sufficient levels to protect the honey bee from the virulent type a , or possibly c variants ( as well as any recombinants between a and c ) . in addition , varroa mites in the devon colonies were controlled using chemical methods , thereby potentially preventing the mites from transmitting type b into the honey bees ; as observed in the swindon apiary ( mordecai et al . , 2015 ) . a recombinant between dwv type a and b has been reported previously ( moore et al . , 2011 ) , and found to be hypervirulent as it was more efficient at replicating than other variants when co - injected directly into the haemolymph of honey bee larvae ( ryabov et al . as it is presently unclear whether the type c variant or the a  c recombinant lead to disease , further direct manipulation experiments such as those carried out by ryabov et al . ( 2014 could be crucial in determining whether the novel type c genome and any of its recombinants are hypervirulent . the large number of dwv type c and a c recombinant reads in the devon samples ( supplementary table s1 and supplementary figure s1 ) suggests that these variants are able to replicate in honey bees . the capacity of dwv to exist as a swarm of recombining variants and occupy large amounts of sequence space may be one of the factors contributing to its capability to maintain a persistent infection . defining dwv as a quasispecies with discrete master variants is comparable with the current categorisation of several other rna viruses . the genetic variability within hepatitis c virus has been classified into four hierarchical strata : genotypes , sub - genotypes , isolates and finally background variation of the quasispecies ( farci and purcell , 2000 ) . a total of seven hepatitis c virus genotypes exist to date , differing from each other by > 15% over their complete coding region ( smith et al . , 2014 ) . similarly , japanese encephalitis virus is currently classified into five genotypes based on sequence identity that differ in replication efficiency in a range of hosts ( han et al . , 2014 ) . in order to further our understanding of dwv and the impact of the virus on honey bees as well as other hosts , a more progressive classification of the dwv quasispecies would be prudent , that is , the new classification of dwv as a quasispecies with at least three master variants that can recombine with each other . bioinformatics pipelines designed for eukaryotes and prokaryotes are not necessarily suited to studying viral systems . for example , pipelines incorporating reference assemblers could hinder the discovery of novel virus variants . the genome sequences presented here are representative of the different variants found within a sample ; however , the vicuna pipeline is not suited to analysing the further level of diversity around these variants , for which diversity models are better suited , such as those developed by wood et al . ( 2014 ) . as well as being of significance to globally important honey bee health , cross - species virus transmission and emergence of new epidemic diseases such as severe acute respiratory syndrome , ebola and influenza are major threats to public health ( parrish et al . , exploring the extent of viral diversity in rna quasispecies , of which dwv may be a suitable model , may offer insight into the mechanisms by which viruses are able to transmit between different hosts as well as how viruses are able to develop resistance to antiviral therapies ( domingo et al . , 2012 ) . further study of the dwv quasispecies may help to explain how dwv is able to exist as multiple variants in many hosts , and may elucidate mechanisms by which it establishes a persistent infection among several hosts but only proves pathogenic in some .
treatment of emerging rna viruses is hampered by the high mutation and replication rates that enable these viruses to operate as a quasispecies . declining honey bee populations 
 have been attributed to the ectoparasitic mite varroa destructor and its affiliation with deformed wing virus ( dwv ) . in the current study 
 we use next - generation sequencing to investigate the dwv quasispecies in an apiary known to suffer from overwintering colony losses . 
 we show that the dwv species complex is made up of three master variants . 
 our results indicate that a new dwv type c variant is distinct from the previously described types a and b , but together they form a distinct clade compared with other members of the iflaviridae . 
 the molecular clock estimation predicts that type c diverged from the other variants 319 years ago . 
 the discovery of a new master variant of dwv has important implications for the positive identification of the true pathogen within global honey bee populations .