BACKGROUND
in addition to producing honey and other hive products, managed honey bee colonies are vital for the pollination of many crops.
pests and diseases that impact honey bees are therefore of economic concern, particularly in light of the long-distance transport of bees, bee products, and beekeeping equipment that can facilitate their spread.
high rates of unexplained colony loss in recent years  <cit>  have heightened these concerns and provided impetus for surveys of geographic and genetic variation in honey bee pathogens  <cit> .
a number of rna viruses have been recognized as major disease agents of honey bees.
in addition to the recognized morbidity and mortality that can result from overt viral infections, lower viral loads that do not produce visible symptoms might nonetheless impact colony health as well.
for example, viral loads are negatively associated with colony survival  <cit>  and there is correlational evidence of synergism among viruses in colony collapse disorder   <cit> .
interactions have also been identified between viruses and other parasites, particularly the ectoparasitic mite varroa destructor <cit> .
recent genomic studies have identified novel viruses that may also be relevant to colony health, but about which little is yet known  <cit> .
additionally, recombination between co-infecting rna viruses may facilitate adaptation or the emergence of novel viruses  <cit> .
it has therefore become a major goal to better understand the diversity of viral species in honey bee populations and how they may interact during co-infection.
understanding the genetic diversity within a viral species is important as well.
rna viruses in particular have high mutation rates, can evolve rapidly, and may exhibit adaptation to specific hosts and/or tissues.
virulence can depend strongly on viral genotype , although host genotype, infection mode, and environment are likely to be important as well.
furthermore, the design and interpretation of studies of host specificity, host adaptation, and virus dispersal depend on an adequate baseline of standing genetic variation and population structure.
this is in part due to a reliance on pcr-based methods for high-throughput detection and quantification of viruses, the validity of which depends on an adequate understanding of sequence variation for primer design.
more fundamentally, the components of population-genetic variance themselves can provide important biological insights on viral demography, mutation rate, patterns of selection, and modes of replication.
in this study, we mined existing illumina sequencing data to evaluate genomic variation within two well-characterized and economically significant rna viruses of honey bees, deformed wing virus   <cit>  and israel acute paralysis virus   <cit> .
both viruses are in the picornavirus-like clade of single-strand rna viruses.
each virus is part of a larger complex of related strains or species that also infect honey bees: dwv and varroa destructor virus  <dig>  are in the iflaviridae; iapv, acute bee paralysis virus , and kashmir bee virus  are in the dicistroviridae.
dwv is common in managed honey bee colonies, is associated with a crumpled wing phenotype in highly infected adults, and, in addition to horizontal and vertical transmission among bees, is vectored by varroa mites  <cit> .
the titer of the replicating  strand of dwv in mites has been shown to be predictive of dwv-associated pathology in bees  <cit> .
iapv is the most recently described honey bee dicistrovirus and can be highly virulent  <cit> .
iapv was implicated in ccd in early work  <cit>  but this association has not recurred in subsequent studies  <cit> .
our results buttress previous studies which found a phylogenetic distinctiveness of iapv in north america relative to other continents  <cit>  and a likely selective sweep underlying reduced genetic variation among dwv isolates  <cit> .
in addition, sliding-window analyses of single nucleotide polymorphisms  indicated that regions under distinct selection pressures are differentiable from the genome background.
snp-diversity spectra differed among samples and may reflect factors such as the infection titer and number of hosts in each sample.
collectively, the data support the use of deep sequencing to investigate population-genetic variation in these viruses, ultimately to correlate viral genotypes with infection outcome.
methods
samples
we analyzed seven illumina cdna sequencing runs, each a single lane on an illumina gaii genome analyzer, with read lengths of 67– <dig> bp.
two samples, termed ccd- and ccd+, are described in  <cit> .
briefly, eight workers each were collected from  <dig> hives diagnosed with ccd and from  <dig> hives that were nominally healthy.
these hives were from temporary commercial apiaries in the united states.
another sample  derived from abdominal tissues  with evident scarring, malformation, or discoloration of the sting gland, pyloric valve, or malpighian tubules.
the weav sample represents six worker bees collected from bee weaver apiaries, navasota, texas.
the larv sample represents ten larvae from a research apiary at the university of georgia.
the brl sample represents  <dig> pupae from the beltsville bee research lab apiary that had been inoculated with a stock solution containing high levels of iapv and dwv .
the varroa sample represents a pool of approximately  <dig>  phoretic female mites that were collected from this same research apiary by dislodging them from adult worker bees in the hive.
these samples, which were collected and sequenced for other research purposes, represent a diverse set of bee life-history stages and different numbers of contributing individuals.
sample preparation and sequencing
details for sample rna extraction varied among the different samples, but all used the trizol reagent and manufacturer's recommended protocol.
cdna synthesis was performed with superscript ii  and primed with combinations of short random oligonucleotides and/or poly-dt primers.
methodological variation is a potentially important cause of variation in coverage across a viral genome in different samples, but ultimately is not germane to the goals of this study, which investigated polymorphism only at those sites with sufficient coverage irrespective of why coverage may have varied among sites or across samples.
sequence analysis
sequencing accessions are aggregated under ncbi bioproject prjna <dig>  and individual accessions are given in table  <dig>  reads were mapped with bwa  <cit>  using a polymorphism limit of 4% and a seed length of  <dig> bp.
consensus sequences were generated from the resulting alignments by identifying the most frequent base .
if a position had zero coverage, the ambiguous character 'n' was used.
indels were not permitted.
consensus-sequence analysis of nucleotide diversity , the number of segregating sites , and fu and li’s d were calculated with windows of  <dig> bases and a 25-base step, using dnasp  <cit> .
phylogenetic trees were created with mega <dig>  <cit>  using the tamura 3-parameter nucleotide model and the neighbor-joining algorithm; alternative methods produced concordant phylogenies.
bootstrap values of trees were based on  <dig>  replicates.
dwv= deformed wing virus, vdv- <dig> = varroa destructor virus- <dig>  iapv = israeli acute paralysis virus, abpv = acute bee paralysis virus, kbv = kashmir bee virus.
*rna from virus-enriched sample, see materials and methods.
snp haplotype diversity was calculated as the gene diversity at each position having at least 10x coverage, using the formula )*2) where pi is the frequency of each of the four nucleotides at that position and n is the total number of reads mapped at that position  <cit> .
haplotype diversity at a site depends only on the frequency of each base and is mathematically equivalent to haplotype diversity for single positions and to expected heterozygosity for diploid populations.
snp haplotype diversity was also calculated for sliding windows of  <dig> bases and a 25-base step.
fst values were calculated for sliding windows of  <dig> bp with a step size of  <dig> bp, using popoolation <dig>  <cit> .
a distribution of fst values for random 100-bp windows was estimated by randomly re-sampling sites to generate  <dig>  artificial intervals for each virus.
RESULTS
samples and read mapping
the seven illumina short-read data sets analyzed are listed in table  <dig>  two runs  were from a survey of worker bees from colonies diagnosed with ccd or nominally healthy controls, respectively.
these samples were collected across a broad geographic range in the u.s. as part of targeted surveys for known bee pathogens  <cit> .
one run  combined three tissues  that exhibited gross pathologies or abnormalities after dissection from worker bees .
another run  was of pupae intentionally inoculated with iapv and dwv at the beltsville bee research laboratory.
the remaining sequencing runs did not have overt disease or a priori expectation of viral infection: a) workers collected from bee weaver apiaries, navasota, texas, u.s. ; b) larvae collected from a research apiary at the university of georgia ; and c) varroa mites collected from the beltsville bee research laboratory.
read mapping was initially done competitively  for dwv, iapv, and their close relatives using parameters described in materials and methods.
this was done to determine the predominant viral species present in each sample and to evaluate whether sequence reads originating from a related species might be attributed incorrectly to the species of interest.
for dwv, reads were mapped to dwv accession nc_ <dig>  and vdv- <dig> accession nc_ <dig> .
for iapv, reads were mapped to iapv accession nc_ <dig>  initially, abpv accession nc_ <dig> , and kbv accession nc_ <dig> .
the number of reads mapping to each accession for each sample is shown in table  <dig>  although vdv- <dig> and recombinants between vdv- <dig> and dwv have been found in honey bee samples from europe  <cit> , we detected zero or a trivial number of reads for vdv- <dig> in all samples, including the varroa sample.
dwv was present in all seven samples, but with low coverage in the ccd- sample.
five samples had reads mapping to iapv, one of which  had appreciable numbers of reads mapping to the related virus kbv, but at a much lower frequency than iapv .
the presence of kbv reads could potentially bias snp estimation for iapv in this sample, but given that the two accessions listed above have a genome-average identity of only  <dig> %, we believe this effect should be negligible.
phylogeny and diversity of viral samples from consensus sequences
consensus sequences of dwv and iapv were generated for each sample in order to estimate their phylogenetic relationship with each other and with additional accessions obtained from genbank.
neighbor-joining, minimum evolution, and maximum parsimony methods produced broadly concordant tree topologies, as did different models of nucleotide mutation .
most dwv sequences from the u.s. clustered apart from the reference genome , which was isolated in italy, at 100% bootstrap support .
the ccd- sample was the exception, but this sample also had the lowest coverage and a number of gaps where no consensus could be inferred.
most u.s. sequences are very similar, as indicated by their short branch lengths, including the consensus sequence derived from the varroa sample.
compared with dwv, the branch lengths of the iapv phylogeny  are substantially longer .
there is also bootstrap support for distinct clades within iapv, in contrast to dwv.
palacios et al.
<cit>  had previously identified comparable genetic structure in iapv; specifically, they inferred three distinct lineages and found that these subdivisions were robust when different genomic regions were examined separately.
even so, all north american sequences cluster together at 100% bootstrap support relative to accessions from other continents.
this was true regardless of whether accession nc_ <dig> , an israeli isolate  <cit>  present in the ncbi genomes database, or u.s. accession eu <dig>  was used as the mapping reference to generate consensus sequences.
we next examined snp variation among sample consensus sequences and genbank whole-genome accessions, to assess how these patterns diverge from their expected values under neutral evolution.
in the absence of recombination, single-molecule viral genomes are expected to behave as one locus with respect to natural selection, showing only stochastic variation among windows.
however, recombination has been reported for both of these viruses  <cit> , so we first investigated whether the set of consensus sequences showed evidence of recombination, using hudson’s test as implemented in dnasp  <cit> .
hudson’s test indicated very low rates of recombination for dwv  but measureable rates of recombination for iapv .
these tests and the results to follow imply that variation in selective pressures among genomic regions can indeed contribute to local deviations in population-genetic parameters.
in fact, we suspect the lack of evidence for recombination in dwv in these samples may be due to the low overall diversity rather than a genuine lack of recombination, given that recombinants between dwv and the related virus vdv- <dig> have been detected  <cit> .
within-sample diversity
the results of the previous section were based on the consensus base at each position in each sample, and did not consider variation segregating within a sample.
within-sample patterns of polymorphism give further power for identifying patterns of natural selection.
in this section, we calculate haplotype diversity at snps within each sample.
haplotype diversity is numerically equivalent to expected heterozygosity in the case of diploid organisms, and is used here instead of π because the latter cannot be readily calculated from short-read data because the linkage of polymorphisms is generally unknown.
however, an ascertainment bias is introduced to snp detection by variation in the depth of sequence coverage within each sample, because the chance of identifying an alternative allele, i.e.
classifying a genomic position as a snp, increases with the number of reads mapped to that position.
we therefore calculated haplotype diversity for sliding-window analysis using only sites with a minimum coverage of ten reads.
to compare the genome-mean values between the two viral species, we further restricted sites to those that also had haplotype diversity >  <dig> , in order to minimize the contribution of low-frequency alleles.
genome scans of fst
while figure  <dig> can provide some hints as to viral regions that are more or less constrained based on mean polymorphism levels, genome scans of fst provide stronger evidence of sites that are either more conserved or more divergent than the genome average.
fst is a sensitive index of allele-frequency divergence that is scaled from zero to one ; because of this normalization, species with low overall polymorphism can nonetheless have high fst among populations , or vice versa.
one complication, however, is that the minimum number of reads per site must be met for all samples analyzed, such that the number of shared snps declines rapidly as the minimum coverage requirement increases.
we therefore used a fairly lenient coverage requirement for our sliding-window estimates: at least  <dig> reads per site and at least 40% of the sites in the window meeting this criterion.
increasing the coverage requirement did not change the qualitative pattern for either genome , but greatly increased the number of windows for which no value could be calculated .
to evaluate the significance of local variation in fst, we randomized genomic positions to generate  <dig>  100-bp windows for each virus.
values less than the  <dig>  percentile or greater than the  <dig>  percentile are below or above the dotted lines in figure  <dig>  respectively.
the two viruses show strikingly different patterns of fst across their genomes.
dwv oscillates frequently between a higher and a lower band of values, such that the histogram  is strongly bimodal, with one mean around  <dig>  and another mean around  <dig> .
this pattern appears to be driven by genetic drift in the brl and weav samples, which have the longest branch lengths in figure  <dig>  low overall diversity compared with the other dwv samples , and pairwise per-window fst values that are much higher on average than among the remaining samples.
excluding these two samples, the mean fst values among the remaining dwv samples are much lower and less variable  than when all seven populations are included .
while diversifying selection is another cause of high fst, the large number of sites with high fst and the similarity of the values indicates that genetic drift, perhaps due to a founder effect, is the more parsimonious interpretation.
it is nonetheless surprising to see such rapid swings in local fst, and it suggests that strongly constrained sites occur in numerous small clusters throughout the genome.
the dotted lines in figure 9a represent the upper and lower  <dig> % of the distribution of randomized fst values, calculated after excluding the brl and weav samples.
the sliding-window fst estimates for the remaining five samples  fall within the extremes of this distribution, suggesting predominantly stochastic variation.
in contrast to dwv, iapv has a unimodal distribution of fst values across the genome and a higher mean value of approximately  <dig>  , although several strong outliers of low fst are apparent.
this pattern suggests that there are only a few windows within which purifying selection is strong and consistent across samples, resulting in local haplotype frequencies that are far more uniform than the genome average.
CONCLUSIONS
in this study, we have performed an initial, random survey of genomic variation in two honey bee viruses, dwv and iapv, of established relevance to honey bee health.
even though the sequencing was originally undertaken for diverse, unrelated reasons, our results demonstrate the utility of population-genomic data for detecting important aspects of viral biology.
more directed sequencing efforts using virus-enriched samples and controlled study designs should provide even greater power to discern evolutionary forces acting on genomic variation.
the data presented provide support for the idea that the spread of the varroa mite has enforced a selective sweep on dwv sequences, by giving high fitness to genotypes that can replicate in and be transferred by the parasite.
dwv replication is known to occur in varroa mites  <cit>  and the introduction of varroa to hawaii was associated with both a great increase in dwv infections and a great decrease in dwv diversity in that state  <cit> .
the latter data strongly implicate selection for vectoring by varroa as a driver of genetic uniformity in dwv.
this hypothesized selective sweep is consistent with the lower diversity and divergence found among populations of dwv relative to iapv, despite comparable levels of within-population variation.
while there is correlational evidence that iapv can be transmitted by varroa mites  <cit> , the lower frequency of iapv relative to dwv and greater genetic distances among iapv samples suggest that varroa is not a frequent vector of iapv in natural populations.
the high read count of dwv in the varroa sample  is consistent with replication in the varroa mite, although the various methods used to generate these data preclude a reliable reference gene for relative quantification across samples or even determining whether the virus was present within tissues other than the gastrointestinal tract.
interestingly, the consensus sequence and polymorphism levels of this sample were very similar to honey-bee derived samples, implying that no large-scale genotypic shift occurs during replication in that host.
more definitive evidence that isolates from either host have equivalent fitness in the other is needed, however.
in the absence of a culture method for dwv, a useful approach would be to analyze viral growth rates and sequence divergence among replicated infections in each host, initiated from natural inoculates from the same or alternate host.
such data would demonstrate a replicable effect of host source on the fitness of the virus.
a similar approach could be taken using inoculates from the few remaining areas of honey-bee cultivation that are free of the varroa mite.
it would also be of interest to investigate dwv diversity in africanized populations, which experience much less parasitism by varroa.
as some honey bee rna viruses have been detected by pcr in other sympatric pollinators  <cit> , indicating the possibility of alternate insect hosts, deep sequencing can also shed light on the process of adaptation to alternative hosts.
in macroscopic eukaryotes that do not disperse widely, high fst at neutral loci is expected to arise from a lack of gene flow between separate interbreeding populations.
this interpretation is dubious for rna viruses, however, which can rapidly reach selection-mutation equilibrium even after single-genome bottlenecks.
thus, high or low fst values more likely indicate selective pressures rather than dispersal rates of virions.
multiple narrow windows of low fst values were strikingly evident among the iapv samples, for example, suggesting strong and concordant selective pressures shaping polymorphism levels at those sites.
even so, the genome-wide pattern of higher fst between dwv sequences from the weav and brl samples relative to other samples suggests that founder effects do occur and can substantially impact the spectrum of genotypes observed.
in fact, phylogenetic analysis of both species shows geographic structuring of variation at large scales, as u.s. samples tend to cluster together relative to samples from other continents.
this continental-scale genetic structure has apparently persisted despite historical trade in bees and/or bee products that likely moved viral particles between continents.
iapv exhibits phylogenetic structure within the u.s. as well, but it is not yet known whether there is also a spatial component of this variation.
both genomes had low among-sample diversity at the 3’ end of the alignments , and an example of strong conservation in an approximately forty-base window was shown in additional file 1: figure s <dig>  the 3’noncoding region of the picornavirus genome is believed to be a major cis-acting regulator of negative-strand rna synthesis by promoting recognition of the replication complex  <cit> .
an important role in genome replication could therefore underlie the low level of nucleotide variation in this region, but infectious clones and cell-culture propagation methods are presently unavailable to support reverse genetic studies of dwv.
interestingly, the highest haplotype diversities at individual iapv snps occurred at the extreme 5’ and 3’ sites ; the latter do not influence figure  <dig> because of alignment gaps for some iapv accessions and consensus sequences.
the greater diversity at these sites could potentially indicate higher rates of mutation associated with the replication of genome ends.
a distinctive feature of dwv is a peak of sequence diversity at the 5’ end of the coding sequence that was observed both among populations  and within populations .
in fact, fu and li’s d was weakly positive upstream of the coding region rather than negative, contrary to what might be expected based on its putative role in binding the host ribosome  <cit> .
it is not clear whether windows that did have significantly negative values of d are experiencing stronger purifying selection than other regions or are simply the extremes of a distribution that has been shifted into predominantly negative values.
it should also be noted that this statistic was applied to consensus sequences rather than haplotypes; the latter would likely provide greater resolution of the extent of among-region variation in selective pressures.
however, our results help define interesting regions to target for amplicon-based deep sequencing with longer or paired-end reads, so that haplotype structure can be retained.
competing interests
the authors declare they have no competing interests.
authors’ contributions
rsc conceived the study, performed population genomic analyses, and drafted the manuscript.
hb and bd participated in the design of the analysis and interpretation of the results.
hb, dv, dw, yc, and je collected samples and contributed sequencing data.
all authors read and approved the final manuscript.
supplementary material
additional file 1: figure s1
a. an alignment of a portion of the 3’ untranslated region of dwv, showing conservation among all consensus sequences and accessions, including vdv- <dig>  the ccd- and ccd+ samples are not shown due to insufficient read coverage in this region.
b. mean haplotype diversity  in the 3’ untranslated region of dwv, showing a lack of within-population variation in the same window  that is marked in panel a.
click here for file
additional file 2: figure s2
histograms of fst values from the sliding-window analysis in figure  <dig>  for dwv, the values are for all seven positive samples including the outlier samples brl and weav.
click here for file
acknowledgments
we thank keith delaplane and jennifer berry for providing larvae from the university of georgia apiary.
dawn lopez prepared samples for sequencing.
this project was supported by the coordinate agricultural project on pollinator health  and the national honey board.

 in addition to producing honey and other hive products, managed honey bee colonies are vital for the pollination of many crops. pests and diseases that impact honey bees are therefore of economic concern, particularly in light of the long-distance transport of bees, bee products, and beekeeping equipment that can facilitate their spread. high rates of unexplained colony loss in recent years have heightened these concerns and provided impetus for surveys of geographic and genetic variation in honey bee pathogens. a number of rna viruses have been recognized as major disease agents of honey bees. in addition to the recognized morbidity and mortality that can result from overt viral infections, lower viral loads that do not produce visible symptoms might nonetheless impact colony health as well. for example, viral loads are negatively associated with colony survival and there is correlational evidence of synergism among viruses in colony collapse disorder. interactions have also been identified between viruses and other parasites, particularly the ectoparasitic mite varroa destructor. recent genomic studies have identified novel viruses that may also be relevant to colony health, but about which little is yet known. additionally, recombination between co-infecting rna viruses may facilitate adaptation or the emergence of novel viruses. it has therefore become a major goal to better understand the diversity of viral species in honey bee populations and how they may interact during co-infection. understanding the genetic diversity within a viral species is important as well. rna viruses in particular have high mutation rates, can evolve rapidly, and may exhibit adaptation to specific hosts and/or tissues. virulence can depend strongly on viral genotype, although host genotype, infection mode, and environment are likely to be important as well. furthermore, the design and interpretation of studies of host specificity, host adaptation, and virus dispersal depend on an adequate baseline of standing genetic variation and population structure. this is in part due to a reliance on pcr-based methods for high-throughput detection and quantification of viruses, the validity of which depends on an adequate understanding of sequence variation for primer design. more fundamentally, the components of population-genetic variance themselves can provide important biological insights on viral demography, mutation rate, patterns of selection, and modes of replication. in this study, we mined existing illumina sequencing data to evaluate genomic variation within two well-characterized and economically significant rna viruses of honey bees, deformed wing virus and israel acute paralysis virus. both viruses are in the picornavirus-like clade of single-strand rna viruses. each virus is part of a larger complex of related strains or species that also infect honey bees: dwv and varroa destructor virus are in the iflaviridae; iapv, acute bee paralysis virus, and kashmir bee virus are in the dicistroviridae. dwv is common in managed honey bee colonies, is associated with a crumpled wing phenotype in highly infected adults, and, in addition to horizontal and vertical transmission among bees, is vectored by varroa mites. the titer of the replicating strand of dwv in mites has been shown to be predictive of dwv-associated pathology in bees. iapv is the most recently described honey bee dicistrovirus and can be highly virulent. iapv was implicated in ccd in early work but this association has not recurred in subsequent studies. our results buttress previous studies which found a phylogenetic distinctiveness of iapv in north america relative to other continents and a likely selective sweep underlying reduced genetic variation among dwv isolates. in addition, sliding-window analyses of single nucleotide polymorphisms indicated that regions under distinct selection pressures are differentiable from the genome background. snp-diversity spectra differed among samples and may reflect factors such as the infection titer and number of hosts in each sample. collectively, the data support the use of deep sequencing to investigate population-genetic variation in these viruses, ultimately to correlate viral genotypes with infection outcome. methods samples we analyzed seven illumina cdna sequencing runs, each a single lane on an illumina gaii genome analyzer, with read lengths of 67– bp. two samples, termed ccd- and ccd+, are described in. briefly, eight workers each were collected from hives diagnosed with ccd and from hives that were nominally healthy. these hives were from temporary commercial apiaries in the united states. another sample derived from abdominal tissues with evident scarring, malformation, or discoloration of the sting gland, pyloric valve, or malpighian tubules. the weav sample represents six worker bees collected from bee weaver apiaries, navasota, texas. the larv sample represents ten larvae from a research apiary at the university of georgia. the brl sample represents pupae from the beltsville bee research lab apiary that had been inoculated with a stock solution containing high levels of iapv and dwv. the varroa sample represents a pool of approximately phoretic female mites that were collected from this same research apiary by dislodging them from adult worker bees in the hive. these samples, which were collected and sequenced for other research purposes, represent a diverse set of bee life-history stages and different numbers of contributing individuals. sample preparation and sequencing details for sample rna extraction varied among the different samples, but all used the trizol reagent and manufacturer's recommended protocol. cdna synthesis was performed with superscript ii and primed with combinations of short random oligonucleotides and/or poly-dt primers. methodological variation is a potentially important cause of variation in coverage across a viral genome in different samples, but ultimately is not germane to the goals of this study, which investigated polymorphism only at those sites with sufficient coverage irrespective of why coverage may have varied among sites or across samples. sequence analysis sequencing accessions are aggregated under ncbi bioproject prjna and individual accessions are given in table reads were mapped with bwa using a polymorphism limit of 4% and a seed length of bp. consensus sequences were generated from the resulting alignments by identifying the most frequent base. if a position had zero coverage, the ambiguous character'n' was used. indels were not permitted. consensus-sequence analysis of nucleotide diversity, the number of segregating sites, and fu and li’s d were calculated with windows of bases and a 25-base step, using dnasp. phylogenetic trees were created with mega using the tamura 3-parameter nucleotide model and the neighbor-joining algorithm; alternative methods produced concordant phylogenies. bootstrap values of trees were based on replicates. dwv= deformed wing virus, vdv- = varroa destructor virus- iapv = israeli acute paralysis virus, abpv = acute bee paralysis virus, kbv = kashmir bee virus. *rna from virus-enriched sample, see materials and methods. snp haplotype diversity was calculated as the gene diversity at each position having at least 10x coverage, using the formula )*2) where pi is the frequency of each of the four nucleotides at that position and n is the total number of reads mapped at that position. haplotype diversity at a site depends only on the frequency of each base and is mathematically equivalent to haplotype diversity for single positions and to expected heterozygosity for diploid populations. snp haplotype diversity was also calculated for sliding windows of bases and a 25-base step. fst values were calculated for sliding windows of bp with a step size of bp, using popoolation. a distribution of fst values for random 100-bp windows was estimated by randomly re-sampling sites to generate artificial intervals for each virus.  samples and read mapping the seven illumina short-read data sets analyzed are listed in table two runs were from a survey of worker bees from colonies diagnosed with ccd or nominally healthy controls, respectively. these samples were collected across a broad geographic range in the u.s. as part of targeted surveys for known bee pathogens. one run combined three tissues that exhibited gross pathologies or abnormalities after dissection from worker bees. another run was of pupae intentionally inoculated with iapv and dwv at the beltsville bee research laboratory. the remaining sequencing runs did not have overt disease or a priori expectation of viral infection: a) workers collected from bee weaver apiaries, navasota, texas, u.s.; b) larvae collected from a research apiary at the university of georgia; and c) varroa mites collected from the beltsville bee research laboratory. read mapping was initially done competitively for dwv, iapv, and their close relatives using parameters described in materials and methods. this was done to determine the predominant viral species present in each sample and to evaluate whether sequence reads originating from a related species might be attributed incorrectly to the species of interest. for dwv, reads were mapped to dwv accession nc_ and vdv- accession nc_. for iapv, reads were mapped to iapv accession nc_ initially, abpv accession nc_, and kbv accession nc_. the number of reads mapping to each accession for each sample is shown in table although vdv- and recombinants between vdv- and dwv have been found in honey bee samples from europe, we detected zero or a trivial number of reads for vdv- in all samples, including the varroa sample. dwv was present in all seven samples, but with low coverage in the ccd- sample. five samples had reads mapping to iapv, one of which had appreciable numbers of reads mapping to the related virus kbv, but at a much lower frequency than iapv. the presence of kbv reads could potentially bias snp estimation for iapv in this sample, but given that the two accessions listed above have a genome-average identity of only %, we believe this effect should be negligible. phylogeny and diversity of viral samples from consensus sequences consensus sequences of dwv and iapv were generated for each sample in order to estimate their phylogenetic relationship with each other and with additional accessions obtained from genbank. neighbor-joining, minimum evolution, and maximum parsimony methods produced broadly concordant tree topologies, as did different models of nucleotide mutation. most dwv sequences from the u.s. clustered apart from the reference genome, which was isolated in italy, at 100% bootstrap support. the ccd- sample was the exception, but this sample also had the lowest coverage and a number of gaps where no consensus could be inferred. most u.s. sequences are very similar, as indicated by their short branch lengths, including the consensus sequence derived from the varroa sample. compared with dwv, the branch lengths of the iapv phylogeny are substantially longer. there is also bootstrap support for distinct clades within iapv, in contrast to dwv. palacios et al. had previously identified comparable genetic structure in iapv; specifically, they inferred three distinct lineages and found that these subdivisions were robust when different genomic regions were examined separately. even so, all north american sequences cluster together at 100% bootstrap support relative to accessions from other continents. this was true regardless of whether accession nc_, an israeli isolate present in the ncbi genomes database, or u.s. accession eu was used as the mapping reference to generate consensus sequences. we next examined snp variation among sample consensus sequences and genbank whole-genome accessions, to assess how these patterns diverge from their expected values under neutral evolution. in the absence of recombination, single-molecule viral genomes are expected to behave as one locus with respect to natural selection, showing only stochastic variation among windows. however, recombination has been reported for both of these viruses, so we first investigated whether the set of consensus sequences showed evidence of recombination, using hudson’s test as implemented in dnasp. hudson’s test indicated very low rates of recombination for dwv but measureable rates of recombination for iapv. these tests and the results to follow imply that variation in selective pressures among genomic regions can indeed contribute to local deviations in population-genetic parameters. in fact, we suspect the lack of evidence for recombination in dwv in these samples may be due to the low overall diversity rather than a genuine lack of recombination, given that recombinants between dwv and the related virus vdv- have been detected. within-sample diversity the results of the previous section were based on the consensus base at each position in each sample, and did not consider variation segregating within a sample. within-sample patterns of polymorphism give further power for identifying patterns of natural selection. in this section, we calculate haplotype diversity at snps within each sample. haplotype diversity is numerically equivalent to expected heterozygosity in the case of diploid organisms, and is used here instead of π because the latter cannot be readily calculated from short-read data because the linkage of polymorphisms is generally unknown. however, an ascertainment bias is introduced to snp detection by variation in the depth of sequence coverage within each sample, because the chance of identifying an alternative allele, i.e. classifying a genomic position as a snp, increases with the number of reads mapped to that position. we therefore calculated haplotype diversity for sliding-window analysis using only sites with a minimum coverage of ten reads. to compare the genome-mean values between the two viral species, we further restricted sites to those that also had haplotype diversity >, in order to minimize the contribution of low-frequency alleles. genome scans of fst while figure can provide some hints as to viral regions that are more or less constrained based on mean polymorphism levels, genome scans of fst provide stronger evidence of sites that are either more conserved or more divergent than the genome average. fst is a sensitive index of allele-frequency divergence that is scaled from zero to one; because of this normalization, species with low overall polymorphism can nonetheless have high fst among populations, or vice versa. one complication, however, is that the minimum number of reads per site must be met for all samples analyzed, such that the number of shared snps declines rapidly as the minimum coverage requirement increases. we therefore used a fairly lenient coverage requirement for our sliding-window estimates: at least reads per site and at least 40% of the sites in the window meeting this criterion. increasing the coverage requirement did not change the qualitative pattern for either genome, but greatly increased the number of windows for which no value could be calculated. to evaluate the significance of local variation in fst, we randomized genomic positions to generate 100-bp windows for each virus. values less than the percentile or greater than the percentile are below or above the dotted lines in figure respectively. the two viruses show strikingly different patterns of fst across their genomes. dwv oscillates frequently between a higher and a lower band of values, such that the histogram is strongly bimodal, with one mean around and another mean around. this pattern appears to be driven by genetic drift in the brl and weav samples, which have the longest branch lengths in figure low overall diversity compared with the other dwv samples, and pairwise per-window fst values that are much higher on average than among the remaining samples. excluding these two samples, the mean fst values among the remaining dwv samples are much lower and less variable than when all seven populations are included. while diversifying selection is another cause of high fst, the large number of sites with high fst and the similarity of the values indicates that genetic drift, perhaps due to a founder effect, is the more parsimonious interpretation. it is nonetheless surprising to see such rapid swings in local fst, and it suggests that strongly constrained sites occur in numerous small clusters throughout the genome. the dotted lines in figure 9a represent the upper and lower % of the distribution of randomized fst values, calculated after excluding the brl and weav samples. the sliding-window fst estimates for the remaining five samples fall within the extremes of this distribution, suggesting predominantly stochastic variation. in contrast to dwv, iapv has a unimodal distribution of fst values across the genome and a higher mean value of approximately, although several strong outliers of low fst are apparent. this pattern suggests that there are only a few windows within which purifying selection is strong and consistent across samples, resulting in local haplotype frequencies that are far more uniform than the genome average.  in this study, we have performed an initial, random survey of genomic variation in two honey bee viruses, dwv and iapv, of established relevance to honey bee health. even though the sequencing was originally undertaken for diverse, unrelated reasons, our results demonstrate the utility of population-genomic data for detecting important aspects of viral biology. more directed sequencing efforts using virus-enriched samples and controlled study designs should provide even greater power to discern evolutionary forces acting on genomic variation. the data presented provide support for the idea that the spread of the varroa mite has enforced a selective sweep on dwv sequences, by giving high fitness to genotypes that can replicate in and be transferred by the parasite. dwv replication is known to occur in varroa mites and the introduction of varroa to hawaii was associated with both a great increase in dwv infections and a great decrease in dwv diversity in that state. the latter data strongly implicate selection for vectoring by varroa as a driver of genetic uniformity in dwv. this hypothesized selective sweep is consistent with the lower diversity and divergence found among populations of dwv relative to iapv, despite comparable levels of within-population variation. while there is correlational evidence that iapv can be transmitted by varroa mites, the lower frequency of iapv relative to dwv and greater genetic distances among iapv samples suggest that varroa is not a frequent vector of iapv in natural populations. the high read count of dwv in the varroa sample is consistent with replication in the varroa mite, although the various methods used to generate these data preclude a reliable reference gene for relative quantification across samples or even determining whether the virus was present within tissues other than the gastrointestinal tract. interestingly, the consensus sequence and polymorphism levels of this sample were very similar to honey-bee derived samples, implying that no large-scale genotypic shift occurs during replication in that host. more definitive evidence that isolates from either host have equivalent fitness in the other is needed, however. in the absence of a culture method for dwv, a useful approach would be to analyze viral growth rates and sequence divergence among replicated infections in each host, initiated from natural inoculates from the same or alternate host. such data would demonstrate a replicable effect of host source on the fitness of the virus. a similar approach could be taken using inoculates from the few remaining areas of honey-bee cultivation that are free of the varroa mite. it would also be of interest to investigate dwv diversity in africanized populations, which experience much less parasitism by varroa. as some honey bee rna viruses have been detected by pcr in other sympatric pollinators, indicating the possibility of alternate insect hosts, deep sequencing can also shed light on the process of adaptation to alternative hosts. in macroscopic eukaryotes that do not disperse widely, high fst at neutral loci is expected to arise from a lack of gene flow between separate interbreeding populations. this interpretation is dubious for rna viruses, however, which can rapidly reach selection-mutation equilibrium even after single-genome bottlenecks. thus, high or low fst values more likely indicate selective pressures rather than dispersal rates of virions. multiple narrow windows of low fst values were strikingly evident among the iapv samples, for example, suggesting strong and concordant selective pressures shaping polymorphism levels at those sites. even so, the genome-wide pattern of higher fst between dwv sequences from the weav and brl samples relative to other samples suggests that founder effects do occur and can substantially impact the spectrum of genotypes observed. in fact, phylogenetic analysis of both species shows geographic structuring of variation at large scales, as u.s. samples tend to cluster together relative to samples from other continents. this continental-scale genetic structure has apparently persisted despite historical trade in bees and/or bee products that likely moved viral particles between continents. iapv exhibits phylogenetic structure within the u.s. as well, but it is not yet known whether there is also a spatial component of this variation. both genomes had low among-sample diversity at the 3’ end of the alignments, and an example of strong conservation in an approximately forty-base window was shown in additional file 1: figure s the 3’noncoding region of the picornavirus genome is believed to be a major cis-acting regulator of negative-strand rna synthesis by promoting recognition of the replication complex. an important role in genome replication could therefore underlie the low level of nucleotide variation in this region, but infectious clones and cell-culture propagation methods are presently unavailable to support reverse genetic studies of dwv. interestingly, the highest haplotype diversities at individual iapv snps occurred at the extreme 5’ and 3’ sites; the latter do not influence figure because of alignment gaps for some iapv accessions and consensus sequences. the greater diversity at these sites could potentially indicate higher rates of mutation associated with the replication of genome ends. a distinctive feature of dwv is a peak of sequence diversity at the 5’ end of the coding sequence that was observed both among populations and within populations. in fact, fu and li’s d was weakly positive upstream of the coding region rather than negative, contrary to what might be expected based on its putative role in binding the host ribosome. it is not clear whether windows that did have significantly negative values of d are experiencing stronger purifying selection than other regions or are simply the extremes of a distribution that has been shifted into predominantly negative values. it should also be noted that this statistic was applied to consensus sequences rather than haplotypes; the latter would likely provide greater resolution of the extent of among-region variation in selective pressures. however, our results help define interesting regions to target for amplicon-based deep sequencing with longer or paired-end reads, so that haplotype structure can be retained. competing interests the authors declare they have no competing interests. authors’ contributions rsc conceived the study, performed population genomic analyses, and drafted the manuscript. hb and bd participated in the design of the analysis and interpretation of the results. hb, dv, dw, yc, and je collected samples and contributed sequencing data. all authors read and approved the final manuscript. supplementary material additional file 1: figure s1 a. an alignment of a portion of the 3’ untranslated region of dwv, showing conservation among all consensus sequences and accessions, including vdv- the ccd- and ccd+ samples are not shown due to insufficient read coverage in this region. b. mean haplotype diversity in the 3’ untranslated region of dwv, showing a lack of within-population variation in the same window that is marked in panel a. click here for file additional file 2: figure s2 histograms of fst values from the sliding-window analysis in figure for dwv, the values are for all seven positive samples including the outlier samples brl and weav. click here for file acknowledgments we thank keith delaplane and jennifer berry for providing larvae from the university of georgia apiary. dawn lopez prepared samples for sequencing. this project was supported by the coordinate agricultural project on pollinator health and the national honey board.
consensus sequences of dwv for each sample showed little phylogenetic divergence, low nucleotide diversity, and strongly negative values of fu and li’s d statistic, suggesting a recent population bottleneck and/or purifying selection. in this study, we mined existing illumina sequence reads to investigate single-nucleotide polymorphisms within two rna viruses of the western honey bee, deformed wing virus and israel acute paralysis virus. apis melliferadeformed wing virusisrael acute paralysis virussingle nucleotide polymorphismpopulation genomics the kakugo strain of dwv fell outside of all other dwv sequences at 100% bootstrap support. fst values were bimodal for dwv, likely reflecting neutral divergence in two low-diversity populations, whereas iapv had several sites that were strong outliers with very low fst. iapv consensus sequences supported the existence of multiple clades as had been previously reported, and fu and li’s d was closer to neutral expectation overall, although a sliding-window analysis identified a significantly positive d within the protease region, suggesting selection maintains diversity in that region. deep sequencing of viruses isolated from infected hosts is an efficient way to measure population-genetic variation and can reveal patterns of dispersal and natural selection. this initial survey of genetic variation within honey bee rna viruses suggests future directions for studies examining the underlying causes of population-genetic structure in these economically important pathogens. all viral rna was extracted from north american samples of honey bees or, in one case, the ectoparasitic mite varroa destructor.
 deep sequencing of viruses isolated from infected hosts is an efficient way to measure population-genetic variation and can reveal patterns of dispersal and natural selection. in this study, we mined existing illumina sequence reads to investigate single-nucleotide polymorphisms within two rna viruses of the western honey bee, deformed wing virus and israel acute paralysis virus. all viral rna was extracted from north american samples of honey bees or, in one case, the ectoparasitic mite varroa destructor.  coverage depth was generally lower for iapv than dwv, and marked gaps in coverage occurred in several narrow regions of iapv. these coverage gaps occurred across sequencing runs and were virtually unchanged when reads were re-mapped with greater permissiveness, suggesting a recurrent sequencing artifact rather than strain divergence. consensus sequences of dwv for each sample showed little phylogenetic divergence, low nucleotide diversity, and strongly negative values of fu and li’s d statistic, suggesting a recent population bottleneck and/or purifying selection. the kakugo strain of dwv fell outside of all other dwv sequences at 100% bootstrap support. iapv consensus sequences supported the existence of multiple clades as had been previously reported, and fu and li’s d was closer to neutral expectation overall, although a sliding-window analysis identified a significantly positive d within the protease region, suggesting selection maintains diversity in that region. within-sample mean diversity was comparable between the two viruses on average, although for both viruses there was substantial variation among samples in mean diversity at third codon positions and in the number of high-diversity sites. fst values were bimodal for dwv, likely reflecting neutral divergence in two low-diversity populations, whereas iapv had several sites that were strong outliers with very low fst.  this initial survey of genetic variation within honey bee rna viruses suggests future directions for studies examining the underlying causes of population-genetic structure in these economically important pathogens. apis melliferadeformed wing virusisrael acute paralysis virussingle nucleotide polymorphismpopulation genomics 