Recurrent evolution of high virulence in isolated populations of a DNA virus
Abstract
Hosts and viruses are constantly evolving in response to each other: as a host attempts to suppress a virus, the virus attempts to evade and suppress the host’s immune system. Here, we describe the recurrent evolution of a virulent strain of a DNA virus, which infects multiple Drosophila species. Specifically, we identified two distinct viral types that differ 100-fold in viral titer in infected individuals, with similar differences observed in multiple species. Our analysis suggests that one of the viral types recurrently evolved at least four times in the past ~30,000 years, three times in Arizona and once in another geographically distinct species. This recurrent evolution may be facilitated by an effective mutation rate which increases as each prior mutation increases viral titer and effective population size. The higher titer viral type suppresses the host-immune system and an increased virulence compared to the low viral titer type.
eLife digest
Animals constantly evolve to protect themselves against viruses, and in turn, viruses evolve to escape their host’s new defenses. As a result, genes involved in this arms’ race are some of the fastest evolving in nature. A better understanding of how host-virus evolution works could help in the search for treatments for many human and animal diseases.
Repetition is one of the gold standard requirements for biological experiments. Watching different groups of animals and viruses evolve under the same conditions makes it possible for researchers to work out whether certain changes are more likely than others. This is easy to do in the laboratory, where conditions can be controlled, but much more complicated to accomplish in the wild. Wild populations are rarely completely isolated, and often face different environmental conditions. One animal-virus pair for which this is not the case is made up of the fly Drosophila innubila, and its virus Drosophila innubila nudivirus. They live in the 'sky islands' of North America, patches of forests surrounded by hundreds of kilometers of desert. These islands are like natural test tubes, isolated ecosystems each with its own separate fly and virus populations and limited gene flow between populations.
To understand how this virus-host pair evolves, Hill and Unckless sequenced the genomes of flies and viruses from four different populations. While the fly genomes did not show evidence of strong differences between populations, the virus genomes did. There were two distinct types of virus, one of which was a lot more effective than the other at infecting flies, possibly because it was better at blocking the fly's immune defenses. Unexpectedly, this virus type had evolved more than once, emerging separately on at least four different occasions. Hill and Unckless suggest that the natural interactions between flies with similar genomes and the virus guide evolution down the same path time and time again.
This work on wild populations contributes to the understanding of the evolution of viruses and their hosts. One question left unanswered is why both types of virus (one more effective at infecting the flies and the other less so) persist in each population when one is better at blocking the fly's immune response? Future work using isolated populations like these could shed more light on the pressures that shape the evolution of viruses and their hosts, potentially helping in the study of human viruses, like HIV.
Introduction
Antagonistic coevolution between hosts and their parasites is nearly ubiquitous across the diversity of life (Burt and Trivers, 2006). As a result, genes involved in immune defense are among the fastest evolving genes in host genomes (Nielsen et al., 2005; Sackton et al., 2007; Enard et al., 2016; Shultz and Sackton, 2019). Viruses are a particular fitness burden on hosts; for viruses to persist within populations, they must successfully invade the host organism, contend with the host-immune system, replicate and then transmit the newly produced particles to a new host (Holmes, 2007; Gifford, 2012). Once successfully established in a population, natural selection acts to modulate the rate the virus propagates relative to its virulence, optimizing the ratio of virulence to transmission (Williams and Nesse, 1991; May and Nowak, 1995; Lipsitch et al., 1996). Due in part to their elevated mutation rate and large population sizes, viruses can accomplish this, and often co-opt or manipulate host-pathways in the process (Burgyán and Havelda, 2011; Davey et al., 2011; Palmer et al., 2019).
Given the strong selective pressures viruses exert on their hosts (and hosts on viruses), it would be useful to assess how the two players behave in replicate evolutionary experiments. Experiments examining the evolution of pathogens and hosts are common in the lab (Perron et al., 2006; Paterson et al., 2010; Bull et al., 2011; Martins et al., 2014; Scanlan et al., 2015) and even in patients (Pennings, 2012; Pennings et al., 2014; Feder et al., 2019), but difficult in natural settings (Souza et al., 2002; Grubaugh et al., 2015), so we are rarely able to observe the coevolutionary dynamics between a virus and its host in replicate populations to examine natural coevolution and identify if evolution repeats itself. These studies are necessary as they frequently help characterize the pathways viruses have evolved to escape or suppress the host-immune system, and whether the same host pathways are common targets of the same virus (Sabin et al., 2010).
Several studies of viruses in natural populations highlight that the same initial virus can follow the same adaptive path in multiple replicate samples. This results in parallel or convergent evolution of a virus better optimized to the host, possible due to the relatively high viral mutation rate and viral population size (Bull et al., 1997; Casino et al., 1999; Crandall et al., 1999; van Mierlo et al., 2012). However, in other systems, the same initial virus can evolve divergently in a host or location-dependent manner. Specifically, virulence and the diversity of mutations can differ dramatically based on differences in the host environment, the extremely strong selective pressures acting on the virus, and the unpredictability of some viruses due to their elevated mutation rate (Martinez-Picado et al., 2002; Cuevas et al., 2003; Real et al., 2005; Grubaugh et al., 2015). Characterizing how viruses adapt to their long-term hosts in naturally structured populations will help broaden and expand our understanding of how hosts and pathogens evolve in response to each other, and how repeatable evolution is in the face of minor environmental differences.
We took advantage of a natural host/DNA virus infection model involving populations separated by hundreds of kilometers to study how viruses evolve with their hosts in replicate populations with limited gene flow. Drosophila innubila is a mycophagous species that lives in montane forests in the ‘Sky islands’ in southwestern North America. They are commonly infected with a double-stranded DNA virus, the Drosophila innubila Nudivirus (DiNV) (Hill and Unckless, 2018). A previous study examined the rates of evolution of DiNV (Hill and Unckless, 2018), finding the envelope and replication machinery to be rapidly evolving in DiNV, suggesting its importance in viral propagation (Hill and Unckless, 2018). Additionally, a viral suppressor of Toll was identified in the genome of a closely related virus and DiNV, suggesting that components of the Toll pathway (perhaps antimicrobial peptides – AMPs) interact with the virus, which these viruses have then evolved to suppress (Palmer et al., 2019). Consistent with this, genes in the Toll pathway and Toll-regulated AMPs are rapidly evolving in D. innubila (Hill et al., 2019). A more thorough examination of the natural evolution of both the host and virus is necessary to understand how DiNV and D. innubila interact with each other beyond this suppressor protein.
We surveyed natural genetic variation in DiNV in four replicate populations to infer its co-evolutionary history with D. innubila. To this end, we (a) estimated how long DiNV has infected D. innubila, (b) inferred the evolutionary history of the host and virus among the populations, including signatures of selection and recombination (c) examined host and virus genetic variation associated with viral titer within individual hosts, and (d) characterized gene expression differences associated with that genetic variation. We identified two viral multilocus genotypes that differ by 11 focal SNPs and found that these viral types are maintained within the same host population and across multiple isolated host populations. These SNPs are tightly linked, likely brought together by a combination of recurrent mutation and recombination, and this linkage appears to be maintained by strong selection. One viral type is associated with 100-fold higher viral titer and increased virulence compared to the other. Further, we found evidence that the SNPs associated with the high titer type have evolved independently in at least three geographically-isolated viral populations infecting the sympatric D. innubila and D. azteca, and a geographically separate population of viruses infecting D. falleni. Together, these results suggest rapid evolutionary dynamics of host–virus interactions, due to recurrent evolution of a highly virulent haplotype that interacts with multiple host pathways.
Results
DiNV segregates for linked variants strongly associated with viral titer
To characterize the evolutionary dynamics of wild Drosophila innubila Nudivirus (DiNV) in its host (D. innubila), we sequenced wild-caught individuals from four populations with the expectation that some (~40% in previous samples) individuals would be infected (Hill and Unckless, 2020). We considered strains to be infected with DiNV if they had at least 10x coverage for 95% of the genome. We confirmed the infection of a subsample of these using PCR for the presence of a DiNV specific gene (Supplementary Table 1). In total, we used sequencing information for 57, 92, 92, and 92 individuals from the Huachucas (HU), Santa Ritas (SR), Chiricahuas (CH), and Prescott (PR) populations with infection rates 26, 44, 71, and 79%, respectively (Supplementary file 1 - Table 1). We also sequenced 35 individual males collected in the Chiricahuas in 2001 (52% infected with DiNV) and 80 individual males collected in the Chiricahua’s in 2018 (Supplementary file 1 - Table 2, 40 infected with DiNV and 40 uninfected determined before sequencing using PCR).
Candidate viral SNPs associated with viral titer with associated genes and the functional category of that gene.
SNP locus | Nearest gene | SNP functional annotation | Nearest gene functional annotation | p-value (FDR-corrected) |
---|---|---|---|---|
G14249T | 19K/PIF-4 | Non-synonymous | Per OS Infectivity factor envelope protein, required for oral infection | 3.49e-15 (4.89e-12) |
C41210T | LEF-4 | Non-synonymous | RNA polymerase subunit for RNA modification | 7.99e-12 (1.12e-09) |
G42389T | gp83 | Upstream | Suspected virulence factor which suppresses Toll activity | 6.73–14 (9.44e-11) |
A59194G | gp51 | Intergenic | Suspected virulence factor | 2.96e-12 (4.15e-09) |
C59275A | PIF-6 | Upstream | Per OS Infectivity factor envelope protein, required for oral infection | 1.65e-18 (2.31e-15) |
T59276C | PIF-6 | Upstream | Per OS Infectivity factor envelope protein, required for oral infection | 1.65e-18 (2.31e-15) |
G66615A | gp19 | Non-synonymous | Suspected virulence factor | 3.49e-15 (4.89e-12) |
C78978T | gp94 | Intergenic | Suspected virulence factor | 5.93e-17 (8.32e-14) |
G78991A | gp94 | Intergenic | Suspected virulence factor | 1.12e-16 (1.57e-13) |
C126118A | ODV-E56-2 | Upstream | Occlusion-derived virus envelope protein required for particle formation | 1.43e-18 (2.01e-15) |
A132593C | Helicase-2 | Non-synonymous | Unwinds DNA and is critical for DNA replication | 1.11e-09 (1.56e-06) |
T140117C | PIF-3 | Upstream | Per OS Infectivity factor envelope protein, required for oral infection | 2.96e-15 (4.15e-12) |
We isolated and sequenced DNA from these samples, quality filtered reads, mapped to the genome, and called genetic variation in the viral genomes to assess the extent of adaptation in each viral population. Consistent with an arms-race between host and virus, envelope and novel virulence (GrBNV-like) genes have a significantly higher proportion of substitutions fixed by adaptive evolution, compared to other viral genes (McDonald-Kreitman-based statistic Direction of Selection [Stoletzki and Eyre-Walker, 2011] and Selection Effect [Eilertson et al., 2012; Figure 4—figure supplement 1], DoS >0, GLM t-value >1.31, p-value<0.05).
Recurrent adaptive evolution in viral proteins known to interact with the host-immune system suggests an arms-race between D. innubila and DINV. To better understand the interactions between host and virus, we performed an association study between viral titer and natural genetic variation in both host and virus. We consider viral titer a proxy for virulence. For each virus-infected individual, we quantified viral titer (as the logarithm of the viral genome coverage normalized to host autosomal genome coverage) and performed an association study across both host and virus variable sites to identify variants significantly associated with viral titer using PLINK (Purcell et al., 2007).
Of 5283 viral SNPs in the 155kbp DiNV genome, 1,403 SNPs are segregating in at least five infected host individuals (five was our minor allele threshold for inclusion in the association study). Of those 1,403 SNPs, 78 are significantly associated with viral titer after multiple testing correction (FDR < 0.01, Significantly associated SNPs p<0.001 over 1000 permutations, Figure 1A). Of these, 16 are less than 2000 bp upstream of the start site of a gene, 18 are coding nonsynonymous, 11 are coding synonymous and 33 are intergenic. The SNP with the absolute lowest p-value for association with viral titer was also the only significantly associated SNP that seemed to segregate within individuals. In fact, the frequency of the derived allele in this nonsynonymous polymorphism in the active site of Helicase-2 shows a negative linear relationship with the log of viral titer (Figure 1—figure Supplement 2, GLM t-value = −20.516, p-value=5.55e-21). However, when we ranked samples by viral titer, we found the derived SNP frequency exhibited an almost bell-shaped curve when plotted against rank (Supplementary Figure 2).

Viral genome-wide association study for DiNV titer in wild D. innubila.
(A) Manhattan plot for each DiNV SNP and the significance of its association with DiNV titer. SNP point, shape and color denotes if they are upstream (black upward arrow), downstream (black downward arrow), intergenic (gray cross), synonymous (blue circle), non-synonymous (red square) or nonsense (green asterisk) mutations. Named SNPs are either part of the significantly associated viral haplotype or, if a smaller size and italicized, are in or near other genes of interest (e.g. Helicase-2). The FDR-corrected p-value cutoff of 0.01 is shown as a dashed line (multiple testing correction for 1403 tests), while the permutation-based genome-wide significance threshold of p=0.01 is shown as a dotted line (based on 1000 permutations). (B) Viral titer for individual wild-caught flies infected with Low and High DiNV haplotypes (containing all 11 High type alleles). The middle bar represents median value, upper and lower bars represent 25th and 75th percentile and whiskers represent a 95% confidence interval. (C) Association between the number of High type SNP variants and the viral titer of a sample. (D) Across five populations, the frequency of the High type is correlated with the frequency of the virus infection. (E) Linkage disequilibrium heatmap between the eleven focal haplotype SNPs pairwise, and with the Helicase-2 SNP. Tiles are colored by the estimated linkage between SNPs, from red (highly linked, r2 = 1) to white (unlinked, r2 = 0).
We also identified a striking association between viral titer and eleven strongly linked polymorphisms found across the DiNV genome (Figure 1A & D, highlighted SNPs, Table 1, Figure 1—figure supplements 1–3, Sig. SNPs). These SNPs were significantly associated in every population when performing the association study on individual populations, or as a single group with population as a covariate. We binned strains with one of the two complete sets of alleles and referred to them as the ‘High Type’ and ‘Low Type’ (Figure 1B, only contains individuals infected by viruses with a complete set of ‘High’ or ‘Low’ Type SNP variants). This multilocus genotype includes three non-synonymous SNPs, five SNPs in the UTRs of known virulence factor genes and three intergenic SNPs (Table 1). Viral titer is, on average, 100-fold higher in individuals infected with the full High type virus compared to the ancestral Low type (Figure 1B). When comparing the ‘High’ and ‘Low’ viral type to Helicase-2 allele frequency, we found that viral titer increases with Helicase-2 SNP frequency in the Low type (and some intermediate types), while titer decreases with Helicase-2 SNP frequency for most intermediate and High types, suggesting some form of negative interaction between the High type SNP variants and the Helicase-2 SNP (Supplementary Figure 2). Consistent with this, we found a negative interaction between Helicase-2 SNP frequency and the presence of High type on viral titer (GLM Log10(titer) ~Helicase-2 SNP * Haplotype: t-value = −7.815, p-value=6.39e-14).
Since these 11 SNPs are almost perfectly linked (Figure 1E), if one SNP were a false positive result, all would be false positives. We sought to address this in two ways. First, these SNPs are significantly associated with viral titer in every population when performing the association test separately or as a group with population as a covariate. Therefore, the likelihood of the same SNP(s) being false positives in three independent tests is quite small (p<1E-09 over 3 iterations of 1000 permutations). Second, we permuted the titer among all individuals within populations and performed the association test (100,000 total permutations). For each permutation, we took the SNP with the lowest p-value from each permutation test and the 10 other SNPs most strongly linked to it. We then compared the viral titer in those with this permuted ‘High’ haplotype to those with the ‘Low’ haplotype. In no case were the permuted differences between High and Low haplotypes as large as those observed from the real data (Figure 1—figure supplement 2A, 100,000 permutations), and the distribution of true p-values diverges dramatically from the permuted null expectation (Figure 1—figure supplement 2B, black). We found this divergence from the null expectation of p-values disappears when including viral haplotype as a covariate (Figure 1—figure supplement 2B, red), suggesting that the High type drives most of the signal seen in the association study.
Though we found few strains with an intermediate number of SNPs (intermediate types), viral titer appears to increase as the number of High type SNPs increases (Figure 1C, GLM t-value = 34.971, p-value=5.912e-16), but the rate of increase slows as the number of High type SNPs increases suggesting diminishing returns (Figure 1C). Some of these polymorphisms are associated with known virulence factors, or are related to the formation of the viral envelope co-opting the host vesicle trafficking system and are rapidly evolving in nudiviruses (e.g. 19K, ODV-E56, PIF-3) (Rohrmann, 2013; Hill and Unckless, 2017; Hill and Unckless, 2018). Additionally, several are associated with genes exclusive to a few nudivirus genomes thought to be novel virulence factors, including gp83, a gene that downregulates Toll-induced antimicrobial peptides (AMPs) and upregulates those induced by IMD (Palmer et al., 2019). Both the Toll and IMD pathway may interact with DNA viruses (Zambon et al., 2005; Costa et al., 2009; Merkling and van Rij, 2013; Ferreira et al., 2014; Lamiable et al., 2016; Palmer et al., 2019).
Among populations, we found a positive correlation between the frequency of the High type and overall DiNV infection frequency (Figure 1D, GLM logistic regression z-value = 6.104, p-value=0.00883), suggesting that the High type may have a higher transmission rate than the Low type, resulting in a larger number of new individuals infected per DiNV-infected individual. We also found both viral types in collections from 2001 and 2017, with the High type significantly more common in the 2017 collection (Fisher Exact Test p-value=0.0167, Figure 1D).
It is possible that mutations could segregate within individuals, but this is relatively rare overall, and we found no evidence of the 11 haplotype SNPs segregating within any individual fly (Supplementary Figure 2). In fact, no individual fly contains more than two significantly associated SNPs segregating within the sample, suggesting that hosts are either infected completely with Low type or High type virus particles.
The high DiNV viral type more effectively suppresses the Drosophila immune system
Given the striking difference in viral titer between High and Low type viruses, we sought to further characterize the differences in infection dynamics between the types, focusing on the differences in expression between viral types. We sequenced mRNA from 80 wild D. innubila males collected in 2018 (Supplementary file 1 - Table 2, 40 infected with DiNV, 40 uninfected) and performed a differential expression analysis between infected and uninfected individuals. Few genes were differentially expressed (DE) between infection states in D. innubila, but these DE genes were enriched for several interesting categories. Specifically, we found IMD-induced antimicrobial peptides (AMPs) were upregulated upon DiNV infection, while one Toll-induced AMP, and several chorion and heat shock protein genes were downregulated (Figure 2—figure supplement 1A). We also compared these results to a laboratory experiment in D. melanogaster of differential expression after infection with a close relative of DiNV (Palmer et al., 2018). These same genes tend to also be differentially expressed in D. melanogaster. Of the 12 genes which are differentially expressed in the same direction in both species, five are AMPs and five are chorion genes (Figure 2—figure supplement 1B). This suggests that even with significant evolutionary divergence in both host and virus, the transcriptional response to infection is similar in both hosts. It was unexpected that chorion proteins are differentially expressed upon infection by DiNV and Kallithea virus (Figure 2, Figure 2—figure supplement 1), especially as DiNV is thought to infect the Drosophila gut (Unckless, 2011) and these samples are male. This suggests DiNV may infect more tissues than just the gut, and that chorion proteins are expressed in the germline of both D. innubila sexes, potentially fulfilling different roles in different sexes/species.

DIfferential expression between two viral haplotypes.
(A) Differential expression of D. innubila and DiNV genes between D. innubila infected with either the Low type or High type DiNV multilocus genotypes. For host genes, the log-fold change of mRNA fragments per million fragments is compared, while for viral genes the log-fold change of viral mRNA fragments per million fragments per viral particle is compared. Genes are colored/labelled by categories of interest, specifically antimicrobial peptides (AMPs), proteins involved in the extracellular matrix and viral proteins. Specific genes of interest, such as Myd88, are also named. The FDR-corrected significance cut-off of 0.01 (10,320 tests) is shown as a dashed line. (B) Expression (in FPKM per viral particle) of gp83 increases with the number of High type SNPs.
To determine how the High and Low type differ in their ability to infect D. innubila, we compared gene expression between D. innubila infected with the two types from this same collection (excluding intermediate types, which also showed a significant difference in viral titer between High type and Low type, p-value=0.0000143). We found 17 host genes and nine viral genes differentially expressed between types (after we controlled for virus copy number as FPKM/titer for viral genes Figure 2A, FDR-corrected p-value<0.01). Specifically, three Toll-regulated immune peptides (IM33, Bomanins BomBC2, and BomT2) and one JAK-STAT regulated immune peptide (Listericin) have reduced expression in High type infected individuals compared to the Low type (Figure 2, Figure 2—figure supplement 2). Viral genes of interest (PIF-3, 19K, gp83) have higher expression per viral particle (FPKM/titer) in the High type compared to the Low type. gp83 also increases in expression per viral particle (FPKM/titer) as the number of High type alleles increases (Figure 2B, t-value = 13.732, p-value=3.36e-15). Together these results suggest that the High type has increased expression of key virulence factors, which in turn, manipulate the expression of host genes involved in immune defense to result in the observed differences in viral titer. Specifically, higher gp83 expression may cause the lower Toll-mediated AMP expression (Palmer et al., 2019). We reasoned that if the High type is suppressing the host Toll pathway, Myd88, the Toll signaling protein upstream of AMPs would also have lower expression. We did find Myd88 expression is lower in strains infected with the High type (Figure 2A, though no significantly so), which in turn might prevent the host from enacting a proper immune response to DiNV infection (Figure 2, Figure 2—figure supplement 1).
Experimental infections recapitulate differences in viral type virulence
To assess if the virulence differs between virus types, we performed experimental infections of D. innubila males using viral filtrate of strains infected with one of the two types of DiNV. Before injections, we used qPCR to identify the viral titer of each sample and dilute all samples to the same viral concentration. In these experiments, as viral titer increases, the survival of infected flies decreases, regardless of viral type (Figure 3—figure supplements 1 and 2, ANOVA residual deviance = 3.536, p-value=2.454e-07, Cox Hazard Ratio z-value >2.227, p-value<0.02592). In both types viral titer also increased for the first 3 days of infection (GLM t-value = 9.817, p-value=3.6e-14). This established that viral titer is a reasonable proxy for virulence and that DiNV is virulent in D. innubila.
For a set of four High type and four Low type-infected individuals, we isolated viruses, diluted to roughly equal concentrations of viral particles and performed infections for replicates of 10 males with microneedles dipped in one of the filtrate samples. Survival is significantly lower for flies infected with High type viruses when compared to either flies pricked with sterile media (Figure 3A, Cox Hazard Ratio z-value = 3.671, p-value=0.000242) or those pricked with Low type virus (Figure 3A, Cox Hazard Ratio z-value = 4.611, p-value=4.02e-06). Flies pricked with Low type virus show a non-significant reduction in survival compared to control flies (Cox Hazard Ratio z-value = 1.353, p-value=0.176). We also measured viral titer over time using qPCR, and found titer increases through time in flies infected with either type (Figure 3B, GLM Log10(titer) ~days + type + vial | strain, days t-value = 9.912, p-value=1.76e-14). Flies infected with High type virus have significantly higher viral titer compared to flies infected with Low type virus (Figure 3B, GLM Log10(titer) ~days + type + vial | strain, type t-value = 3.934, p-value=0.000211). These results suggest that the higher viral titer observed in the High Type is also associated with higher virulence.

Effect of viral type in experimental infections.
(A) Survival curves of D. innubila infected with high and low viral types compared to control flies pricked with sterile media, for 15 days post-infection. Survival 5 days post-infection separated by strain is shown in Figure 3—figure supplement 2. (B) qPCR copy number of viral p47 relative to tpi in D. innubila infected with DiNV filtrate of high and low types.
DiNV types are under strong selection in the host
We next tested whether genes likely to be involved in host/virus interaction show signs of recurrent natural selection. Using McDonald-Kreitman based statistics for estimating the proportion of substitutions fixed by selection (McDonald and Kreitman, 1991; Stoletzki and Eyre-Walker, 2011; Eilertson et al., 2012), we tested whether genes that are associated with the High and Low types exhibited different signatures of natural selection compared to other viral genes. We calculated the Selection Effect, the proportion of substitutions fixed by adaptive evolution, weighted by the total number of substitutions in the genome (Eilertson et al., 2012). We found that genes associated with SNPs in the initial association study for viral titer, which defined the High and Low types (listed in Table 1 and including 19K, PIF-3 and LEF-4) have significantly higher rate of substitutions being fixed due to selection than background genes across all categories (Figure 4¸ type-associated genes versus all other, t-value = 2.718, p-value=0.00068). When separating the polymorphisms from the substitutions, we found that envelope genes have a significant excess of functional substitutions per site compared to other genes (Figure 4—figure supplement 2A, GLM t-value = 3.62, p-value=0.00107), while genes of unknown function have a significant excess of non-synonymous polymorphisms (Figure 4—figure supplement 2A, GLM t-value = 2.33, p-value=0.02241) and a deficit of non-synonymous substitutions (Figure 4—figure supplement 2A, GLM t-value = −3.894, p-value=9.93e-05), which will lower the average selection effect of background genes, suggesting that the relative excess adaptation seen in some genes is in part due to an excess of variation in genes on unknown function.

Genes implicated in host/virus interaction are rapidly evolving by positive selection in the Chiricahua population.
Difference in selection effect for viral and host gene categories of interest from nearby background genes (average shown as 0, the dashed line), as indicated by the proportion of substitutions fixed by adaptation, weighted by mutations in SnIPRE (Eilertson et al., 2012). Genes that have associated SNPs from the association study are highlighted in red, while genes which are differentially expressed upon infection, or between viral types are labelled in blue. All association study hits are also differentially expressed and labelled in red. Genes of interest are named.
We also performed an association study using the host polymorphism and found 13 significantly associated SNPs, after controlling for the viral haplotype (Figure 4—figure supplement 3, p<0.01), but found no significant enrichments (p-value>0.05) or genes of interest (e.g. those involved in Toll signaling or antiviral pathways). We therefore looked for enrichments in genes associated with the top 100 significantly associated SNPs versus all other genes and found piRNA genes enriched (GO enrichment = 3.44, p-value=0.035, FDR-corrected). Given that siRNA genes are not highly ubiquitously expressed in D. innubila (Hill et al., 2019), that DiNV reduces host fecundity (Unckless, 2011), and that a close relative of DiNV infects the host ovaries (Palmer et al., 2018), DiNV could interact with chorion proteins during oogenesis, and piRNAs could be suppressing DiNV (Lewis et al., 2018). Consistent with the arms race model, host genes we suspect are interacting with DiNV (such as the association study hits, AMPs, chorion genes, piRNA genes, and extracellular genes) show elevated levels of substitutions fixed by selection compared to background genes in D. innubila (Figure 4 and Figure 4—figure supplement 1, GLM p-value<0.05) (Hill and Unckless, 2020). Finally, differentially expressed chorion genes, extracellular genes and AMPs have significantly more adaptive substitutions than non-differentially expressed genes in the same categories (Figure 4, blue dots, differentially expressed versus all other T-test: D. innubila t-value = 4.755, p-value=0.000671). When separating the polymorphisms from the substitutions, we found an excess of functional substitutions per site in AMPs (Figure 4—figure supplement 2B, GLM t-value = 4.776, p-value=1.81e-06), driving their excess of adaptive substitutions compared to other genes. Overall, these results suggest strong selection is acting on both the host to suppress viral activity and the virus to escape this suppression.
Recombination in DiNV may facilitate the evolution of the high type
We were interested in examining the genetic variation of DiNV to determine the relationship of the host and the two putative types of DiNV. To determine the appropriate approaches for measuring these patterns, we first need to determine the effective rate of recombination in the virus. Though recombination is necessary for proper nudivirus replication (Kelly, 1982; Kamita et al., 2003; Rohrmann, 2013), often these recombination events will be between nearly identical viral particles resulting in no detectable signature of recombination across the genome (Rohrmann, 2013). For recombination to leave its signature in genetic variation, two divergent viruses must coinfect the same cell – we refer to this as effective recombination. It is unclear how common such recombination is in nudiviruses.
We used three methods to examine rates of effective recombination across the DiNV genome: First we used GARD to identify the number of recombination breakpoints across samples (Kosakovsky Pond et al., 2006). Second, we screened for recombination events by finding all four combinations of alleles between two SNPs. Third we calculated the linkage disequilibrium pairwise between all SNPs (Shin et al., 2006). We found recombination is relatively common in DiNV (Figure 1—figure supplement 3), with 307 potential recombination events genome-wide in recent history (based on combinations of alleles across our samples), and that the genome has relatively low linkage disequilibrium (Figure 1—figure supplement 3, Figure 5—figure supplements 1 and 2). In DiNV, the eleven SNPs significantly associated with viral titer (the High and Low types), are spread across the ~155 kilobase pair genome, yet are nearly perfectly linked to each other but not to other SNPs (Figure 1—figure supplement 3, Figure 5—figure supplements 1 and 2).
We wanted to examine if it was possible that the High type was generated by recombination of SNPs onto the same background, under the assumption that if the High type was formed via recombination we would expect to find recombination events each side of each significantly associated SNP. Using GARD and the four-allele test, six of the significantly associated SNPs show evidence of a recombination event on one side of the SNP, while three have evidence of a recombination event on each side of the significantly associated SNP (Figure 5—figure supplement 3, Supplementary Data). This suggests that recombination could have aided in the formation of the multilocus genotype, by recombination between two intermediate types to form the High type. However, we found no evidence of recombination between individuals in different populations (e.g. no recombinant haplotypes that are ½ CH and ½ PR), suggesting little movement of one or more significantly associated SNPs between populations. We also found no evidence of recombination between the complete High and complete Low types. This suggests that the initial SNPs for the High type have evolved on separate backgrounds in all three populations, though recombination between intermediate strains may have allowed for the formation of the complete High type, or for the generation of other intermediate types. For this to happen, each SNP may have recurrently evolved in each population, even if not sequentially.
The high viral type of DiNV evolved repeatedly in three D. innubila populations
We next sought to understand the evolutionary origin of the two putative types. Given that both types are found in all populations surveyed (Figure 1D) with no evidence of recombination between populations, we hypothesized that this could occur one of three ways: First, the derived haplotype was present ancestrally and has been maintained since before geographic isolation occurred. Second, the derived haplotype evolved following geographic isolation and has spread via migration between locations. Third, the derived haplotype has recurrently evolved in each location.
To distinguish between these possibilities and determine the timeframe of divergence, we used the site frequency spectrum of silent DiNV polymorphism to estimate effective population size backwards in time for all populations (Liu and Fu, 2015). We found that the three populations (CH, HU and SR) expanded from a single viral particle (Ne = 1) to millions of particles during the last glacial maximum (30–100 thousand years ago) when D. innubila settled its current range (Figure 5—figure supplement 4; Hill and Unckless, 2020). This supports a single invasion event during a host-range change for each location. PR appears to expand between 1 and 10 thousand years ago, suggesting a much more recent bottleneck during the range expansion in PR (Figure 5—figure supplement 4; Hill and Unckless, 2020).
We aligned genomic regions containing SNPs to two related nudiviruses, Kallithea virus and Oryctes rhinoceros Nudivirus (OrNV) (Wang et al., 2008; Hill and Unckless, 2018; Palmer et al., 2018). The High type alleles are not present in either Kallithea or OrNV, and are not found in short read information for wild D. melanogaster infected with Kallithea virus (Webster et al., 2015), suggesting they are derived in DiNV.
We generated consensus DiNV sequences for each infected D. innubila individual and created a whole-genome phylogeny to infer geographic diffusion of samples using BEAST2 (Bouckaert et al., 2014). We then performed ancestral reconstruction of the presence of the High type across the phylogeny using APE (Paradis et al., 2004). Our samples grouped as three populations (with HU and SR forming one population) and, consistent with our expectation, the Low type is the ancestral state (Figure 5A). Interestingly we found PR is nested within HU/SR, based on five SNPs which segregate in the Low type HU/SR but are fixed in all PR samples (Figure 5A), consistent with the expansion of DiNV north over time. Surprisingly, the High type appears to have evolved repeatedly and convergently within each population, forming separate groups within each population (Figure 5A). The High type also clusters within each population in a phylogeny or principal component analysis of all viral SNPs (Figure 5—figure supplement 1C), and when repeating these analyses while excluding the eleven focal SNPs.

The evolution and maintenance of two viral types.
(A) Phylogeographic reconstruction of the spread of DiNV through D. innubila, rooted on the Kallithea virus reference sequence, including a reconstruction of the High type evolution (with strains containing all 11 High type variants shown in purple, strains with an intermediate number of high type variants are shown in pink, and strains with no high type variants are shown in black). Branches are colored when the SNPs found in the background for each High haplotype are present in the population, showing that the background differs per population. Black branches show states where branch tips do not contain all the shared High/Low population specific background SNPs. (B) Order of mutations in the viral haplotype appearing in each population. Apart from three mutations the order is consistent between locations. SNP order has been bootstrapped across multiple sample phylogenies, SNPs with bootstrap support of >95% are labelled with a *. (C) The number of generations needed for ‘High titer’ mutations to evolve in simulated populations, given that each mutation increases the mutation rate. The number of generations between each mutation appearing decreases as titer increases.
We next surveyed each background SNP (e.g. SNPs not associated with the High or Low type) to determine if the general background supports one of the three outlined ways in which the High type evolved and spread in each location. We grouped SNPs by their presence in just the High type or Low type (supporting a single origin and spread by migration) or if they were unique to a single population but shared between the High and Low types. In total, 341 SNPs (24% of SNPs surveyed) are unique to a single population yet are still shared between both High and Low types (Figure 5 and Figure 5—figure supplement 1), compared to 23 SNPs (including the 11 High type SNPs) shared between locations but exclusive to High type samples. This is consistent with the lack of recombination between individuals in different populations, suggesting the High type associated SNPs are derived recurrently in each population on different backgrounds. Consistent with this, using TreeTime we found that the eleven significantly associated SNPs could have recurrently evolved across our samples, after accounting for recombination and mutation rate (Kosakovsky Pond et al., 2006; Sagulenko et al., 2018).
We identified 341 SNPs (161 for CH, 127 for HU and SR, and 53 for PR), are present in all High type samples of a single population but a variable proportion of Low types for that population (between 19 and 94%) and are unique to that population. This pattern fits with the High type recurrently evolving on a single background (a different background in each location), supporting recurrent evolution of the High type, and against two ancestrally maintained viral types. The population-specific background SNPs are spread throughout the DiNV genome, with little evidence of recombination with the High type SNPs, making it unlikely that these SNPs recombined onto different backgrounds, and were instead present when the High type first evolved in each population (Figure 5—figure supplements 1 and 2). We found little evidence of gene conversion producing this background signature (Figure 5—figure supplement 3). While we identified recombination events between background SNPs in the Low type, we found no evidence of this occurring due to their fixed states in the High type, likely because the High type has swept to higher frequencies within individuals before recombination can occur.
The elevated DiNV mutation rate associated with increased titer allows the high type to recurrently evolve in each population
Though there is strong linkage between the High type SNPs, they are not perfectly associated with each other (Figure 5—figure supplements 1 and 2). Using this slight disassociation and APE (Paradis et al., 2004), we performed ancestral reconstruction of SNP origins in each population (assuming recurrent evolution) and found that, excluding three variable SNPs, the evolution of these SNPs was a similar order in each population (Figure 5B).
The recurrent evolution of these mutations in a specific order is feasible as the mutation rate is high and the waiting time between mutations should decrease as titer increases. Logically, if the wait time for the appearance of a beneficial mutation is 1/(2Ne*μ*s), increasing Ne (with titer) should decrease the wait time for the appearance of the High type mutations (Gillespie, 2004). To determine if this recurrent evolution is plausible in our estimated timeframe (~10,000 years), we simulated viral populations using a discrete susceptible-infectious model implemented in deSolve (Soetaert et al., 2010) using estimated baculovirus mutation rates, recombination rates, ranges of viral titer taken our samples and estimated population sizes for each viral population (parameters described in the methods). In this model we used an effective mutation rate scaled to viral titer, considering the mutation rate per particle, so total mutations per generation increase with viral titer. For simplicity, we limited these simulations to the first five mutations as these are mutations of largest effect and appear to be necessary for the remaining six mutations to appear. We also included epistasis which reduces the increase in titer for each subsequent mutation (), as a similar reduction is seen in our samples (Figure 1C). The simulations suggest that waiting time for the first mutation that increases titer is highly variable between replicates but usually occurs within 1000 generations (~200 years at most, assuming the virus is transmitted from adult hosts to larval hosts via feces, with a 5 host generations per year as a conservative minimum estimate of viral generations per year, in >99.93% of replicates, Figure 5C). The average wait time for each subsequent mutation decreases monotonically (GLM t-value = -2.389, p-value = 0.03686). In most cases, the next mutation appears in the background of the previous high titer mutation (Figure 5C) due to the elevated effective mutation rate and increased basic reproduction number (R0). In 34 of 1000 simulations, when a mutation does appear on a different background, recombination facilitates the generation of the full complement of mutations. Given that effective mutation rate is much higher than the effective recombination rate following the appearance of the first mutation (as recombination does not scale with titer in nudiviruses and baculoviruses Kamita et al., 2003), this likely accounts for the mutation rates dominance in our simulations, though recombination likely does play a role in recombination between intermediate strains. The accumulation of mutations occurs at close to a geometric (approximately exponential) rate. Additionally, the standard deviation of time wait times also decreases with each new mutation (GLM t-value = -2.441, p-value = 0.04241), increasing the certainty that the entire multilocus genotype will appear in a population rapidly once the initial mutations appear. This chain reaction of adaptation could easily facilitate the repeated evolution of the virulent High type independently in three populations, with all eleven mutations fixing in a population within 6000 generations (~1200 years maximum) in all replicates (3372 generations on average, ~675 years maximum), a plausible amount of time given our estimated timeframe.
Both viral types are found in two other Drosophila species and have also evolved in a geographically distinct population
Since we found two types of DiNV are present in all D. innubila populations, and that other species are infected with DiNV (Unckless, 2011), we hypothesized that another species could be a reservoir for the less virulent Low type. We chose to study D. azteca from the Chiricahuas since it is frequently infected with DiNV (~33% infection), overlaps with D. innubila, and is genetically divergent (40–60 million years) which could mean a very different genetic interaction between host and virus (Unckless, 2011). We also examined DiNV-infected D. falleni (a close relative of D. innubila with nonoverlapping geographic range, collected in Athens, Georgia) as an outgroup. We sequenced 36 D. azteca and 56 D. falleni. Both viral types are present in both additional species, but the High type is rare in D. azteca (Figure 6B). The High type has a significantly higher titer than the Low type in both cases (Figure 6A,D. azteca GLM t-value = 6.71, p-value=0.0056, D. falleni GLM t-value = 8.12, p-value=0.000371). Viral titer is not significantly different across species for either High or Low type (Figure 6A, GLM t-value = −1.351, p-value=0.179). We also found the D. azteca samples cluster with CH D. innubila samples and contain the CH background SNPs (Figure 5—figure supplement 1C), suggesting no differentiation in the virus infecting different species. D. falleni DiNV, on the other hand, clusters completely separately from the other samples, likely due to its geographic separation. However, the D. falleni samples still have a derived cluster of High type virus, suggesting a fourth independent evolution of the High type in Georgia. Despite the lack of divergence between viruses infecting the two species in Arizona, a lower proportion of the D. azteca population is infected with DiNV, and the High Type is less common than the Low type DiNV (Figure 6B). Perhaps, even though the relative differences in titer are preserved between the two species, the Low Type is favored in D. azteca because this reduced virulence leads to a greater basic reproductive rate for the virus in D. azteca. Thus, the two types of the virus may be maintained in both host species because though they have become specialized to maximize fitness in one host, transmission between host species could lead to their continued presence in both hosts.

Titer and frequency of different DiNV types infecting different species.
(A) Viral titer for CH samples of D. azteca, D. falleni and D. innubila infected with High and Low type DiNV. The D.innubila data presented here is a reconstruction of Figure 1B. (B) Proportion of D. azteca 2017, D. falleni and D. innubila 2017 CH population infected with High and Low type DiNV.
We repeated our association study in D. azteca, D. falleni, and each D. innubila population separately. In all cases we found the 11 High type SNPs are associated with higher viral titer (GLM t-value >4.28, p-value>0.0001 in all cases). This was not the case for the Helicase-2 SNP (despite its presence in D. falleni) nor any other SNPs shared between populations. After controlling for the High type, we found no other significant DiNV SNPs in D. azteca associated with viral titer. For DiNV infecting D. falleni, we found 478 significantly associated SNPs (FDR-corrected p-value<0.01), though none of them with as large an effect as the High type associated SNPs.
Discussion
Viruses are constantly evolving not just to propagate within a host, but also to optimize their infection across hosts. This optimization involves the relationship between their ability to infect an individual and to transmit to others, tempered by the pathogenic effects on infected hosts caused by viral activity (May and Nowak, 1995; Lipsitch et al., 1996; Alizon and van Baalen, 2008). Since the host is also evolving in response to the virus, an evolutionary arms-race often ensues (Dawkins and Krebs, 1979; Kaltz and Shykoff, 1998; Daugherty and Malik, 2012). DNA viruses have large genomes and often recombination, placing them as a somewhat transitionary pathogen between RNA viruses, bacteria and eukaryotic pathogens and parasites. Here, to work towards expanding our understanding of the co-evolution of viruses and their hosts, we examine the population dynamics of Drosophila innubila Nudivirus (DiNV), a DNA virus infecting D. innubila (Unckless, 2011). Within our set of viral samples, we found two DiNV types which differ by 11 SNPs (Figure 1, named High and Low types). One haplotype (the High type) is associated with higher viral titer, likely due to an increased manipulation of the host-immune system and increased expression of viral factors (Figure 2).
The derived High type has likely recurrently evolved in each population since the last glacial maximum (~10,000 years ago). The two types appear to not co-infect individuals, and mutations appear in a similar order as if navigating an epistatic fitness landscape (Dobzhansky, 1937; Kondrashov et al., 2002; Gavrilets, 2004). Finally, despite the higher titer and transmission rate of the High type compared to the Low type, we found the two types in all sampled populations. Thus, a fundamental question is: why are there two haplotypes? Possible explanations are (a) that the two haplotypes are neutral and coexist due to genetic drift, (b) we caught High type amid a selective sweep, or (c) the two haplotypes are adaptively maintained due to a tradeoff.
Given the immense differences in titer and survival between High and Low types and the rest of the preponderance of evidence presented above, we find it implausible that the two haplotypes are associated with equal viral fitness and therefore evolving under a model of genetic drift (Gillespie, 2004).
It is possible we caught the High type DiNV in the middle of a selective sweep in each population. Differing basic reproduction numbers (R0) would result in changes in the ratio of types over time (as seen between 2001 and 2017, Figure 1D). As we only have two time points, we could be witnessing a selective sweep of the High type spreading to fixation (Nielsen, 2005), with recombination causing the observed differences in the background. As we found the High type appears to have evolved recurrently, it would be unlikely that we have caught intermediate sweeps in all four populations sampled (Figure 1D). Using the frequency of the High type between 2001 and 2017 CH samples, we can calculate the selection coefficient for the High type if increasing at an exponential rate (which is likely if mid-selective sweep given the intermediate frequency at both time points). D. innubila is active during the monsoon season if Arizona (late July to early September, 10 weeks with a generation time of 2–3 weeks) (Patterson and Stone, 1949). If we assume five generations per year (the minimum number of viral generations per year, limited by the host generation time) and an increase among infected individuals from 52% on 2001 to 71% in 2017, the selection coefficient = 0.0055, which suggests the High type would take ~350 years maximum to fix in a population once it has arisen (from a frequency 1/Ne to one in the viral population). Given the coalescence time of the two types is close to the expansion time of the two viruses (2–30 thousand years), this does not fit with our results, suggesting the two types are more likely to be adaptively maintained and not amid a selective sweep. Further, if the High type was sweeping, it would be remarkable for us to catch these sweeps occurring in all four populations sampled, given the lack of gene flow between populations, our estimated time to fixation of 350 years, and estimated initial infection ~30 thousand years ago.
Finally, two viral types could be maintained because of a tradeoff between the two virus types. We can imagine several scenarios for such a tradeoff. First, within a single population of genetically identical hosts, two viral types could coexist if they had equal basic reproduction numbers (R0) (Nowak and May, 1994; May and Nowak, 1995). The basic reproduction number in a simple epidemiology model is the ratio of the transmission rate of the pathogen to its virulence. So, while we’ve observed that virulence is higher for the High type virus compared to the Low type, it may be also true that transmission is also higher in the High type (inferred from Figure 1D) and those differences are perfectly balanced creating equal R0. We find this unlikely but a useful starting point for considering the adaptive maintenance of the two types. A more likely second scenario is one where in a single population of a given species, different host genotypes favor different viral type and that the interaction between optimal R0 among the viral genotypes and susceptibility of the host genotypes creates a stable equilibrium for both host genotypes and viral genotypes (Anderson and May, 1981). These host genotypes could also be different host species. In fact, in the Chiricahuas, D. innubila was more likely to be infected with the High type virus and D. azteca was more likely to be infected with the Low type virus (Figure 6) even though the two species are sympatric and the differences in viral titer between the High and Low types were qualitatively similar. Therefore, R0 might be higher for the High type in D. innubila and higher for the Low type in D. azteca with some transmission between host species. These community dynamics could maintain both viruses in both species. The second broad category of tradeoffs involves environmental heterogeneity, which could also create the conditions for the maintenance of both viral types. Furthermore, differing environmental conditions (such as differences in population densities) could affect the frequencies in each location (Anderson and May, 1981). Finally, different viral types could also preferentially infect different tissues, and have different titers due to limitations of the given infected tissues. Previous work in nudiviruses identified nudivurses that infect different tissues including the gut, fat body and ovaries (Jackson et al., 2005; Burand et al., 2012; Palmer et al., 2018). Two DiNV strains infecting different tissues may also help explain the lack of recombination between High and Low types, as viral particles cannot exchange information if they segregate to different tissues. However, this would ignore the fact that we found no evidence of coinfection and viruses infecting two tissues could more easily coinfect.
It is possible that the High type has not recurrently evolved in each population and instead the two viral types are maintained by selection in the face of ongoing gene conversion. However, this would not explain the shared background variants which are fixed in the High type and population exclusive, or the evidence of recurrent evolution in the assembled phylogeny (Figure 5). It is entirely possible however that within each population, selection is maintaining the High type due to the fitness deficits of intermediate types, while gene conversion erodes the association of the eleven haplotype SNPs, which is possible given the recombination events found around the significantly associated SNPs (Figure 3—figure supplement 1).
Studies of different viruses have found that they can undergo speciation-like events through the accumulation of genetic incompatibilities (Matsubara and Otsuji, 1978; Rokyta and Wichman, 2009; Meyer et al., 2016). In these cases, the two viral types cannot recombine and produce viable viral progeny. This may be occurring in DiNV if infection by one viral type limits the ability of the other type to infect the same cell – akin to prezygotic reproductive isolation (Coyne and Orr, 2004). Alternatively, coinfection could happen but generate inviable recombinant particles (Meyer et al., 2016) – akin to postzygotic reproductive isolation. For example, the fourth mutation might only be favored when in the presence of the derived allele of the third mutation because the third mutation causes a conformational change in protein A that allows the conformational change in protein B caused by the fourth mutation to produce infective viral particles (Orr, 1995). The fact that 19K is known to form a complex with other PIF proteins and other membrane proteins suggests it could play the central role in the formation of this incompatibility system (Wang et al., 2007; Rohrmann, 2013). If incompatibilities are present between types, they are likely in rapidly evolving vital nudivirus genes (Hill and Unckless, 2017; Hill and Unckless, 2018). As well as being rapidly evolving across baculoviruses and nudivurses, we found these genes also have high rates of adaptive evolution (Figure 4) and are associated with SNPs that distinguish the High and Low DiNV strains (Figure 1), further supporting the idea that the few genes are key to nudivirus replication are targets of host-suppression and locked in arms-races across the entire nudivirus/baculovirus phylogeny (Hill and Unckless, 2017).
DNA viruses such as DiNV have complicated replication cycles and large genomes. This makes them a sort of evolutionary intermediate between RNA viruses (small genomes, high mutation rates) and eukaryotes (large genomes, low mutation rates) and tangential to bacteria and archaea (intermediate genomes, low recombination rates) (Rohrmann, 2013). However, adaptation appears to occur through changes in a few key proteins (Hill and Unckless, 2017; Hill and Unckless, 2018). Here we found the evolution of two competing viral types that differ in variants near these few key genes, which cause the strains to differ in virulence and titer. Overall our results suggest that the high mutation rates and extremely high levels of selection can result in the repeated and convergent evolution of novel host-virus interactions. Additionally, we found that these host-virus interactions for large DNA viruses can be much more complicated than previous models suggest (Dolan et al., 2018; Feder et al., 2019).
Materials and methods
Fly collection, DNA isolation, and sequencing
Request a detailed protocolIn this study, we used previously collected and sequenced D. innubila (Hill and Unckless, 2020). Briefly we collected these flies across the four mountainous locations in Arizona between the 22nd of August and the 11th of September 2017. Specifically, we collected at the Southwest research station in the Chiricahua mountains (~5400 feet elevation, 31.871 latitude −109.237 longitude), Prescott National Forest (~7900 feet elevation, 34.540 latitude −112.469 longitude), Madera Canyon in the Santa Rita mountains (~4900 feet elevation, 31.729 latitude −110.881 longitude) and Miller Peak in the Huachuca mountains (~5900 feet elevation, 31.632 latitude −110.340 longitude). Baits consisted of store-bought white button mushrooms (Agaricus bisporus) placed in large piles about 30 cm in diameter, at least five baits per location. We used a sweep net to collect flies over the baits in either the early morning or late afternoon between one and three days after the bait was set. Flies were sorted by sex and species at the University of Arizona and were flash frozen at −80°C before being shipped on dry ice to the University of Kansas in Lawrence, KS. During these collections we also obtained D. azteca which we also sorted by species and sex and flash frozen. D. falleni were collected using a similar method in the Smoky Mountains (~6600 feet elevation) in Georgia in 2017 by Kelly Dyer, these flies were then sorted at the University of Georgia in Athens GA and shipped on dry ice to the University of Kansas in Lawrence, KS.
For collected CH D. innubila, D. falleni and D. azteca, we attempted to assess the frequency of DiNV infection using PCR, looking for amplification of the viral gene p47. Using primers from Unckless, 2011, P47F: 5′–TGAAACCAGAATGACATATATAACGC and P47R: 5′–TCGGTTTCTCAATTAACTTGATAGC. We used the following conditions: 95°C 30 s, 55°C 30 s, 72°C 60 s per cycle for 35 cycles. We chose the cutoff for viral infection based these PCR infection results. We compared the PCR positives with the fold coverage in the genome and found that fly samples with 10x coverage of the viral genome for 90% of the viral genome also had PCR positive results, so we chose this as our cutoff.
We sorted 343 D. innubila flies, 60 DiNV positive D. falleni and 40 DiNV positive D. azteca which we then homogenized and used to extract DNA using the Qiagen Gentra Puregene Tissue kit (USA Qiagen Inc, Germantown, MD, USA). We prepared a genomic DNA library of these 343 DNA samples using a modified version of the Nextera DNA library prep kit (~350 bp insert size, Illumina Inc, San Diego, CA, USA) meant to conserve reagents. We sequenced the D. innubila libraries on two lanes of an Illumina HiSeq 4000 run (150 bp paired-end) (Data to be deposited in the SRA). We sequenced the D. falleni and D. azteca libraries on a separate run of a lane of an Illumina HiSeq 4000 (150bp paired-end).
For 80 male Drosophila innubila collected in 2018 (indicated in Supplementary file 1 - Table 2), we split the sample homogenate in half, isolated DNA from half as described above and isolating RNA using the Direct-zol RNA Microprep protocol (R2061, ZymoResearch, Irvine, CA, USA). We then polyA-selected on these samples to isolate mRNA and prepared a cDNA library for each of these 80 RNA samples using a modified version of the Nextera TruSeq library prep kit meant to conserve reagents and sequenced these samples on a NovaSeq NS6K SP 100SE (100 bp single end). We also sequenced DNA for these samples, with DNA isolated and prepared as above, also sequenced on a NovaSeq NS6K SP 100SE (100 bp single end).
Sample filtering, mapping, and alignment
Request a detailed protocolFollowing sequencing, we removed primer and adapter sequences using cutadapt (Martin, 2011) and Scythe (Buffalo, 2018) and trimmed all data using Sickle (-t sanger -q 20 l 50) (Joshi and Fass, 2011). We masked the D. innubila reference genome (Hill et al., 2019), using D. innubila TE sequences and RepeatMasker (Smit and Hubley, 2008; Smit and Hubley, 2013). We then mapped short reads to the masked genome and the Drosophila innubila Nudivirus genome (DiNV) (Hill and Unckless, 2018) using BWA MEM (Li and Durbin, 2009) and sorted using SAMtools (Li et al., 2009). Following this we added read groups, marked and removed sequencing and optical duplicates, and realigned around indels in each mapped BAM file using GATK and Picard, 2020 (http://broadinstitute.github.io/picard; McKenna et al., 2010; DePristo et al., 2011). We considered lines to be infected with DiNV if at least 95% of the viral genome is covered to at least 10-fold coverage. This most significantly overlaps with our CH D. innubila PCR results (χ2 = 71.791, p-value=2.392e-17). Additionally, index-switching at a rate of 0.0001 from the highest titer sample to an uninfected sample would still leave the uninfected with ~5.2-fold coverage of the viral genome, and so still considered as uninfected.
We then filtered for low coverage and mis-identified species by removing individuals with low coverage of the D. innubila genome (less than 5-fold coverage for 80% of the non-repetitive genome), and individuals we suspected of being misidentified as D. innubila following collection. This left us with 318 D. innubila wild flies with at least 5-fold coverage across at least 80% of the euchromatic genome, of which 254 are infected with DiNV (Supplementary file 1 - Table 1). We also checked for read pairs which were split mapped between the DiNV genome and the D. innubila genome using SAMtools.
For D. falleni we used a previously generated D. innubila genome with D. falleni variants inserted (Hill et al., 2019). We masked the genome with Repeatmasker (Smit and Hubley, 2013) and mapped short reads to the masked genome, the repeat sequences and the DiNV genome using BWA MEM and SAMtools (Li and Durbin, 2009; Li et al., 2009). Then, as with D. innubila we filtered for low coverage and mis-identified species by removing individuals with low coverage (less than fivefold coverage for 80% of the non-repetitive genome) leaving us with 56 D. falleni samples infected with DiNV.
For D. azteca, we downloaded the genome from NCBI (Accession: GCA_005876895.1) which we then called repeats from with RepeatModeler (Smit and Hubley, 2008). We masked the genome with Repeatmasker (Smit and Hubley, 2013) and mapped short reads to the masked genome, the repeat sequences and the DiNV genome using BWA MEM and SAMtools (Li and Durbin, 2009; Li et al., 2009). As with D. innubila we then filtered for low coverage and mis-identified species by removing individuals with low coverage of the D. azteca genome (less than 5-fold coverage for 80% of the non-repetitive genome), which left us with 37 D. azteca samples infected with DiNV. We then called DiNV variation using LoFreq (Wilm et al., 2012).
Calling nucleotide polymorphisms across the population samples
Request a detailed protocolFor the 318 sequenced samples with reasonable coverage, for host polymorphism, we used the previously generated multiple strain VCF file, generated using a standard GATK HaplotypeCaller/BCFTools pipeline. We used LoFreq (Wilm et al., 2012) to call polymorphic viral SNPs within each of the 254 DiNV-infected samples, following filtering using BCFtools to remove sites below a quality of 950 and a frequency less than 5%. We then merged each VCF to create a multiple strain VCF file, containing 5,283 SNPs in the DiNV genome. The LoFreq VCF (Wilm et al., 2012) output contains estimates of the frequency of each SNP in DiNV in each sample, to confirm these frequencies, in SAMtools (Li et al., 2009) we generated mPileups for each sample and for SNPs of interest (related to viral titer), we counted the number of each nucleotide to confirm the estimated frequencies of these nucleotides at each position in each sample. To confirm that there are no coinfections of types, we also subsampled samples and randomly merged Low and High type samples and again generated mPileup files, for SNPs of interest we again counted the number of each nucleotide at each position and confirmed these matched our expected counts in the merged files. We then compared these artificial coinfections to actual samples to confirm the presence or absence of coinfections, finding no samples consistent with coinfections. We then used SNPeff to identify the annotation of each SNP and label synonymous and non-synonymous (Cingolani et al., 2012). We extracted the synonymous site frequency spectrum to estimate the effective population size backwards in time using StairwayPlot (Liu and Fu, 2015).
Using the viral VCF and the genetics r package we calculated r2 measure of linkage disequilibrium. We visualized the linkage between the eleven focal SNPs, and between the eleven focal SNPs and 1000 random SNPs using LDheatmap (Shin et al., 2006). We also sorted r2 scores by the types of SNPs the measure is between, if it is between the focal type SNPs, the SNPs found on each populations High type background, and other SNPs.
For 100,000 permutations, we randomized the viral titer associated with each strain. For each permutation, we binned individuals into high and low artificial haplotypes based on the allele for the SNP at 126118 (the most significantly associated SNP) and the 10 SNP most strongly linked to this SNP. Finally, we found the difference between the two artificial bins and compared this to the difference between the viral titer of the two haplotypes. We then counted the proportion of the 100,000 permutations with a difference lower than the true difference.
Identifying signatures of adaptive evolution using McDonald-Kreitman-based tests
Request a detailed protocolWe filtered the total viral VCF with annotations by SNPeff and retained only non-synonymous (replacement) or synonymous (silent) SNPs. We also mapped reads from Kallithea virus (dS ~ 0.133) and Orcytes rhinoceros Nudivirus (OrNV, dS ~ 0.279) to the DiNV genome to identify the ancestral state at each site and polarize SNPs to specific branches (Hill and Unckless, 2018). Specifically, we sought to determine sites which are derived polymorphic in DiNV and which are substitutions fixed on the DiNV branch of the phylogeny. After removing singletons, we used the raw counts of fixed and polymorphic silent and replacement sites per gene to estimate McDonald-Kreitman-based statistics, specifically direction of selection (DoS) (McDonald and Kreitman, 1991; Smith and Eyre-Walker, 2002; Stoletzki and Eyre-Walker, 2011). We also used these values in SnIPRE (Eilertson et al., 2012), which reframes McDonald-Kreitman based statistics as a linear model, taking into account the total number of non-synonymous and synonymous mutations occurring in user defined categories to predict the expected number of these substitutions and calculate a selection effect relative to the observed and expected number of mutations (Eilertson et al., 2012). We calculated the SnIPRE selection effect for each gene using the total number of mutations on the chromosome of the focal gene.
For the host, we repeated this process using the SNPeff annotated VCF in the SnIPRE pipeline to identify signatures of adaptation in the host genome.
We then calculated the difference in each statistic between each gene and the median of all other viral genes, to identify how much that gene deferred from the genomic background. We repeated this for the host genome, limiting the comparison to nearby genes (within 100kbp on the same chromosome).
Identifying differentially expressed genes between DiNV-infected and uninfected Drosophila innubila
Request a detailed protocolFor 100 male Drosophila innubila collected in 2018 (indicated in Supplementary file 1 - Table 2), we homogenized each fly separately in 100 µL of PBS. We then split the sample homogenate in half, isolated DNA from half as described above and isolating RNA using the Direct-zol RNA Microprep protocol (R2061). Using the isolated DNA, we tested each sample for DiNV using PCR for P47 as described previously, using 40 DiNV-infected samples and 40 uninfected samples. We then prepared a cDNA library for each of these 80 RNA samples using a modified version of the Nextera TruSeq library prep kit meant to conserve reagents and sequenced these samples on a NovaSeq NS6K SP 100SE (100 bp single end). We also sequenced DNA for these samples, with DNA isolated and prepared as above, also sequenced on a NovaSeq NS6K SP 100SE (100 bp single end). The DNA sequenced here was mapped as described above, with variation called as described above for other DNA samples.
Following trimming and filtering the data as described in the methods, we mapped all mRNA sequencing data to a database of rRNA (Quast et al., 2013) to remove rRNA contaminants. Then we mapped the short read data to the masked D. innubila genome and DiNV genome using GSNAP (-N 1 -o sam) (Wu and Nacu, 2010). We estimated counts of reads uniquely mapped to D. innubila or DiNV genes using HTSEQ (Anders et al., 2015) for each sample. Using EdgeR (Robinson et al., 2010) we calculated the counts per million (CPM) of each gene in each sample and counted the number of samples with CPM > 1 for each gene. We find that over 70.3% of genes have a CPM > 1 in at least 70 samples. For the remaining genes, we find these genes are expressed in all samples of a subset of the strains (e.g. DiNV uninfected, DiNV-infected, DiNV-high infected, DiNV-low infected). This supports the validity of the annotation of D. innubila, given most genes are expressed in some manner, and suggests our RNA sequencing samples show expression results consistent with the original annotation of the D. innubila genome.
We attempted to improve the annotation of the D. innubila genome to find genes expressed only under infection. We extracted reads that mapped to unannotated portions of the genome and combined these for uninfected samples, samples infected with High type DiNV and samples infected with Low type DiNV as three separate samples. We then generated a de novo assembly for each of these three groups using Trinity and Velvet (Schulz et al., 2012; Haas et al., 2013). We then remapped these assemblies to the genome to identify other transcripts and found the consensus of these two for each sample. Using the Cufflinks pipeline (Ghosh and Chan, 2016), we mapped reads to the D. innubila genome and counted the number of reads mapping to each of these putative novel transcript regions, identifying 15,676 regions of at least 100 bp, with at least one read mapping in at least one sample. Of these, 717 putative genic regions have at least 1 CPM in all 80 samples, or in all samples of one group (DiNV uninfected, DiNV-infected, DiNV-low infected, DiNV-high infected). We next attempted to identify if any of these genes are differentially expressed between types, specifically between uninfected strains and DiNV-infected strains, and between Low-type infected and High-type infected strains. Using a matrix of CPM for each putative transcript region in each sample, we calculated the extent of differential expression between each type using EdgeR (Robinson et al., 2010), after removing regions that are under expressed, normalizing data and estimating the dispersion of expression. We find that 26 putative genes are differentially expressed between infected and uninfected types, and 69 putative genes are differentially expressed between High and Low types. We took these regions and identified any homology to D. virilis transcripts using blastn (Altschul et al., 1990). We find annotations for 37 putative genes are either expressed in all samples, or differentially expressed between samples. Of the 14 putative genes expressed in all samples, nine have the closest blast hit to an rRNA gene, and five have hits to unknown genes. For 23 differentially expressed putative genes with blast hits, three genes are like antimicrobial peptides (IM1, IM14, IM3), these genes are significantly downregulated upon infection, like other Toll-regulated AMPs, and have significantly lower expression in strains infected with High type DiNV compared to Low types. The remaining 20 genes all have similarity to genes associated with cell cycle regulation, actin regulation and tumor suppression genes.
Identifying genes associated with viral titer in Drosophila innubila
Request a detailed protocolAs the logarithm of viral titer was normally distributed (Shapiro-Wilk test W = 0.05413, p-value=0.342), we used PLINK (Purcell et al., 2007) to associate nucleotide polymorphism to logarithm of viral titer in infected samples. We first generated a relationship matrix for the viral samples using PLINK. This kinship matrix allows us control for pairs of individuals who, on average, have similar alleles and which, on average, have dissimilar values. We then fit a linear model in PLINK including population, sex, Wolbachia presence, the date of collection and the distance matrix for relationship of each sample (inferred using PLINK, shown as relationship[strain]).
Before performing the association study, we removed host SNPs found in fewer than five samples and merged SNPs in perfect linkage within 10kbp of each other, filtering down from 5283 viral polymorphisms to 1403 viral polymorphisms. We then identified associations between the logarithm of viral titer and the frequency of the viral polymorphism in each individual sample, resulting in the following model:
Where hs = host sex, p=population of collection, w = Wolbachia presence, dc = date collected, relationship[strain]=distance matrix. Following model fitting, performed a posthoc analysis to determine if any clumped SNPs were significant and if they should be assessed separately, we also identified covariates which seemed to show little or no effect on viral titer (p-value>0.1) using an ANOVA in R (R Development Core Team, 2013), and removed these, refitting the model. This was done step-wise, resulting in the following model:
Following this we also performed an association study using PLINK (Purcell et al., 2007) in the host, using previously called host variation, and considering viral haplotype as an additional covariate. Before we performed this analysis, we clumped SNPs that are strongly linked (r2 >0.99) within 10kbp of each other.
Whereas previous, and with vh = viral haplotype. Using GOrilla (Eden et al., 2009), we found no gene categories enriched in the significant SNPs (Figure 4—figure supplement 3).
We repeated this analysis for DiNV variants in D. azteca and D. falleni separately. We performed the association study twice, first using the original model, then including viral haplotype as an additional covariate.
Because we are performing multiple correlated tests, we next determined the genome-wide significance threshold for the association between a SNP and the viral titer by permutation. We permuted the viral titer information randomly across samples and repeated the genome-wide association studies as described above. We recorded the minimum p-value across the genome from 1000 permutations to generate a null distribution. To identify if the SNPs involved in the High type are more likely to appear together than expected, we again used a permutation test to identify the frequency that these combinations occur across 1000 permutations, where a P of 0.01 = 7.987812e-10.
Estimating viral titer using qPCR
Request a detailed protocolFollowing the identification of the viral haplotype associated with viral titer we sought to determine the effect of viral haplotypes in actual infections. For 20 samples with fly homogenate, we determine the viral titer and haplotype following filtration with a 0.22 µM filter.
We performed qPCR for the viral gene p47 (Forward 5-TCGTGCCGCTAAGCATATAG-3, Reverse 5-AAAGCTACATCTGTGCGAGG-3) on 1 µL of fly filtrate per sample and compared the estimated Cq values across three replicates to estimated viral copy number to confirm viral concentration (protocol: 2 min at 95°C, 40 cycles of 95°C for 30 s and 59°C for 20 s, followed by 2 min at 72°C). Following this we diluted samples to similar Cq values, relative to the sample with the highest Cq value. We confirmed this by repeating qPCR with p47 primers of 1 µL of each sample.
For each filtrate sample we performed infections on 30 D. innubila males 4–5 days following emergence using pricks with sterile needles dipped in viral filtrate. We recorded survival of each fly each day and removed dead flies. Finally, we took samples 1, 3, and 5 days post infection and measured viral copies of p47 relative to tpi at each time point.
Phylogeography of DiNV infection and the evolution of the different viral haplotypes
Request a detailed protocolFor each DiNV-infected D. innubila sample, we reconstructed the consensus DiNV genome infecting them using GATK AlternateReferenceMaker and the VCF generated for each strain (McKenna et al., 2010; DePristo et al., 2011). We used BEAST2 with the SkyLine package (Stadler et al., 2013) to build the phylogeny of DiNV genomes using 100 million iterations with a burn in of 1 million, sampling every 1000 trees (Bouckaert et al., 2014). After generating a set of phylogenies in BEAST2, we checked that convergence had occurred across the samples using Tracer and generated a final consensus phylogeny using TreeAnnotator (Bouckaert et al., 2014). For each High type SNP, we reconstructed the evolution of the SNP as character states across the BEAST2 phylogenies to infer the SNP appearance order. We used the all different rates (ARD) discrete model in APE (Paradis et al., 2004) on 1000 randomly sampled converged BEAST2 phylogenies to infer SNP states at each branch and calculated the bootstrap support for the order of appearance for each mutation in each population. To independently identify recurrent mutations across the DiNV phylogeny, we also used TreeTime (Sagulenko et al., 2018). Finally, we created a matrix of SNPs present in at least two viral samples and used this matrix in a principal component analysis in R (R Development Core Team, 2013), labeling each sample by their viral type in the PCA (based on the presence or absence of the significant 19K SNP).
Simulating the evolution of the high and low viral haplotypes
Request a detailed protocolWe sought to simulate the infection of DiNV in D. innubila when considering the evolution of a high titer viral haplotype, specifically if two viral types can be maintained against each other at stable frequencies, and if the high viral haplotype with five shared mutations could evolve recurrently in the given time period given realistic parameters. We used the R package DeSolve (Soetaert et al., 2010) to simulate infection dynamics in an SI model. We did not include the resistance class, under the assumption that flies won’t live long enough to clear the infection. Therefore, the proportion of population infected per generation is as follows:
Where = infection parameter (e.g. the increased likelihood an infected individual spreading its infection before dying) and = virulence parameter (e.g. the increased likelihood an infected individual has of dying before it can spread its infection). Infected = The proportion of the population infected with DiNV.
Susceptible = the proportion of the uninfected population.
We attempted to assess if the ‘high titer’ viral haplotype could possibly evolve recurrently in each population in the time scale seen in our findings. We again used the SI model, this time discrete with population sizes set to 1 million individuals, based on StairwayPlot estimates (Liu and Fu 2015), starting with 1 infected individual (1/Ne). We reasoned that if the wait time between beneficial mutations is 1/(2Ne*u*s), and increasing viral titer also increases Ne, then the increased titer also decreases the wait time between the appearance of mutations. For each infected individual in the population, we tracked the viral titer and also recorded the presence of absence of five mutations, with each mutation increasing the titer of infection, but each further mutation having successively smaller increases in viral titer (), representing the epistatic interaction of high titer associated mutations seen in the viral haplotype. Viral titer did not increase if a mutation occurs before the preceding mutation was present occurred, to matching the order of appearance of mutations seen in DiNV. In all simulations, we assumed 5 viral generations per year, based on the assumption of 5 host generations per year being the minimum number of viral generations (likely much higher), and so useful to conservatively estimate the amount of time required for the recurrent evolution of the High type.
In the simulations, we multiplied the infection parameter, mutation parameter and virulence parameter by viral titer, under the assumption that viral titer increases both infection and death rate, and the mutation rate is per viral particle (Maeda et al., 1993; Kamita et al., 2003). We considered a per site mutation rate of 10−6 mutations per generation, based on estimated baculovirus mutation rate (Rohrmann, 2013; Chateigner et al., 2015), with a specific mutation rate for each of the high titer variants of 3.3e-7 (1e-6 mutations per base per generation * 1/3 for the correct mutation) * viral titer, where each mutation occurs independently (so all five could arise in one generation at a rate of (3.3e-7)5 = 4.12e-33). We then simulated populations in replicate 1000 times for 100,000 generations with a starting infection frequency of 10% for the ‘low titer’ haplotype, recording the frequency of the virus in a population, the frequency of the haplotype and the time that each ‘high titer’ mutation reaches high enough frequency to escape stochastic behavior and behave deterministically under selection (Gillespie, 2004). We also factored in recombination between each variant site a rate of 2% per site-window per generation (Kamita et al., 2003). Specifically, we randomly paired viral genomes and in 2% of pairs we randomly recombined the variant combinations. We did not increase recombination rate as titer increased as this was not observed in nature (Kamita et al., 2003).
To estimate the possible selection coefficient for DiNV in the CH population, we assumed an exponential distribution and five host generations per year (80 host generations between 2001 and 2017). We then solved the following equation:
Where P2017 = the frequency of the High type among viral samples in 2017 (71%), P2001 = the frequency of the High type among viral samples in 2001 (52%), s = the selection coefficient and t = the number of generations (80). We then used this estimated selection coefficient in the same equation to find the number of generations (t) to go from 1/2Nes (0.0000714, assuming an Ne of 1000000) to fixation (0.99):
Experimental infections of Drosophila innubila with DiNV
Request a detailed protocolWe chose D. innubila samples infected with DiNV and with sequenced genomes, four infected with the High type DiNV and four infected with the Low type. For these samples we estimated their viral copy number per host genome as described previously. We used qPCR on p47 and tpi to find the differences in Cq values to calculate the concentration of each sample relative to the lowest concentration sample and diluted 50 µL of filtrate for each sample to match the concentration of each sample to the samples with the lowest titer (IPR07). For a separate 50 µL of the IPR01 sample, we performed 1 in 10 serial dilutions to give 45 µL of filtrate at full concentration, 1 in 10 concentration, 1 in 100 concentration and 1 in 1000 concentration. Using these sets of samples (matched titer and serial dilutions) we next performed experimental infections.
We transferred 50 D. innubila (of roughly equal sex ratio) to new food and let them lay eggs for 1 week, following this we collected male offspring aged 2–5 days for experimental infections.
Across four separate days in the mornings (between 9am and 11am), we infected the collected male flies with each sample. For flies in batches of 10, we performed pricks with microneedles dipped in the prepared viral filtrate. For each day we also had two control replicates of 10 flies pricked with microneedles dipped in sterile media. Following infections, we checked on each vial of 10 flies one-hour post infection and removed dead flies (likely killed by the needle instead of the virus). We also checked each vial each morning for 15 days, removing dead flies (freezing to determine the viral titer), and flipping flies to new food every 3–4 days. Checking at 10am each day, we recorded the day that each fly died, what filtrate they had been infected with, and what replicate/infection day set they belonged to. We next looked for differences in survival over time compared to sterile wound controls using a Fit proportional hazards regression model in R (R Development Core Team, 2013; Kassambara et al., 2017), considering titer, viral isolate (nested in haplotype) and replicate/vial as co-variates (day of death ~ [titer or haplotype | strain] + infection date).
For a second set of experimental infections (performed as described above, stabbed with diluted filtrate from different strains), we also removed three living flies 1 hr, 1 day, and 5 days post infection. Using qPCR, we found the difference in p47 log-Cq and tpi log-Cq to estimate the viral copy number for each sample over time.
Data availability
Sequencing data have been deposited on the NCBI SRA under the study accession: SRP187240 Genomes used in this study are available at the following accessions: Drosophila innubila - GCF_004354385.1 Drosophila innubila Nudivirus - GCF_004132165.1 Drosophila azteca - GCA_005876895.1.
-
NCBI Sequence Read ArchiveID SRP187240. Drosophila Sky Island sequencing.
-
Dryad Digital RepositoryDrosophila Sky Island data analysis.https://doi.org/10.5061/dryad.2fqz612mh
-
figshareID innubila_population_genomics/87662. Innubila population genomics.
References
-
Transmission-virulence trade-offs in vector-borne diseasesTheoretical Population Biology 74:6–15.https://doi.org/10.1016/j.tpb.2008.04.003
-
Basic local alignment search toolJournal of Molecular Biology 215:403–410.https://doi.org/10.1016/S0022-2836(05)80360-2
-
The population dynamics of microparasites and their invertebrate hostsPhilosophical Transactions of the Royal Society of London 291:451–524.https://doi.org/10.1098/rstb.1981.0005
-
BEAST 2: a software platform for bayesian evolutionary analysisPLOS Computational Biology 10:e1003537.https://doi.org/10.1371/journal.pcbi.1003537
-
SoftwareScytheScythe.
-
Viral suppressors of RNA silencingTrends in Plant Science 16:265–272.https://doi.org/10.1016/j.tplants.2011.02.010
-
SoftwareSpeciationSpeciation.
-
Evolution of RNA virus in spatially structured heterogeneous environmentsJournal of Evolutionary Biology 16:456–466.https://doi.org/10.1046/j.1420-9101.2003.00547.x
-
Rules of engagement: molecular insights from host-virus arms racesAnnual Review of Genetics 46:677–700.https://doi.org/10.1146/annurev-genet-110711-155522
-
How viruses hijack cell regulationTrends in Biochemical Sciences 36:159–169.https://doi.org/10.1016/j.tibs.2010.10.002
-
Arms races between and within speciesProceedings of the Royal Society of London B: Biological Sciences 205:489–511.https://doi.org/10.1098/rspb.1979.0081
-
SoftwareGenetics and the origin of speciesGenetics and the origin of species.
-
Mechanisms and concepts in RNA virus population dynamics and evolutionAnnual Review of Virology 5:69–92.https://doi.org/10.1146/annurev-virology-101416-041718
-
SnIPRE: selection inference using a poisson random effects modelPLOS Computational Biology 8:e1002806.https://doi.org/10.1371/journal.pcbi.1002806
-
Evolutionary dynamics in structured populations under strong population genetic forcesG3: Genes, Genomes, Genetics 9:3395–3407.https://doi.org/10.1534/g3.119.400605
-
Analysis of RNA-Seq data using TopHat and cufflinksMethods Mol Biol 1374:339–361.https://doi.org/10.1007/978-1-4939-3167-5_18
-
Viral evolution in deep time: lentiviruses and mammalsTrends in Genetics 28:89–100.https://doi.org/10.1016/j.tig.2011.11.003
-
SoftwarePopulation Genetics: A Concise GuidePopulation Genetics: A Concise Guide.
-
The genome of Drosophila innubila Reveals Lineage-Specific Patterns of Selection in Immune GenesMolecular Biology and Evolution 36:1405–1417.https://doi.org/10.1093/molbev/msz059
-
The dynamic evolution of Drosophila innubila NudivirusInfection, Genetics and Evolution 57:151–157.https://doi.org/10.1016/j.meegid.2017.11.013
-
Oryctes virus--time for a new look at a useful biocontrol agentJournal of Invertebrate Pathology 89:91–94.https://doi.org/10.1016/j.jip.2005.03.009
-
SoftwareSickle: A sliding window, adaptive, quality-based trimming tool for fastQ filesSickle: A sliding window, adaptive, quality-based trimming tool for fastQ files.
-
High-frequency homologous recombination between baculoviruses involves DNA replicationJournal of Virology 77:13053–13061.https://doi.org/10.1128/JVI.77.24.13053-13061.2003
-
SoftwareSurvminer: Drawing Survival Curves Using'ggplot2'Survminer: Drawing Survival Curves Using'ggplot2'.
-
Baculovirus ReplicationJournal of General Virology 63:1–13.https://doi.org/10.1099/0022-1317-63-1-1
-
GARD: a genetic algorithm for recombination detectionBioinformatics 22:3096–3098.https://doi.org/10.1093/bioinformatics/btl474
-
Pan-arthropod analysis reveals somatic piRNAs as an ancestral defence against transposable elementsNature Ecology & Evolution 2:174–181.https://doi.org/10.1038/s41559-017-0403-4
-
The sequence alignment/Map format and SAMtoolsBioinformatics 25:2078–2079.https://doi.org/10.1093/bioinformatics/btp352
-
Exploring population size changes using SNP frequency spectraNature Genetics 47:555–559.https://doi.org/10.1038/ng.3254
-
Coinfection and the evolution of parasite virulenceProceedings. Biological Sciences 261:209–215.https://doi.org/10.1098/rspb.1995.0138
-
Beyond RNAi: antiviral defense strategies in Drosophila and mosquitoJournal of Insect Physiology 59:159–170.https://doi.org/10.1016/j.jinsphys.2012.07.004
-
Molecular signatures of natural selectionAnnual Review of Genetics 39:197–218.https://doi.org/10.1146/annurev.genet.39.073003.112420
-
Superinfection and the evolution of parasite virulenceProceedings R. Soc 255:81–89.https://doi.org/10.1098/rspb.1994.0012
-
Population genetics of speciation: the evolution of hybrid incompatibilitiesGenetics 139:1805–1813.
-
Induction and suppression of NF- κB signalling by a DNA virus of DrosophilaJournal of Virology 93:e01443.https://doi.org/10.1128/JVI.01443-18
-
Standing genetic variation and the evolution of drug resistance in HIVPLOS Computational Biology 8:e1002527.https://doi.org/10.1371/journal.pcbi.1002527
-
Experimental evolution of resistance to an antimicrobial peptideProceedings of the Royal Society B: Biological Sciences 273:251–256.https://doi.org/10.1098/rspb.2005.3301
-
PLINK: a tool set for whole-genome association and population-based linkage analysesThe American Journal of Human Genetics 81:559–575.https://doi.org/10.1086/519795
-
The SILVA ribosomal RNA gene database project: improved data processing and web-based toolsNucleic Acids Research 41:D590–D596.https://doi.org/10.1093/nar/gks1219
-
SoftwareR: A Language and Environment for Statistical ComputingR Foundation for Statistical Computing, Vienna, Austria.
-
SoftwareBaculovirus Molecular BiologyBaculovirus Molecular Biology.
-
Genic incompatibilities in two hybrid bacteriophagesMolecular Biology and Evolution 26:2831–2839.https://doi.org/10.1093/molbev/msp199
-
Innate antiviral immunity in DrosophilaCurrent Opinion in Immunology 22:4–9.https://doi.org/10.1016/j.coi.2010.01.007
-
Dynamic evolution of the innate immune system in DrosophilaNature Genetics 39:1461–1468.https://doi.org/10.1038/ng.2007.60
-
TreeTime: maximum-likelihood phylodynamic analysisVirus Evolution 4:vex042.https://doi.org/10.1093/ve/vex042
-
LDheatmap : An R Function for Graphical Display of Pairwise Linkage Disequilibria Between Single Nucleotide PolymorphismsJournal of Statistical Software 16:jss.v016.c03.https://doi.org/10.18637/jss.v016.c03
-
SoftwareRepeatModeler Open-1.0RepeatModeler Open-1.0.
-
SoftwareRepeatMasker Open-4.0RepeatMasker Open-4.0.
-
Solving differential equations in R : package deSolveJournal of Statistical Software 33:1–24.https://doi.org/10.18637/jss.v033.i09
-
Estimation of the neutrality indexMolecular Biology and Evolution 28:63–70.https://doi.org/10.1093/molbev/msq249
-
The dawn of darwinian medicineThe Quarterly Review of Biology 66:1–22.https://doi.org/10.1086/417048
Decision letter
-
Diethard TautzSenior Editor; Max-Planck Institute for Evolutionary Biology, Germany
-
Dieter EbertReviewing Editor; University of Basel, Switzerland
-
Sebastian GagneuxReviewer; University of Basel, Switzerland
In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.
Acceptance summary:
This study reports on a phenotypic and genetic polymorphism of DiNVirus in natural populations of Drosophila. The authors present a series of experiments and assessments to understand how the polymorphism evolved and what implication it has for the host and conclude that the observed polymorphism arose multiple times independently and that it is maintained in a polymorphic state. The results are very clear and convincing and provide an excellent example for the power of natural selection in shaping host-parasite interactions.
Decision letter after peer review:
Thank you for submitting your article "Recurrent evolution of two competing haplotypes in an insect DNA virus" for consideration by eLife. Your article has been reviewed by three peer reviewers, one of whom is a member of our Board of Reviewing Editors, and the evaluation has been overseen by Diethard Tautz as the Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: Sebastian Gagneux (Reviewer #3).
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
As the editors have judged that your manuscript is of interest, but as described below that additional experiments are required before it is published, we would like to draw your attention to changes in our revision policy that we have made in response to COVID-19 (https://elifesciences.org/articles/57162). First, because many researchers have temporarily lost access to the labs, we will give authors as much time as they need to submit revised manuscripts. We are also offering, if you choose, to post the manuscript to bioRxiv (if it is not already there) along with this decision letter and a formal designation that the manuscript is "in revision at eLife". Please let us know if you would like to pursue this option. (If your work is more suitable for medRxiv, you will need to post the preprint yourself, as the mechanisms for us to do so are still in development.)
Your manuscript presents an exciting finding of the evolution and biology of a natural viral pathogen of Drosophila. You suggest that two distinct haplotypes of the virus evolved multiple times in isolated host populations. These two haplotypes are characterised by 11 SNPs and associated phenotypic traits (differences in viral titer). The independent evolution of such a complex trait seems rather unusual, making this nice example of evolution of host-parasite interactions.
The three reviewers of this manuscripts have a lot of praise for the study. However, they also raise a number of important points that needs to be addressed. Most important, the substantive technical points raised by reviewer 2 need careful consideration. This includes explaining and resolving the inconsistencies in the analysis and the right choice of the methods. Reviewer 1 pointed out problems with the simulation. Below are the detailed reviews. A final decision will only be possible when the problems in the analysis are solved.
Reviewer #1:
This study reports on a phenotypic and genetic polymorphism of DiNVirus in Drosophila. The study is very appealing as this study system is a natural system (unlike D. melanogaster) with a well understood ecology and biogeography. The authors go through a series of experiments and assessments to understand how the polymorphism evolved and what implication it has for the host. The main conclusion is that the observed polymorphism arose multiple times independently and that it is maintained in a polymorphic state. The mechanism for this maintenance is not entirely clear. A simulation model is used to bring some light into this puzzle. For the most part, the study is solid and well carried out. Below are a number of points that may help the authors to present their material more clearly.
The Introduction is not to the point. The Introduction moves among various aspects of host-parasite interactions without giving the reader an idea where this is going. At various places topics are raised that then later dismissed. For example: the first two paragraphs are about host response to viruses. In the third paragraph it become system specific, but now switches mainly to the DiNV system. It is not clear where one is going here. The Introduction (sixth paragraph) is summed up with a very generic phrase without much perspective on what is going to follow. At this place I still do not know what this paper will be about. A much more targeted Introduction is necessary. What are the questions? What was driving this research? What are the hypotheses?
With 11 SNPs spread across the chromosome and obligate recombination every round of replication, the proportion of viral offspring with sub-optimal multi-SNP-genotypes must be huge. How can the right genotypes be maintained? Discuss.
The authors run computer simulations (called SIR models, but they seem actually to be SI models) to support their ideas about virus evolution. The simulation regarding the accumulation of mutation is fine and gives quantitative support for presented evolutionary scenario.
A second simulation is about the competition of the two viral types. This is by far the weakest part of the manuscript. I am not convinced about the value of this simulation. The outcome can be predicted from the assumptions of the model, several of them are very speculative. Without data on transmission over the course of the infection, the simulations are not very helpful. I suggest to leave this out. It is ok to speculate about this in the Discussion, but it makes the Results part heavy and less strong. Also, the associated figure (Figure 7) is hard to understand.
The authors stress at multiple place that it is likely that the increased virulence of the high type is traded-off against transmission. The evidence for this is much weaker than the strength of these statements suggests. I strongly suggest to tone this down.
Reviewer #2:
This manuscript presents an exciting analysis of the evolution and biology of a natural viral pathogen of Drosophila. In outline: the authors claim that two distinct haplotypes of the virus co-occur in multiple host populations, that these 'evolved' independently (i.e. separate origins) through convergence, and that they have alternative 'life histories' (high titer and virulence, versus low titer and virulence), which may permit their coexistence.
The practical experiments and sequencing data are substantial and appear sound, and the work is likely to be of interest to a very broad host-pathogen audience. However, while some of these headline claims are well-supported, I have a number of serious concerns about the analyses and interpretation of others, and a few key methodological details are missing. The analyses would need substantial revision, or at least additional checks, before I would be convinced by the story.
1) Is there any possibility of circularity in defining the 'high' and 'low' haplotypes?
Two 'types' are defined based on 11 linked SNPs that are described as having 'high' or 'low' titer phenotypes. First, this requires a more robust approach to define 'significance' as (if there is any LD) tests are not independent, and sample sizes relative to predictive SNPs are small. Phenotypes (titers) should be permuted across genotypes (within populations) a thousand times and the analyses re-run for each permutation, then the tails of this distribution used to define significance, as commonly done for e.g. the DGRP. Second, this feels like it is dangerously circular: supposing all 11 SNPs were false positives, and arbitrary haplotypes erroneously defined based on them? Post hoc analysis of the difference between the two haplotypes could still show Figure 1B, because those differences were what drove the (erroneous) detection of the SNPs. So would Figure 1C, since that result necessarily follows from the definition of the haplotypes.
To convince me that the 'two haplotypes' interpretation is real, for each of the randomisation replicates one would need to define a 'high' and a 'low' haplotype based on the any randomisation spurious 'significant' SNPs and re-run the rest of Figure 1 to demonstrate that the separation between the real haplotypes is greater than that between those that result from permutation tests. I appreciate that this will require some computing time, but the authors already have all of the code, so it shouldn't be more than an afternoon of 'hands on' time.
Since the haplotypes overlap in phenotype space (Figure 1B) and not all have all SNPs (Figure 1C), how were intermediate haplotypes (with <11 SNPs) assigned to high or low?
2) Did the haplotypes evolve independently three or four times?
Of the 11 SNPs defining the haplotype, three are non-synonymous, five are in the UTRs of known virulence genes, and three are intergenic SNPs. Is it really credible that a specific base change has arisen and been selected for independently in each of four populations, at each of the 11 sites? i.e. always A->C being beneficial at a site, but never A->G at that site? Even for the non-coding ones? This is an extraordinary claim that would require extraordinary evidence, for which the rates of recombination and gene conversion seem at the heart of the argument.
In some places the authors appear to assume (or assert) that recombination is frequent, while in others they assume it is completely absent, and nowhere do they explicitly test for it. If it is common, then none of the tree-based analyses can be used. If it is absent, then it is very hard to explain the patterns of diversity or LD in Figure 1—figure supplement 3, and some of the simulations may be inappropriate.
The tree-based analyses seem to be the basis of the major claim that these haplotypes arose independently multiple times, and that the order of the mutations arising was similar. Any tree analysis absolutely requires the absence of recombination or gene conversion, and this needs to be explicitly tested for here (e.g. using GARD, or possibly if diversity is low by inferring the ARG using tsinfer). The failure of LD to decay with distance hints that recombination is absent, but gene conversion over very short distances could still shuffle mutations between haplotypes. However, assuming its branch lengths are in mutations, the shape of the tree in Figure 6A (many short tips below large crowns) shouts 'recombination' to me. Where did this tree come from? It is not ultrametric, so if it was inferred with BEAST this is not standard output.
If the authors do have evidence that recombination is absent, can they confirm that the other analyses (population size history from the SFS; simulations) also assumed zero recombination? Although zero recombination would make me worry even more about the p-values in the GWAS discovery of the SNPs that define the haplotypes – but permutation tests would help deal with that.
I am really confused about their view of recombination, because a couple of times they imply recombination may be common, for example by noting that recombination is required in Nudivirus replication – but then this would invalidate the tree-based analyses? Plus, mechanistic recombination is irrelevant without coinfection infection, since recombination only between identical haplotypes does not have any effect. The observation that no co-infections occur suggests recombination should be very rare – but that's not what the tree looks like. The sixth paragraph of the Discussion suggests that not enough thought has been given to the likelihood, or implications, of recombination.
Finally, if there is no recombination, then 'fully derived' haplotypes (those with all 11 SNPs) can only be as old as the youngest of the 11 SNPs, unless the SNPs arose multiple time within populations as well. Is that compatible with the patterns of diversity and the timescale? One check would be to mask those 11, re-infer the tree, and check that the clades are still monophyletic. If these arose by re-current mutations due to strong selection, they could be warping the tree – this might be akin to the problem of recurrent drug or MHC selected mutations in HIV trees, where the known sites are excluded before analysis.
If we believe the two haplotypes are real, a much more credible 'story' to me would be two distinct and potentially old haplotypes maintained by selection in the face of a low level of ongoing gene conversion (or recombination). Is there some aspect of the data that this does not fit? This seems just as good a 'story', just as interesting, and just as publishable.
3) The MK-like analyses
Although a relatively minor part of the story, the MK analysis is extremely unclear. In part this is because it doesn't appear in the Materials and methods.
i) What software was used? The supporting data implies VCFtools can do this (which surprised me) but the method it uses was not explained. Was SNIpre fitted with the original code? The text seem to imply that raw counts were used, which would not be suitable unless Ks was <0.3
ii) Although divergence is required, we're not told how (or from what species) it was estimated. If Ks>0.3, then something more than raw counts should certainly be used. If Ks>0.8, I would strongly advise against doing any sort of MK analysis.
iii) Figure 4 gives 'difference from background'. What is meant by difference from background? But if the background is a signal of strong constraint, could a positive signal here mean relaxed constraint rather than positive selection. What is the evidence that it's positive selection rather than relaxed constraint?
iv) In Figure 4—figure supplement 1 the populations are presented separately, and differences among them are interpreted as differences in positive selection. But surely the divergence number must be massive, larger than the polymorphism number, so almost all of this variation is due to differences in constraint affecting the Pn/Ps ratio. Or is the method one that uses high frequency derived SNPs as evidence of positive selection? If so, how was ancestral state identified? This is probably not possible to do reliably if Ks>0.3. We need to see some raw numbers, and a lot more detail on the methods.
Reviewer #3:
This is an interesting piece of work on the co-evolution of the Drosophilainnubila Nudivirus (DiNV) and its host. The authors analyzed several natural populations of D. innubila, some of which were partially infected with DiNV, and simultaneously characterized both the host and the infecting viral pathogen using a combination of DNA and RNA sequencing and various complementary analytical approaches. Using a GWAS approach, they discovered a viral variant with high virulence (High Type) that differed by 11 strongly linked SNPs from low virulence variants. They demonstrate that the High Type associates with a higher viral titer and increased host mortality, and also validate these findings using experimental infection assays. Using a transcriptomic approach, they show that the High Type overexpresses genes known to be linked to viral virulence, and this correlated with the under expression of host genes involved in antiviral immunity, indicating that the increased virulence of the High Type is at least partially due to the inhibition of host defense mechanisms. They further show that loci associated with differences in virulence were under strong selection for adaptation, particularly genes involved in the viral envelope and virulence proteins. Based on their reconstruction of the most likely evolutionary histories of the High Type in the different host populations, they further conclude that the High Type emerged multiple times independently, and yet, the High Type did not outcompete the Low virulence variants in these populations. This might indicate varying trade-offs between virulence and transmission in the High and Low Types that vary across these populations. Finally, they show that the same phenomenon for High Type evolution of DiNV can also be observed in other Drosophila species. I have just a few comments:
1) Throughout the manuscript, the authors switch back and forth between the present and past tense, which seems awkward from a stylistic point of, I'd suggest to stick to the past tense through-out.
2) The fact that the exact same 11 SNPs evolve multiple times independently in mostly the same order is interesting. The authors note that an expected alternative could be different SNPs emerging in the same genes, but they don't discuss the potential mechanism, by which the exact same SNPs seem to be preferred instead.
3) The authors found little evidence of mixed infection with both the High and Low Types and conclude potential in-compatibility. Please expand on the potential mechanisms of this.
4) Related to the above comment, the authors observe no "hybrids" of High and Low Types again, suggesting "incompatibility". Please discuss the difference between this genetic/genomic versus ecological incompatibility referred to above, as well as the potential link between these two types of incompatibilities.
5) Please rephrase the first sentence of the Discussion (what is the deference between "to better infect" and "optimizing the infection"?).
https://doi.org/10.7554/eLife.58931.sa1Author response
Reviewer #1:
[…] Below are a number of points that may help the authors to present their material more clearly.
The Introduction is not to the point. The Introduction moves among various aspects of host-parasite interactions without giving the reader an idea where this is going. At various places topics are raised that then later dismissed. For example: the first two paragraphs are about host response to viruses. In the third paragraph it become system specific, but now switches mainly to the DiNV system. It is not clear where one is going here. The Introduction (sixth paragraph) is summed up with a very generic phrase without much perspective on what is going to follow. At this place I still do not know what this paper will be about. A much more targeted Introduction is necessary. What are the questions? What was driving this research? What are the hypotheses?
We have rewritten the Introduction following your suggestions, thank you for the comments regarding this. We have moved sections around and rewritten parts so there is a better flow between paragraphs. We have also better provided our initial hypotheses and how these led to our findings. Specifically, we use the first paragraph to describe the host/virus coevolution and how viruses suppress the host, we then outline why people study viruses and why the DiNV system is ideal for viral study (with details of the system). Finally, we use the last paragraph to discuss the experiments we performed and a basic outline of our results.
With 11 SNPs spread across the chromosome and obligate recombination every round of replication, the proportion of viral offspring with sub-optimal multi-SNP-genotypes must be huge. How can the right genotypes be maintained? Discuss.
We have further expanded upon this in both the Results and the Discussion. In the Results we highlight that though there is obligate recombination, the chance that two strains with differing genotypes infect the same cell in an actual organism is very low, so we cannot know the actual rate or recombination, but it is likely lower than the high rate suggested and would lead to a lower proportion of intermediate strains. In the Discussion we also cover this and highlight that one possible difference between the High and Low strains is that they could infect different tissues (as is seen between other nudivirus strains). If the two types preferentially infect different tissues, then the chance of them co-infecting the same cell is impossible and would result in our observed absence of recombination between the High and Low type. We also discuss incompatibilities that could occur in suboptimal genotypes which could be the cause of their absence in our survey, or other factors which could drive the maintenance of the full haplotypes.
The authors run computer simulations (called SIR models, but they seem actually to be SI models) to support their ideas about virus evolution. The simulation regarding the accumulation of mutation is fine and gives quantitative support for presented evolutionary scenario.
A second simulation is about the competition of the two viral types. This is by far the weakest part of the manuscript. I am not convinced about the value of this simulation. The outcome can be predicted from the assumptions of the model, several of them are very speculative. Without data on transmission over the course of the infection, the simulations are not very helpful. I suggest to leave this out. It is ok to speculate about this in the Discussion, but it makes the Results part heavy and less strong. Also, the associated figure (Figure 7) is hard to understand.
We have now removed the second set of simulations and associated figures regarding the possibility of virulence/transmission equilibrium, as we agree it is the weakest part of the manuscript and is unlikely, as it is an unstable equilibrium. We now discuss the idea of a virulence/transmission equilibrium in the Discussion, however, we also discuss all the possible things that could have resulted in us finding two viral types (incomplete sweep, trade-off, optimized to different hosts).
The authors stress at multiple place that it is likely that the increased virulence of the high type is traded-off against transmission. The evidence for this is much weaker than the strength of these statements suggests. I strongly suggest to tone this down.
We have rewritten our Discussion to tone this down, in our Results we describe a relationship we find between increasing titer and increasing virulence. Beyond this, in the Discussion we examine multiple possible causes for the two maintained strains. We hypothesize that the two strains could be neutral and are drifting, that the High type is partway through a sweep, that a trade-off allows them to exist alongside each other or that each virus type is optimized to different host species/genotypes (with migration between two types). Overall, we think we have removed the majority of our Results and Discussion where talk about the trade-off, as it is not the only suggested explanation.
Reviewer #2:
[…] The practical experiments and sequencing data are substantial and appear sound, and the work is likely to be of interest to a very broad host-pathogen audience. However, while some of these headline claims are well-supported, I have a number of serious concerns about the analyses and interpretation of others, and a few key methodological details are missing. The analyses would need substantial revision, or at least additional checks, before I would be convinced by the story.
1) Is there any possibility of circularity in defining the 'high' and 'low' haplotypes?
Two 'types' are defined based on 11 linked SNPs that are described as having 'high' or 'low' titer phenotypes. First, this requires a more robust approach to define 'significance' as (if there is any LD) tests are not independent, and sample sizes relative to predictive SNPs are small. Phenotypes (titers) should be permuted across genotypes (within populations) a thousand times and the analyses re-run for each permutation, then the tails of this distribution used to define significance, as commonly done for e.g. the DGRP. Second, this feels like it is dangerously circular: supposing all 11 SNPs were false positives, and arbitrary haplotypes erroneously defined based on them? Post hoc analysis of the difference between the two haplotypes could still show Figure 1B, because those differences were what drove the (erroneous) detection of the SNPs. So would Figure 1C, since that result necessarily follows from the definition of the haplotypes.
There are two important considerations here so we will address them separately. We address the false negative problem below.
Defining significance since tests are not independent: Based on reviewer comments, we have significantly revised our approach at examining linkage disequilibrium and hope this clarifies some of the concerns. Given the signatures of LD, we did perform permutations and now include the permutation threshold in Figure 1A (slightly less stringent than the corrected P-value with a FDR of 0.01). So this was an important improvement but did not qualitatively change the result.
To convince me that the 'two haplotypes' interpretation is real, for each of the randomisation replicates one would need to define a 'high' and a 'low' haplotype based on the any randomisation spurious 'significant' SNPs and re-run the rest of Figure 1 to demonstrate that the separation between the real haplotypes is greater than that between those that result from permutation tests. I appreciate that this will require some computing time, but the authors already have all of the code, so it shouldn't be more than an afternoon of 'hands on' time.
The haplotype association problem: This is an important concern – if the first SNP is a false positive, they are all false positives. We think we have three convincing lines of evidence that these are not false positive SNPs. First, we have now performed the association study in the five populations (3 innubila, 1 azteca and 1 falleni) independently. Given the viral population structure among these populations (except maybe the sympatric innubila/azteca), these represent 5 independent tests and the haplotypes are significantly associated in all five populations. Thus a false negative would have to be a false positive 5 (or 4 if we don’t consider sympatric innubila/azteca as separate) times. Even if we consider our FDR threshold of 0.01 (these SNPs have lower P values), the likelihood of a false positive in 5 independent tests is 0.01^5 = 10^-10 or 1 in 10 billion (if 4 independent tests, that value drops to 10^-8 or 1 in 100 million). Second, we have followed the reviewer’s suggestions (we think) and permuted the titer associated with each strain 100,000 times, and each time binned strains by the most significant SNP and its 10 most closely linked SNPs, then find the difference in titer between two types. We find that the difference between the high and low type is significantly higher than by random chance for all SNPs, and no random combinations of alleles show the titer increase seen in Figure 1B. We have included the distribution of permuted differences between the artificial high and low types in Figure 1—figure supplement 2, with the true difference shown as a dotted line.
Since the haplotypes overlap in phenotype space (Figure 1B) and not all have all SNPs (Figure 1C), how were intermediate haplotypes (with <11 SNPs) assigned to high or low?
We have now changed the manuscript to consider intermediate strains to be a third subset of DiNV and have excluded them from analyses comparing High to Low type, including Figure 1B.
2) Did the haplotypes evolve independently three or four times?
As we have clarified in the manuscript, the haplotype has recurrently evolved three times to our knowledge in the Sky Islands, and at least once more in another location.
Of the 11 SNPs defining the haplotype, three are non-synonymous, five are in the UTRs of known virulence genes, and three are intergenic SNPs. Is it really credible that a specific base change has arisen and been selected for independently in each of four populations, at each of the 11 sites? i.e. always A->C being beneficial at a site, but never A->G at that site? Even for the non-coding ones? This is an extraordinary claim that would require extraordinary evidence, for which the rates of recombination and gene conversion seem at the heart of the argument.
We agree it is extraordinary that specific base changes always arising is strange, though this is not uncommon for viruses, especially given the low, but ubiquitous levels of codon usage bias seen in viruses. Additionally, it is possible that specific nucleotide changes are necessary for regulatory differences. We have discussed this in the manuscript and have also searched for gene conversion and find no evidence of it between the significant SNPs and surrounding SNPs (Figure 5—figure supplement 3). Our argument for increased effective mutation rate coupled to titer also increases the likelihood of this occurring, as we discuss in the manuscript.
In some places the authors appear to assume (or assert) that recombination is frequent, while in others they assume it is completely absent, and nowhere do they explicitly test for it. If it is common, then none of the tree-based analyses can be used. If it is absent, then it is very hard to explain the patterns of diversity or LD in Figure 1—figure supplement 3, and some of the simulations may be inappropriate.
We now address recombination in the second section of the manuscript (after the associations). The reviewer is correct that establishing and being clear about the recombination landscape is crucial for the rest of the manuscript. We have inserted also clarifications in the manuscript in the Discussion. In summary, while recombination events occur anywhere in the genome and frequently, the effective recombination rate is not as extremely high as we may have intimated, as most recombination events are between identical genomes. Recombination between two distinct haplotypes (including the low and high titer haplotypes) would require that the two types infect the same cell and this may be rare even if a host is superinfected. Even when recombination occurs between the two haplotypes, we hypothesize that recombinant genomes may have incompatibilities and so may not be able to leave the host cell. Based on these comments we have also changed our Discussion in the manuscript and have more thoroughly searched for evidence of recombination events between different SNPs. We find no evidence of recombination between populations, meaning each SNP must recurrently evolve in each population, however the final haplotype may have formed by bringing the SNPs onto the same background.
The tree-based analyses seem to be the basis of the major claim that these haplotypes arose independently multiple times, and that the order of the mutations arising was similar. Any tree analysis absolutely requires the absence of recombination or gene conversion, and this needs to be explicitly tested for here (e.g. using GARD, or possibly if diversity is low by inferring the ARG using tsinfer). The failure of LD to decay with distance hints that recombination is absent, but gene conversion over very short distances could still shuffle mutations between haplotypes. However, assuming its branch lengths are in mutations, the shape of the tree in Figure 6A (many short tips below large crowns) shouts 'recombination' to me. Where did this tree come from? It is not ultrametric, so if it was inferred with BEAST this is not standard output.
The reviewer is correct, we have now used GARD to identify recombination events and find evidence of recombination events between significantly associated SNPs and background SNPs. We also see evidence of recombination in our r2 analyses between background SNPs but not our significant SNPs. We have also regenerated the phylogeny in BEAST2 using BDsky/skyline (now in Figure 5, appropriate given the recombination findings). We have also generated a phylogeny in BEAST2 considering recombination and get a very similar phylogeny to that seen now in Figure 5. In all cases we find recurrent evolution of each high type SNP in each population, but the combinations of SNPs may either be due to recurrent mutation or recombination. As we find no evidence of recombination between populations, even if the full complement of 11 SNPs has not recurrently evolved on the same background, we still find evidence of recurrent recombination of each mutation in each population.
If the authors do have evidence that recombination is absent, can they confirm that the other analyses (population size history from the SFS; simulations) also assumed zero recombination? Although zero recombination would make me worry even more about the p-values in the GWAS discovery of the SNPs that define the haplotypes – but permutation tests would help deal with that.
To identify recombination events and measure associations of alleles, we have performed the following tests/used the following tools:
– GARD
– R2 in r
– 4 allele test between SNP combinations
We have also factored in recombination in the following analyses and still find support for recurrent evolution of each SNP (though SNPs may have recombined to form the final High Type):
– BEAST2
– TreeTimes
– DeSolve Simulations
– StairwayPlot
We find no evidence of recombination between populations, e.g. we find no recombinant haplotypes which are ½ CH and ½ PR ; however we find recombination within populations. The results shown in Figure 1—figure supplements 2 and 3, suggests recombination is on average quite common. Based on the GARD results and BEAST2 analysis, recombination could have occurred between intermediate types in each population, which may have generated the full High type. However, for this to occur each mutation must recurrently occur in each population on the intermediate type.
I am really confused about their view of recombination, because a couple of times they imply recombination may be common, for example by noting that recombination is required in Nudivirus replication – but then this would invalidate the tree-based analyses? Plus, mechanistic recombination is irrelevant without coinfection infection, since recombination only between identical haplotypes does not have any effect. The observation that no co-infections occur suggests recombination should be very rare – but that's not what the tree looks like. The sixth paragraph of the Discussion suggests that not enough thought has been given to the likelihood, or implications, of recombination.
We think the confusion is the difference between actual recombination (the physical process) and effective recombination (the genetic signature of recombination between two distinct haplotypes). Coinfection is required for effective recombination but not actual recombination. We expect recombination between neutral genotypes to be common if coinfections are common, but do not see recombination between the High and Low types, suggesting either they do not co-infect or recombinants between the two types are inviable. We have included this in the Discussion and reworked our discussion of recombination throughout.
Finally, if there is no recombination, then 'fully derived' haplotypes (those with all 11 SNPs) can only be as old as the youngest of the 11 SNPs, unless the SNPs arose multiple time within populations as well. Is that compatible with the patterns of diversity and the timescale? One check would be to mask those 11, re-infer the tree, and check that the clades are still monophyletic. If these arose by re-current mutations due to strong selection, they could be warping the tree – this might be akin to the problem of recurrent drug or MHC selected mutations in HIV trees, where the known sites are excluded before analysis.
If we believe the two haplotypes are real, a much more credible 'story' to me would be two distinct and potentially old haplotypes maintained by selection in the face of a low level of ongoing gene conversion (or recombination). Is there some aspect of the data that this does not fit? This seems just as good a 'story', just as interesting, and just as publishable.
In all cases, the phylogeny has been generated without the 11 focal SNPs but was also generated with the SNPs and is nearly identical in both cases, we have now mentioned this in the manuscript. The idea of the two old haplotypes being maintained in the face of gene conversion is interesting, and we now discuss this in the manuscript. However, the idea of long-term maintenance does not fit with any analyses we have performed, as High types do not cluster together. Additionally, we find the High type is present in the geographically distinct Drosophila falleni population and clusters completely separately from the innubila samples.
We do find some evidence of recombination events occurring around the significantly associated SNPs, so it is possible that multiple intermediate strains (with different complements of SNPs) have recurrently evolved in each population and have recombined to generate the full High type, or that the High type has evolved once in each population and is eroded as you described. We address this possibility in both the Results and the Discussion. We now have a dedicated section of our Results regarding detecting recombination and inferring if the High type SNPs were brought together by recombination. We are unable to conclude if the High type was formed exclusively by recombination (likely both sequential mutation and recombination play a role), but find no evidence of recombination between populations, which suggests that even if recombination brings the SNPs together onto the same background, the SNPs must have recurrently evolved in each population.
3) The MK-like analyses
Although a relatively minor part of the story, the MK analysis is extremely unclear. In part this is because it doesn't appear in the Materials and methods.
i) What software was used? The supporting data implies VCFtools can do this (which surprised me) but the method it uses was not explained. Was SNIpre fitted with the original code? The text seem to imply that raw counts were used, which would not be suitable unless Ks was <0.3
ii) Although divergence is required, we're not told how (or from what species) it was estimated. If Ks>0.3, then something more than raw counts should certainly be used. If Ks>0.8, I would strongly advise against doing any sort of MK analysis.
iii) Figure 4 gives 'difference from background'. What is meant by difference from background? But if the background is a signal of strong constraint, could a positive signal here mean relaxed constraint rather than positive selection. What is the evidence that it's positive selection rather than relaxed constraint?
iv) In Figure 4—figure supplement 1 the populations are presented separately, and differences among them are interpreted as differences in positive selection. But surely the divergence number must be massive, larger than the polymorphism number, so almost all of this variation is due to differences in constraint affecting the Pn/Ps ratio. Or is the method one that uses high frequency derived SNPs as evidence of positive selection? If so, how was ancestral state identified? This is probably not possible to do reliably if Ks>0.3. We need to see some raw numbers, and a lot more detail on the methods.
We apologise, the MK-based analyses methods were removed during the writing of the manuscript accidentally, as our results have been reconfigured as 2 manuscripts from one large manuscript, we have added them back. VCFtools was not used, we instead used the original SnIPRE code and SNPeff SNP assignments. We have also explained what we mean by difference from the background, specifically that we found the average statistical measure for nearby genes and compared each focal gene to that, to see how much the focal gene deferred from its surrounding average. We have attempted to answer all your points and questions in this new section of the Materials and methods. Average dS from the outgroups used was 0.133 and 0.279 between DiNV-Kallithea and DiNV-OrNV respectively, which we feel is enough divergence to adequately perform the analysis and is not too much to result in all differences being due to the Pn/Ps ratio. We also used these comparisons to identify the ancestral state, based on the Kallithea and OrNV variant. We have now included our measures of Dn/Ds and Pn/Ps as well as the estimated Selection Effect for our genes of interest as a supplementary figure. In this case we find that envelope proteins and AMPs have an excess of Dn/Ds per site compared to other genes, which is driving their elevated selection effect, while an excess of Pn/Ps and a deficit of Dn/Ds in unknown function DiNV genes is likely driving most other differences. We have chosen to use MK-based tests as opposed to Dn/Ds to identify if the changes are due to positive selection rather than relaxed constraint, as Dn/Ds is weighted by Pn/Ps to identify the proportion of substitutions fixed by selection.
Reviewer #3:
[…] 1) Throughout the manuscript, the authors switch back and forth between the present and past tense, which seems awkward from a stylistic point of, I'd suggest to stick to the past tense through-out.
We apologise for this, we have corrected this inconsistency.
2) The fact that the exact same 11 SNPs evolve multiple times independently in mostly the same order is interesting. The authors note that an expected alternative could be different SNPs emerging in the same genes, but they don't discuss the potential mechanism, by which the exact same SNPs seem to be preferred instead.
We have expanded our discussion of the compatibility, both from an evolution perspective and from a functional perspective. First, we hypothesize that the significantly associated SNPs have epistatic effects and thus must occur in a particular order to traverse the fitness landscape form Low type to the High type. In line with this, very few genes show signatures of adaptation in DiNV, which could further limit what beneficial mutations fix. We also discuss the different conditions which could lead to two types of DiNV being found in each population. We also hypothesize some methods these incompatibilities could occur, such as a change in the 19K/PIF-6 protein structure which limits membrane access for viral particles without this variant and could be rescued by a matching change in another protein. Or these two changes together are neutral but in combination prevent membrane access for the other viral type. We also find evidence of recombination events around the significantly associated SNPs, which could easily facilitate the recombination of variants onto the same background, without needing to traverse a lower fitness genotype combination.
3) The authors found little evidence of mixed infection with both the High and Low Types and conclude potential in-compatibility. Please expand on the potential mechanisms of this.
We have expanded our discussion of the compatibility, both from an evolution perspective and from a functional perspective, we have also highlighted similar studies where diverged viruses are unable to coninfect similar systems or are unable to produce viable recombinant particles.
4) Related to the above comment, the authors observe no "hybrids" of High and Low Types again, suggesting "incompatibility". Please discuss the difference between this genetic/genomic versus ecological incompatibility referred to above, as well as the potential link between these two types of incompatibilities.
We have included this in our expanded section in the Discussion. Briefly we discuss how the incompatibility may not be caused by a functional effect, as most SNPs are upstream of genes, but could be due to an effect of the expression of the proteins in combination on viral particle formation, host fitness or even viral ability to transmit between hosts. We do not have zero intermediate types, but we have explained that these could be steps between two fitness peaks as the few intermediates without low fitness.
5) Please rephrase the first sentence of the Discussion (what is the deference between "to better infect" and "optimizing the infection"?).
We have attempted to clarify that we mean that viruses are evolving according to the transmission/virulence trade-off and so are not just evolving to propagate within hosts but to transmit to others, which usually requires the virus not kill the host before transmission can occur.
https://doi.org/10.7554/eLife.58931.sa2Article and author information
Author details
Funding
KU CMADP (P20 GM103638)
- Tom Hill
- Robert L Unckless
K-INBRE (P20 GM103418)
- Tom Hill
National Institutes of Health (R00 GM114714)
- Robert L Unckless
National Institutes of Health (R01 AI139154)
- Robert L Unckless
National Science Foundation (DEB-1737824)
- Tom Hill
- Robert L Unckless
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
This work was completed with helpful discussion from Justin Blumenstiel, Joanne Chapman, John Kelly, Stuart MacDonald, Andrew Mongue and Carolyn Wessinger. We would especially like to thank Maria Orive, Kelly Dyer and Paul Ginsberg for helpful feedback in the writing of the manuscript and framing of the discussion. Collections were completed with assistance from Todd Schlenke, Paul Ginsberg, Kelly Dyer, Brandon Cooper, John Jaenike and the Southwest Research Station. We thank Brittny Smith and Jenny Hackett at the KU CMADP Genome Sequencing Core (NIH Grant P20 GM103638) and K-INBRE Bioinformatics Core for assistance in genome isolation, library preparation, sequencing and computational resources. This work was supported by a K-INBRE postdoctoral grant to TH (NIH Grant P20 GM103418). This work was also funded by NIH Grants R00 GM114714 and R01 AI139154 to RLU. D. falleni collection was funded by NSF grant DEB-1737824. Supplementary Data 1-9 is available in the following data dryad folder: https://doi.org/10.5061/dryad.2fqz612mh. (doi:10.5061/dryad.2fqz612mh). Additional data regarding D. innubila population genomics is available in the following FigShare folder: https://figshare.com/projects/innubila_population_genomics/87662.
Senior Editor
- Diethard Tautz, Max-Planck Institute for Evolutionary Biology, Germany
Reviewing Editor
- Dieter Ebert, University of Basel, Switzerland
Reviewer
- Sebastian Gagneux, University of Basel, Switzerland
Version history
- Received: May 14, 2020
- Accepted: October 28, 2020
- Accepted Manuscript published: October 28, 2020 (version 1)
- Version of Record published: November 24, 2020 (version 2)
Copyright
© 2020, Hill and Unckless
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics
-
- 1,430
- Page views
-
- 172
- Downloads
-
- 7
- Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading
-
- Evolutionary Biology
5-Methylcytosine (5mC) and DNA methyltransferases (DNMTs) are broadly conserved in eukaryotes but are also frequently lost during evolution. The mammalian SNF2 family ATPase HELLS and its plant ortholog DDM1 are critical for maintaining 5mC. Mutations in HELLS, its activator CDCA7, and the de novo DNA methyltransferase DNMT3B, cause immunodeficiency-centromeric instability-facial anomalies (ICF) syndrome, a genetic disorder associated with the loss of DNA methylation. We here examine the coevolution of CDCA7, HELLS and DNMTs. While DNMT3, the maintenance DNA methyltransferase DNMT1, HELLS, and CDCA7 are all highly conserved in vertebrates and green plants, they are frequently co-lost in other evolutionary clades. The presence-absence patterns of these genes are not random; almost all CDCA7 harboring eukaryote species also have HELLS and DNMT1 (or another maintenance methyltransferase, DNMT5). Coevolution of presence-absence patterns (CoPAP) analysis in Ecdysozoa further indicates coevolutionary linkages among CDCA7, HELLS, DNMT1 and its activator UHRF1. We hypothesize that CDCA7 becomes dispensable in species that lost HELLS or DNA methylation, and/or the loss of CDCA7 triggers the replacement of DNA methylation by other chromatin regulation mechanisms. Our study suggests that a unique specialized role of CDCA7 in HELLS-dependent DNA methylation maintenance is broadly inherited from the last eukaryotic common ancestor.
-
- Developmental Biology
- Evolutionary Biology
Gene expression has been employed for homologizing body regions across bilateria. The molecular comparison of vertebrate and fly brains has led to a number of disputed homology hypotheses. Data from the fly Drosophila melanogaster have recently been complemented by extensive data from the red flour beetle Tribolium castaneum with its more insect-typical development. In this review, we revisit the molecular mapping of the neuroectoderm of insects and vertebrates to reconsider homology hypotheses. We claim that the protocerebrum is non-segmental and homologous to the vertebrate fore- and midbrain. The boundary between antennal and ocular regions correspond to the vertebrate mid-hindbrain boundary while the deutocerebrum represents the anterior-most ganglion with serial homology to the trunk. The insect head placode is shares common embryonic origin with the vertebrate adenohypophyseal placode. Intriguingly, vertebrate eyes develop from a different region compared to the insect compound eyes calling organ homology into question. Finally, we suggest a molecular re-definition of the classic concepts of archi- and prosocerebrum.