{'fraction_non_alphanumeric': 0.03154285714285714, 'fraction_numerical': 0.02343896103896104, 'mean_word_length': 4.910106840230873, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Innate immune pathways act synergistically to constrain RNA virus evolution in Drosophila melanogaster

Host-pathogen interactions impose recurrent selective pressures that lead to constant adaptation and counter-adaptation in both competing species. Here, we sought to study this evolutionary arms-race and assessed the impact of the innate immune system on viral population diversity and evolution, using Drosophila melanogaster as model host and its natural pathogen Drosophila C virus (DCV). We isogenized eight fly genotypes generating animals defective for RNAi, Imd and Toll innate immune pathways as well as pathogen sensing and gut renewal pathways. Wild-type or mutant flies were then orally infected with DCV, and the virus was serially passaged ten times via reinfection in naïve flies. Viral population diversity was studied after each viral passage by high-throughput sequencing, and infection phenotypes were assessed at the beginning and at the end of the evolution experiment. We found that the absence of any of the various immune pathways studied increased viral genetic diversity while attenuating virulence. Strikingly, these effects were observed in a range of host factors described as having mainly antiviral or antibacterial functions. Together, our results indicate that the innate immune system as a whole, and not specific antiviral defense pathways in isolation, generally constrains viral diversity and evolution.


Introduction
Interaction between hosts and pathogens trigger defense and counter-defense mechanisms that often result in reciprocal adaptation and coevolution of both organisms 1 . Empirical evidence of such arms-races involving both species can be drawn from genome-wide analysis of hosts and pathogens and in experimental evolution settings. For example, evolutionary analysis of mammalian genomes has revealed evidence of host-virus coevolution between different retroviruses and antiviral factors 2,3 , and in plants, host resistance genes and virulence genes encoded by pathogens have been found to co-evolve 4 . Likewise, between bacteria and their infecting bacteriophages, experimental co-evolution studies resulted in the occurrence of genetic variants in both a bacterial lipopolysaccharide synthesis gene and the phage tail fiber gene which binds to lipopolysaccharide during adsorption 5 . In nematodes and their pathogenic bacteria, the number of toxin-expressing plasmids varies during adaptation to the host 6 .
In insects, analyses of sequences within and between Drosophila species have shown evidence of adaptive evolution in immunity related genes [7][8][9][10] . In a study that deep sequenced siRNAs from mosquitoes infected with West Nile virus, it was found that the regions of the viral genome more intensively targeted by RNA interference (RNAi) contained a higher number of mutations than genomic regions less affected by this pathway, suggesting that this antiviral defense mechanism imposes a selective pressure on the viral population 11 . Similar observations on the selective pressure imposed by the RNAi pathway on viral evolution have been made in plant-and human-infecting viruses [12][13][14][15][16] . Drosophila melanogaster is a well-studied insect model to decipher virus-host interactions and therefore the impact of host antiviral immunity on viral diversity and evolution. Different Drosophila immune pathways and mechanisms are involved in antiviral defense 17,18 . As is the case for all invertebrates, defense against pathogens in Drosophila relies on innate immunity, which constitutes the first and only defense against microbes. Innate immunity is characterized by the recognition of pathogen derived molecules, called pathogen-associated molecular patterns (PAMPs), by host encoded receptors (pathogen recognition receptors -PRRs), which leads to a rapid defense response.
The RNAi mechanism is known to play a central role in Drosophila antiviral defense, mainly through the action of the small interfering (si) RNA pathway [19][20][21][22] . Antiviral RNAi is triggered by virtually all insect-infecting viruses, resulting in targeting of the viral genome in a sequence-specific manner to control infection. Several other pathways have antiviral properties in flies, but their roles against viruses seem to be virus specific. The Toll and Imd (Immune deficiency) pathways, originally described to be involved in antibacterial and antifungal responses, have been shown to play a role in antiviral defense against Drosophila C virus (DCV), cricket paralysis virus (CrPV), Drosophila X virus, Nora virus, and Flock house virus [23][24][25][26] . The Janus kinase signal transducers and activators of transcription (JAK-STAT) pathway can be activated upon DCV or CrPV infection in flies, triggering the expression of antiviral factors 27,28 .
DCV, a positive-sense single stranded RNA virus from the genus Cripavirus within the Dicistrioviridae family and Picornavirales order 29 , is a well characterized natural pathogen of the fruit fly that can be found in laboratory and wild populations 30 . As for many other Drosophila-infecting viruses, defense against DCV depends on the joint action of different innate immune pathways and mechanisms. RNAi, Toll and Imd pathways, but also the protein encoded by the gene Vago, play a role in the defense against this virus 20, [24][25][26][27][31][32][33] .
DCV is thought to be naturally acquired by ingestion 30,34,35 . For orally acquired pathogens, the digestive tract, and the gut in particular, represents the first host defense barrier. Despite many studies using oral bacterial infections 36 , the role of gut-specific antiviral responses in Drosophila is not fully understood. Responses triggered against bacterial pathogens in the gut include the production of reactive oxygen species and antimicrobial peptides, as well as tissue repair and regeneration mechanisms 37 . Furthermore, the maintenance of gut homeostasis after tissue damage caused by pathogenic bacteria relies on the activity of JAK-STAT and epidermal growth factor receptor (EGFR) pathways, amongst others [37][38][39] . In the hallmark of viral infections, a role of the Imd and extracellular-signal-regulated kinase (ERK) pathways in the antiviral response in the gut has been suggested 24,40 . It is important to note that, like many other RNA viruses with error-prone polymerases and fast replication kinetics, DCV exists as large populations composed of a cloud of genetically related mutant variants known as viral quasispecies or mutant swarms 41 . Viral mutant swarms constitute a dynamic repertoire of genetic and phenotypic variability that renders great adaptability.
In this work, we leveraged the vast knowledge on antiviral mechanisms, the extensive genetic tool-box available for D. melanogaster, the intrinsic variability of the DCV mutant swarm, and the great depth power of next generation sequencing to study the impact of innate immune pathways on viral diversity and evolution. We aimed to determine not only if each pathway has a specific impact on the selective pressure imposed to DCV mutant swarms, but also their relative impact. In addition, we investigated possible links between selected viral variants (viral function) and specific defense mechanisms. Our results with infections in flies defective for several immune pathways show that the host genotype has an impact on viral genetic diversity regardless of the immune pathway being affected and this is accompanied by an attenuation of the virulence along evolutionary passages. We also describe complex mutation dynamics, with several examples of clonal interference in which increases in frequency of adaptive mutations have been displaced by other mutations of stronger effect that arose in different genetic backgrounds. Overall, our results highlight that innate immune pathways constrain RNA virus evolution and further demonstrate that antiviral responses in Drosophila are likely polygenic.

Production of fly mutant lines for innate immune pathways
To determine the impact of the innate immune system on virus population diversity and evolution, we selected fly lines with impaired function in genes belonging to most of the Drosophila innate immune pathways: RNAi, Toll and Imd. We selected genes involved both upstream and downstream of the immune pathways, such as receptors or ligands that trigger the immune response, and effectors of the response (Figure 1a). For the RNAi pathway, Dicer 2 (Dcr-2) and Argonaute 2 (Ago-2); for the Toll pathway, the ligand of Toll receptor Spätzle (spz), and the NF-κB transcription factor Dorsal-related immunity factor (Dif); for the Imd pathway, the NF-κB transcription factor Relish (Rel). We also added to the study the host factor Vago (Vago), that is upregulated during viral infections in a Dicer 2-dependent manner. Because DCV is orally acquired, and to explore the impact of gut homeostasis on the antiviral response, a mutant line for Epidermal growth factor receptor (Egfr), a gene involved in gut epithelium renewal, was also included in our panel. With the exception of Egfr and Dif, all of the selected genes were previously described to play an antiviral role against DCV infection [19][20][21][23][24][25] . It is important to mention that, in contrast to the RNAi antiviral mechanism that relies on the direct interaction between the components of the RNAi pathway and the viral genome, the molecular mechanisms underlying the antiviral responses mediated by Toll, Imd, and Vago in Drosophila remain largely unknow.
To reduce genetic variation due to differences in genetic background, mutant flies were isogenized prior to beginning viral evolution experiments. Homozygous loss-of-function lines for Dcr-2 (Dcr-2 L811fsX and Dcr-2 R416X ), Ago-2 (Ago-2 414 ), spz (spz 2 ), Dif (Dif 1 ), Rel (Rel E20 ), and Vago (Vago ΔM10 ) and a hypomorphic mutant line for Egfr (Egfr t1 ) were produced in the same genetic background by crossing parental lines at least 10 times to w 1118 flies. Infection phenotypes of the newly produced fly lines were characterized by following their survival after inoculation with DCV by intrathoracic injection (Supplementary Figure 1a). As previously described, Dcr-2 L811fsX/L811fsX , Dcr-2 R416X/R416X and Ago-2 414/414 mutants infected with DCV died faster than w 1118 flies 20,21 , as well as Vago ΔM10/ΔM10 mutants 33 . Toll pathway mutants spz 2/2 and Dif 1/1 and Imd pathway mutant Rel E20/E20 were less sensitive to DCV infection than w 1118 flies as they died later than w 1118 flies (Supplementary Figure 1a); however, these mutants maintained the previously observed increased susceptibility to infection by Gram + and Gram -bacteria, respectively (Supplementary Figure 1b and 1c). No difference in virus-induced mortality was found between w 1118 and Egfr t1/t1 mutant flies (Supplementary Figure 1a). This set of isogenic mutant flies with contrasting phenotypes to DCV infection provided us with the host model system to perform the viral evolution experiment.

Experimental DCV evolution
To study the impact of innate immune pathways on virus population diversity and evolution, DCV from a viral stock was serially passaged (P = 1 to P = 10) in w 1118 flies and in the isogenic innate immune deficient fly lines (Figure 1a and 1b). DCV population diversity was studied after each passage by next generation sequencing (NGS) and DCV virulence was analyzed at the beginning and at the end of the evolution experiment.
To follow viral infection during the course of the experiment, viral load (TCID 50 ) was determined by end point dilution and prevalence (percentage of flies positive for TCID 50 ) was calculated for all passages in individual flies from DCV contaminated cages. We found that for most fly genotypes and for both biological replicates, DCV infection prevailed along the 10 viral passages (Extended Data Figure 1a and  Europe PMC Funders Author Manuscripts both biological replicates, while viral load in the other fly genotypes remained relatively stable (Durbin-Watson test in the range 1.5 -2.5) for at least one of the biological replicates (Extended Data Figure 1b). The negative strand of the DCV genome was detected in P = 10 in all genotypes and biological replicates, confirming that active viral replication occurred for the duration of the evolution experiment (Extended Data Figure 1c). Whether remnants of non-replicating virus remained in the fly surface was not assessed. Of note, the DCV stock was experimentally introduced to the system only once, to start the P = 1.
To assess the impact that fly genotype, biological replicate, and viral passage has on viral loads, the log-transformed TCID 50 values from each fly genotype (Extended Data Figure  1d) were fitted to the generalized linear model ( the differences in viral load among experimental blocks observed for a particular host genotype also depended on the evolutionary passages, although once again the effect could be considered as minor η P 2 = 1.49% . Next, we evaluated whether differences exist in viral load between immune competent (w 1118 ) and the different mutant fly genotypes. In all eight cases, DCV accumulated to significantly higher levels in the immune deficient flies than in the wild-type flies (p < 0.001), with the smallest significant difference corresponding to viral populations replicating in Rel E20/E20 and Dif 1/1 and the largest to those replicating in Egfr t1/t1 and Dcr-2 R416X/R416X (Extended Data Figure 1d).
Overall, these results show that in both immune competent (w 1118 ) and immune deficient flies, DCV oral infection was maintained along passages and confirm that mutant flies are more permissive to DCV infection.

Viral nucleotide diversity increases in the absence of a fully functional immune response
To look into the selective pressure imposed by the Drosophila innate immune pathways on DCV population variation and dynamics, we analyzed virus genome diversity after each passage. Half of the population of infected flies was used to sequence the full-length DCV genome by NGS (Figure 1b and 1c). The viral stocks used to start the experiment, S2 DCV stock and DCV stock, were also sequenced (see Methods section). Analysis for the NGS data was performed using the computational pipeline Viral Variance Analysis (ViVan) 42 . Sequence coverage was at least 8,000 reads per position on the genome. To determine the error rate of the sequencing procedure, including library preparation, four sequencing technical replicates of the S2 DCV stock were used (Supplementary Figure  2). An allele frequency threshold of 0.0028 was used for all subsequent analyses based on variant detection and frequency correlation between technical replicates (see Methods section). We next calculated the site-averaged nucleotide diversity (π) on all polymorphic sites (n = 1869) across the full-length viral genome and present in the full dataset ( Figure 2), with the aim of determining if the lack of activity of a given innate immune pathway had an impact on viral population genetic diversity, in terms of size of the viral mutant swarm.
First, we asked if there was any difference in DCV population diversity and dynamics between the different fly genotypes along the complete evolution experiment. To answer this question, we analyzed if the host genotype, viral passages, biological replicate, or the interactions between these factors had an impact on the evolution of viral population diversity, considering the full-length DCV genome, across all passages. We found that only the fly genotype had a statistically significant impact on π (χ 2 = 25.545, 8 d.f., p = 0.001) (  Figure 2 and Supplementary Table 1).
Next, we wondered if the general differences observed in viral nucleotide diversity, between fly genotypes were associated with a particular viral genomic region (i.e., if a determined viral function was affected during the evolution experiment) ( Figure 1c). Of note, the intergenic region internal ribosome entry site (IGR IRES) was not included in the analysis because its lack of genetic variation prevented us from determining its nucleotide diversity value. We found that the fly genotype had a statistically significant effect on the nucleotide diversity found in each DCV genomic region (χ 2 = 27.178, 8 d.f., p < 0.001), which further differed between each specific viral genomic region (χ 2 = 11.698, 8 d.f., p = 0.008). As a second-order interaction, an effect of the fly genotype and the biological replicate was found (χ 2 = 16.314, 8 d.f., p = 0.038) ( Table 1). Comparison of viral genetic diversity within the genomic regions allowed us to distinguish three main groups: group 1 (less diversity): 3'UTR; group 2: 5'UTR IRES; and group 3 (more diversity): ORF1 and ORF2 (Extended Data Figure 2 and Supplementary Table 1).
Finally, we wondered if viral diversity evolved from the starting viral stock (DCV stock) in each fly genotype. π present in P = 1, P = 5 and P = 10 was compared between fly genotypes and with the diversity present in the DCV stock. We found that pairwise comparisons of viral nucleotide diversity present in each fly genotype in P = 1, between each other and versus DCV stock, yield no statistically significant difference (p = 1.000) (Supplementary Table 1). In P = 5 viral diversity was reduced only in w 1118 (Group 1/2; p = 0.026 and p = 0.032) compared to the starting viral stock (Extended Data Figure 2 and Supplementary Table 1). In P = 10, viral nucleotide diversity present in w 1118 (Group 1; p = 0.032 and p = 0.041), spz 2/2 (Group 1; p = 0.020 and p = 0.025), Dif 1/1 (Group 1; p = 0.005 and p = 0.006) and Rel E20/E20 (Group 1/2; p = 0.046) mutant flies was reduced when compared to DCV diversity from the DCV stock (Extended Data Figure 2 and Supplementary Table 1).
Altogether, the results show that the absence of a fully functional immune system results in an increase of viral population diversity that remains constant along passages. They also show that the coding regions of the virus are more prone to accumulate variation than the non-coding regions where regulatory elements are present.

Viral population diversity derives from preexisting standing genetic variation
Next, we examined if the levels of viral diversity observed in DCV populations from innate immune mutants compared to the w 1118 line were accompanied with the fixation of particular genetic changes in the mutant swarms, and whether (i) whether these changes can be associated with fitness effects and (ii) whether potentially adaptive mutations arose in response to particular immune responses. To do so, we estimated the selection coefficients for each SNP using their variation in frequency across evolutionary time ( Figure 3 and Extended Data Figure 3), using a classic population genetics approach 43 (Table 2). Thirtysix SNPs yielded significant estimates of selection coefficients (this number reduces to 10 if a stricter FDR correction is applied; Table 2). Twenty-one of them were already detected in the ancestral S2 DCV stock, henceforth a maximum of 15 new SNPs might have arisen during the evolution experiment. Estimated selection coefficients for all these SNPs ranged between −0.304 per passage (synonymous mutation RdRp/C5713U) and 1.204 per passage (VP2/G6311C nonsynonymous change R16P), with a median value of 0.286 per passage (interquartile rank = 0.265). Nine mutations were observed in more than one lineage (range 2 -7 times), with synonymous mutations VP3/U7824C appearing in seven lineages of six different host genotypes and mutation 5'UTR/A280U in five lineages of five host genotypes ( Table 2). These nine SNPs were all present in the S2 DCV stock. Indeed, the frequency of SNPs among evolving lineages is significantly correlated with their frequency in the ancestral S2 DCV stock (Pearson's r = 0.401, 36 df, p = 0.013), but not with their measured fitness effect (r = −0.091, 36 df, p = 0.588).
An interesting question is whether the fitness effects associated with each of these nine SNPs were the same across all genotypes or, conversely, whether fitness effects were host genotype-dependent. To test this hypothesis, we performed one-way ANOVA tests comparing fitness effects (Table 2) across the corresponding host genotypes. In all cases, significant differences were observed (F ≥ 15.637 and p ≤ 0.001, and ≥ 93.99% of total observed variance in fitness effects explained by true genetic differences among host genotypes), supporting the notion that fitness effects are indeed host-genotype dependent.
A pertinent example is the case of the synonymous mutation VP3/U7824C, which was the most prevalent mutation (F 6,45 = 158.862, p < 0.001, 99.37% of genetic variance).
In this case, a post hoc Bonferroni test shows that host genotypes can be classified into three groups according to the fitness effect of this SNP. In genotypes Dcr-2 R416X/R416X and Rel E20/E20 , the mutation has a deleterious effect (on average, −0.2260 per passage); in genotypes Egfr t1/t1 and Vago ΔM10/ΔM10 , the mutation is moderately beneficial (on average, 0.1257 per passage; and in genotypes w 1118 and Ago-2 414/414 , the mutation had a strong beneficial effect (on average, 0.502 per passage).
As shown in Figure 3 and Extended Data Figure 3a, some SNPs show a strong parallelism in their temporal dynamics, suggesting they might be linked into haplotypes. This is particularly relevant for mutations shown in Table 2. To test this possibility, we computed all pairwise Pearson correlation coefficients between mutation frequencies along evolutionary time. The results of these analyses are shown in Extended Data Figure 3b to 3k as heatmaps.
Again, as an illustrative example, we discuss here the case of the viral population BR2 evolved in Ago-2 414/414 (Extended Data Figure 3d).  Figure 3b). These two cases, as well as populations BR1 evolved in Rel E20/E20 , BR2 evolved in spz 2/2 and BR1 and BR2 evolved in vago ΔM10/ΔM10 illustrate some examples of haplotypes (Extended Data Figure 3f, 3e, 3h, and 3i). Other viral populations, especially those evolved in Egfr t1/t1 flies, show much more complex patterns (Extended Data Figure 3j and 3k) in which haplotypes change over time by acquiring de novo mutations.
When mapping the 36 SNPs found to have significant estimates of selection coefficients in the viral genome (Table 2 and Extended Data Figure 4), we found that two mapped to the 5'UTR IRES, twelve to ORF1, one to the IGR IRES, 20 to ORF2, and one to the 3'UTR.
Of the twelve mutations observed in ORF1, which encodes the non-structural proteins, four mapped to the 3C viral protease and five to the RdRp. Only one of these mutations in the 3C viral protease was non-synonymous. Of the 20 mutations in ORF2, which encodes the viral structural proteins, eight mapped to VP2, five to VP3, and seven to VP1. These correspond to the three majors predicted DCV capsid proteins.
Taken together, these results show that viral population diversity over these ten in vivo passages mainly derived from preexisting standing genetic variation in the ancestral DCV population. Furthermore, temporal dynamics of population diversity were linked to the fly genotype in which the virus evolved.

DCV virulence decreases along passages in the absence of immune pathways
Finally, we wondered if DCV virulence varied among each lineage in the different fly genotypes. Infectious DCV stocks were produced from viral passages P = 1 and P = 10 and from all fly genotypes. Because the viral evolution experiment was performed by DCV orofecal transmission, we first evaluated DCV virulence by feeding w 1118 flies with DCV stocks derived from P = 1 or P = 10; survival was evaluated from each fly genotype. We found that only a small proportion of flies (between 5% and 20%) succumbed to DCV infection, and no statistically significant differences in mortality were found between mock and virus infected flies, regardless of viral passage or fly genotype (Supplementary Figure  3). This is in agreement with previously published works showing that DCV oral infections are cleared in w 1118 flies 31 . We next decided to investigate the evolution of virulence by intrathoracic inoculation of DCV stocks. We found that w 1118 flies were less sensitive to viral infection when inoculated with DCV stocks derived from P = 10 since their median survival time was longer than those inoculated with stocks from P = 1 for most DCV stock origins (Figure 4a and Supplementary Table 2). Notable exceptions were DCV stocks from BR2 of Vago ΔM10/ΔM10 mutant flies, for which w 1118 flies were more sensitive to P = 10 than to P = 1, and stocks from BR1 of spz 2/2 and BR2 of Egfr t1/t1 mutant flies, for which no difference in median survival time after infection with DCV between P = 1 and P = 10 was detected.
A fundamental question in evolutionary biology is the role that past evolutionary events may have in the outcome of evolution 44 . If ongoing evolution is strongly contingent with past evolutionary events, ancestral phenotypic differences should be retained to some extent, while if other evolutionary forces such as selection and stochastic events (i.e., mutation and genetic drift) dominate, then ancestral differences can be eroded and, in the extreme case, even fully removed. Here, we observed significant differences in the performance of the ancestral DCV across the eight host genotypes. To test whether these differences are still observable in the evolved population, we compared the median survival time (Figure 5a and Supplementary Table 2) for DCV populations isolated at the beginning of the evolution experiment P = 1 and at the end P = 10 ( Figure 4b). Under the null hypothesis of strong historical contingency, it is expected that data will fit to a regression line of slope 1 and intercepting the ordinate axis at 0. However, if ancestral differences have been removed, data would fit significantly better to a regression line with a slope smaller than one and with an intercept greater than zero 44 . Figure 4b shows the data and its fit to the null hypothesis (solid black line) and the alternative hypothesis (dashed red line). A partial F-test shows that adding an intercept to the regression equation significantly improves the fit (F 1,16 = 28.437, p < 0.001), thus supporting the notion that ancestral differences among host genotypes have been removed by the action of subsequent adaptation, that is, the fixation of beneficial mutations.

Discussion
In this work we aimed at determining the overall impact of innate immunity on viral evolution. Based on the arms-race hypothesis, we speculated that if a given host defense mechanism imposes a specific selective pressure on a particular pathogen function, the absence of this defense mechanism would result in the relaxation of the selective constraint, which would in turn be detectable in the pathogen at the genomic and phenotypic levels.
We found that viral population diversity evolved differently according to each fly genotype; however, viral population diversity mostly derives from ancestral standing genetic variation (i.e., few "new" mutations were selected). Our results further confirm the polygenic nature of antiviral responses; there is not a specific, main immune defense mechanism against a particular virus, but instead a repertoire of defense mechanisms that are triggered after infection and that might interact with each other.
Our results are compatible with a pervasive presence of clonal interference. In the absence of sexual reproduction, clonal interference is the process by which beneficial alleles originated in different clades within a population compete with each other, resulting in one of them reaching fixation. Subsequently, the outcompeted beneficial allele may appear in the new dominant genetic background and, assuming no negative epistasis among both loci, become fixed. As a consequence, beneficial mutations may fix sequentially, thus slowing down the rate of adaptation 45 . Given their large effective population size and high mutation rates, viral populations are expected to contain considerable amounts of potentially beneficial standing variation, making them prone to clonal interference. Indeed, it has been previously shown to operate in experimental populations of vesicular stomatitis virus adapting to cell cultures 46,47 , in bacteriophage ϕX174 populations adapting to harsh saline environments 48 , in tobacco etch virus adapting to novel plant host species 49 , among HIV-1 escape variants within individual patients 50 , and also at the epidemiological level among influenza A virus lineages diversifying antigenically 51 . In our own results, clonal interference can be observed in populations BR1 evolved in Dcr-2 L811fsX/L811fsX , BR1 evolved in Ago-2 414/414 , BR1 evolved in spz 2/2 , BR2 evolved in Rel E20/E20 , and BR2 evolved in Vago ΔM10/ΔM10 . All of these viral populations share similar patterns in which some beneficial allele (or haplotypes) rose in frequency, reached a peak at some intermediate passage, then declined in frequency and were finally outcompeted by a different beneficial mutation (or haplotype) that had lower initial frequency. For example, the nonsynonymous mutation VP2/G6931A (A223T) appeared de novo in population BR1 evolved in spz 2/2 , and outcompeted several mutations likely linked in a haplotype ( Figure 3). Tightly linked to clonal interference is the concept of leap-frogging 52 , in which the beneficial mutation that ends up dominating the population is less genetically related to the previously dominant haplotype than to the common ancestor of both ( Figure 3). The VP2/G6931A mutation illustrates this case well, as it appeared in a genetic background that was minoritarian rather than in the dominant one. Likewise, the mutation VP2/G6311C (R16P), observed in BR1 evolved in w 1118 flies, appeared in a low frequency genetic background different from the most abundant one in previous passages. Finally, the haplotype containing five different mutations observed in BR2 evolved in spz 2/2 became dominant in frequency after P = 6, outcompeting two other mutations that were dominating the population until then.
The existence and fixation of haplotypes along our evolution experiment deserves further discussion. Linked mutations generate three possible interference effects 53 . First, all mutations might contribute additively, or may be involved in positive epistasis, to the fitness of the haplotype as a whole, thus increasing its chances to become fixed. Second, hitchhiking and genetic draft may occur, by which deleterious or neutral alleles are driven to fixation along with a linked beneficial allele. Third, there may be background selection by which the spread of a beneficial allele is impeded, or at least delayed, owing to the presence of linked deleterious alleles. For instance, we can hypothesize that haplotype VP3/ U7824C-VP1/C8227U-VP1/C8424U, which swept to fixation in population BR2 evolved in Ago-2 414/414 , may represent a case of genetic draft: two synonymous mutations, potentially neutral, linked to a nonsynonymous one that may be the actual target of selection. Yet, the lack of an infectious clone for DCV does not allow us to test this hypothesis.
Some of the mutations we found to be associated with positive selection coefficients were synonymous changes (Table 3). However, equating synonymous mutations with neutral mutations in compacted RNA genomes has proven to be misleading 54,55 . Selection operates at different levels of a virus's infection cycle, and not all these levels necessarily depend on the amino acid sequence of encoded proteins. For instance, a lack of matching between virus and host codon usages would slowdown translational speed and efficiency 56 ; mutations affecting the folding of regulatory secondary structures at noncoding regions would affect the interaction with host and viral factors and thus impact the expression of downstream genes (e.g., mutations 5'UTR/A280U, IGR/A6108G and 3'UTR/U9163A all with significant fitness effects - Table 3) 57 ; or evasion from antiviral RNAi defenses by changing the most important relevant sites in the target of siRNAs 12,13 .
It is interesting to observe that viral diversity in mutants for antiviral RNAi, whose mode of action relies on a direct interaction with the viral genome, did not display increased diversity when compared to mutants from the other immune pathways. One could expect that the release of the selective pressure that RNAi exerts on the virus genome may allow for the appearance of mutations in the viral suppressor of RNAi. Nonetheless, we did not observe such a change, possibly because the RNAi suppressor in DCV shares the first 99 amino acids of the RdRp 58,59 and mutations would affect polymerase activity. The antiviral action of the other immune pathways remains still unknown and may even be indirect; for example the known roles of Imd, Toll, and Egfr pathways in controlling fly microbiota 38,39 might possibly affect the prevalence of virus infections. In this regard, it is important to highlight that the diversity of DCV in the Dif 1/1 mutant (Toll pathway, already described not to have an impact on DCV defense 25 ), was indistinguishable from w 1118 , pointing to the specificalthough uncharacterized -antiviral functions of these other immune pathways.
Another consideration when interpreting our results is the nature of the virus stock used. This virus stock has been maintained for years in Drosophila S2 cells. The observation that viral population diversity decreased along passages in w 1118 flies, highlights the strength of the selective forces that constrain the virus from adapting to a new environment. During the successive passages, in the absence of a given immune response, the capacity of the virus to evolve will be determined by a combination of two factors: the adaptation to the new environment (constraining factor) and the lack of immune response (relaxation factor). Because DCV replication is significantly increased in immune deficient mutants, the potential for population diversification is higher. This effect is clearly observed in w 1118 flies where the virus is "only" adapting to the new environment and DCV populations evolved in w 1118 flies show less variation than all other lineages. Future experimental evolution studies using viral stocks derived from flies, instead of cell cultures, are warranted to address this topic.
In a study published recently 60 , Navarro et al. used Arabidopsis thaliana and turnip mosaic virus to carry out experimental virus evolution assays with a similar design to ours. In their work, the authors used plant mutants compromised in their antiviral response (more permissive to viral infection) or with an enhanced antiviral response (less permissive to viral infection) and allowed the virus to evolve for 12 passages. Similarly, to what we found in the D. melanogaster -DCV system, the authors showed that viral population evolutions dynamics, as well as viral loads, depend on host genotype. Interestingly, a reduction of ancestral genetic variation regardless of the immune pathway affected was also clearly observed, in agreement with our observations.
Taken together, our results point to the concerted action of the different immune pathways to limit viral evolution. Response to infection does not simply consist of activating immune pathways, it also encompasses a broad range of physiological consequences including metabolic adaptations, stress responses and tissue repair. Critically, upon infection, the homeostatic regulation of these pathways is altered. However, such alterations do not always result in increased disease severity and in fact can even lead to improved survival (or health) despite active virus replication.

Fly strains and husbandry
Flies were maintained on a standard cornmeal diet (Bloomington) at a constant temperature of 25 °C. All fly lines were cleaned of possible chronic infections (viruses and Wolbachia) as described previously 61

Virus production and titration
The S2 DCV stock used to start the experiment was prepared in S2 cells. Cells were maintained in Schneider culture medium and at 25 °C and observed daily. Cells were harvested when cytopathic effects were detected, then frozen at −80 °C, thawed on ice and centrifuged for 15 min at 15,000 g at 4 °C. The supernatant was recovered, aliquoted and stored at −80 °C. Viral stocks were tittered in S2 cells, determined using the end-point dilution method and expressed as fifty-percent tissue culture infective dose (TCID 50 ) 68 .
To produce the DCV stocks from passages P = 1 and P = 10 from the evolution experiment half of the population of flies infected with DCV from each fly genotype (approx. 250 flies) was homogenized in 1 × PBS, homogenates were frozen at −80 °C, then thawed on ice, centrifuged to discard the tissue debris, supernatant was recovered and filtered to discard bacteria contamination, then aliquoted and stored at −80 °C. Viral stocks were tittered in S2 cells using the end-point dilution method and expressed as fifty-percent tissue culture infective dose (TCID 50 ).

Viral and bacterial infections and survival analysis
To characterize the isogenized fly lines, 4 to 5 day old female flies were intrathoracically injected with a Nanoject II apparatus (Drummond Scientific) with 50 nl of the pathogen suspension. For DCV infections, a suspension of 10 TCID 50 units of DCV in 10 mM Tris buffer, pH 7 was used. An injection of the same volume of 10 mM Tris, pH 7 served as a mock-infected control. Infected flies were kept at 25 °C, transferred into fresh vials every 2 days and number of dead flies was scored daily. For bacterial infections, 50 nl of suspensions in 1× PBS buffer, pH 7, of optical density (OD) = 10 for Enterococcus faecalis, and of OD = 200 for Erwinia carotovora carotovora 15 (Ecc15) were used. An injection of the same volume of 1× PBS buffer served as a mock-infected control. Flies infected with E. faecalis were kept at 25 °C, and flies infected with Ecc15 were kept at 29 °C. Flies were transferred into fresh vials every 2 days and number of dead flies was scored daily.

Virus experimental evolution
To produce the starting DCV stock (DCV stock) 5 to 6 days old w 1118 female flies were intrathoracically injected with 100 TCID 50 of DCV from a stock produced in S2 Drosophila cells (S2 DCV stock) or mock infected. At 4 dpi, N = 90 DCV infected flies (DCV stock) were placed in cages containing fresh medium, left for 3 days and then removed to place in this DCV-or mock-contaminated cages N = 500 5-to-6-day old wild type or mutant flies (males and females). Flies were allowed to feed ad libitum for 3 days (oral inoculation period), then moved to a clean cage for 1 day, and further placed into a new clean cage and left for 4 days, when they were harvested (P = 1). A new group of flies was then placed into the contaminated cages. This procedure was repeated 10 times (10 DCV Passages, P = 1 to P = 10) and replicated twice (biological replicates BR1 and BR2). The total amount of flies from each passage, fly genotype, and biological replicate was collected and randomly divided in halves (approx. 250 flies), one half was used to extract total RNA and produce the NGS libraries and the other half to produce viral stocks to evaluate DCV virulence.

Characterization of infection during passages
Individual flies from each passage were anesthetized and homogenized in 100 ml of 1' PBS buffer. The tubes containing the homogenates were centrifuged for 5 min at 15,000×g at 4 °C to discard the tissue debris. The supernatant was recovered and used to determine viral load (TCID 50 ) by end point dilution and prevalence (percentage of flies positive for TCID 50 ) for each fly genotype, viral passages and biological replicate.
For statistical analyses, TCID 50 data were transformed as T = log(TCID 50 + 1) and then fitted to a generalized linear model in which fly genotype (G) and BR (B) were treated as orthogonal factors. G was considered as a fixed effects factor whereas B was considered as a random effects factor. Evolutionary passage (P) was introduced in the model as a fixed effects covariable. We also considered second and third order interactions between the two factors and the covariable. The model equation thus reads: T ijk (P ) τ + P + G i + B j + (P × G) i + (P × B) j + (G × B) ij + (P × G × B) ij + ε ijk .
Where T ijk (P) is the transformed TCID 50 observed for a particular titration assay k of BR j of fly genotype i, τ represents the grand mean value and ε ijk stands for the error assumed to be Gaussian distributed at every P. The significance of each term in the model was evaluated using a likelihood ratio test that follows a χ 2 probability distribution. The magnitude of the effects was evaluated using the η P 2 statistic (proportion of total variability in the traits vector attributable to each factor in the model; conventionally, values of η P 2 ≥ 0.15 are considered as large effects). These analyses were done using SPSS version 27 (IBM, Armonk, NY).

Detection of negative strand DCV RNA by strand-specific RT-qPCR (ssRT-qPCR)
To determine the amount of negative strand DCV RNA present in the viral stocks produced from each fly genotype in P = 10, S2 DCV stock, and DCV stock, total RNA was extracted from the DCV stocks produced from P = 10 (all fly genotypes, both biological replicates) and from the DCV stocks used to start the experiment. ssRT-qPCR was performed with these RNA samples essentially as described 35 . We used 800 ng of RNA to perform reverse transcription with SuperScript II reverse transcriptase (Invitrogen) according to the manufacturer's instructions, with the exception that primer annealing occurred at 70 °C and cDNA synthesis occurred at 50 °C for 30 minutes. Reverse transcription was performed using a forward primer containing a non-target tag sequence (DCV_tag_F: AATTCAAGCTCGTCTTCCTCGAGGCTGTGTTTGCGCGAAG) A standard curve was produced by reverse transcription of a tenfold dilution series (from 10 8 to 10 3 copies per reaction) of in vitro transcribed RNA corresponding to a portion of the full-length negative strand DCV RNA. Following reverse transcription, cDNA was diluted 1:10 and used for for qPCR with the Luminaris Color HiGreen low ROX qPCR Master Mix (Thermo Scientific) according to the manufacturer's instructions. A forward primer containing the non-target tag sequence (Tag_qPCR_F: AATTCAAGCTCGTCTTCCTCG) and a a DCV-specific reverse primer (DCV_qPCR_R: AATGGCAAGCGCACACAATTA) were used for qPCR.

RNA extraction, cDNA synthesis and NGS library production
To produce the NGS libraries from the evolution experiment, half of the total population of flies infected with DCV from each fly genotype, viral passage and biological replicates (approx. 250 flies) was used. To produce the NGS libraries from the viral stock from S2 cells (S2 DCV stock), two different aliquots of the stocks were used. To produce the NGS libraries from the DCV stock (virus infecting w 1118 female flies used to contaminate the cages to start the evolution experiment), half of the population of the infected flies (approx. 800 flies: N = 90 flies/cage ' 9 fly genotypes ' 2 BR) was used. In all cases, total RNA was extracted using TRIzol reagent (Invitrogen) following the manufacturer's instructions and the final concentration was determined using a NanoDrop ND-1000 Spectrophotometer. 300 ng of total RNA were used to produce the cDNA using oligo ( Sequencing of DCV populations from Dif 1/1 mutant flies from P = 4 and P = 6 from BR1 and P = 8 from BR2 did not work.

Genetic diversity analyses
Variant frequency threshold-To determine the error rate of the sequencing procedure, including library preparation, four sequencing technical replicates from S2 DCV stock were used (Supplementary Figure 3a). First, pairwise comparison was done to identify the variant frequency threshold above which at least 95% of the variants were detected in both considered replicates (highest detection threshold = 0.0028). All variants above detection threshold were then correlated between each technical replicate to ensure good correlation between reported frequency values: the Pearson correlation coefficient between the detected frequency for variants was r ≥ 0.982 for all pairwise correlation (p < 0.001). The R packages used for these analysis were somewhere else described [69][70][71][72] .
Nucleotide diversity (π)-Nucleotide diversity of the viral population was computed using the following formula 73 : with D, the sequencing depth and p the frequency of the minority variant at each nucleotide site. For diallelic SNV, π ranges from 0 to 0.5 (both alleles at equal frequency). In the subsequent analyses, π was averaged over all polymorphic nucleotide sites of the DCV genome of each sample 74 . A site was considered polymorphic if at least one sample showed the presence of a nucleotide variant at said position of the DCV genome. Log 10transformed site-averaged π values were then compared between fly genotypes (orthogonal factor), biological replicates (orthogonal factor), passages (continuous variable) and genomic regions (orthogonal factor) and their interactions using a generalized linear model. The significance of each term in the model was evaluated using a likelihood ratio test that follows a χ 2 probability distribution.
Estimation of relative mutational fitness effects-We have followed the classic population genetics method described in Hartl and Clark (1989) 43 . In short, lets xl(t) be the frequency of a mutant allele (SNP) at genomic position l and passage t and, therefore, 1 − X l (t) the frequency of the wild-type allele. It holds that Haplotype inference was done using two different statistical approaches. First, by assessing the similarity between temporal dynamics of all possible pairs of loci. To this end, Pearson partial correlation coefficients (controlling for passages) were computed and their significance level corrected for multiple tests of the same null hypothesis using Benjamini and Hochberg (1995) 75 false discovery rate (FDR) method. Correlation coefficient matrices were visualized as heatmaps in which more similar alleles were clustered together. Second, we confirmed the results from the first method using the longitudinal variant allele frequency factorization problem (LVAFFP) method as implemented in CALDER 76 . LVAFFP generates spanning trees of a directed graph constructed from the variant allele frequencies.
The output of CALDER was used as input of TimeScape 77 to generate the Muller plots that illustrate the ancestry of mutations and haplotypes along the evolution experiment ( Figure  3).  (P = 1) DCV virulence. The dashed red line represents the linear regression, and the black line represents the expected relationship under the null hypothesis of ancestral differences in DCV virulence which are maintained after evolution despite noise introduced by random events (mutation and drift). Mongelli  Analysis of the impact of each experimental variable on the evolution of DCV nucleotide diversity (mean log 10 (π) per site).
The site-averaged nucleotide diversity (π) on all polymorphic sites (n = 1869) across the full-length viral genome was determined, the log10-transformed π values were fitted to the generalized linear model (GLM) and the impact of the variables determined by an Analysis of Deviance (Type III tests). BR: biological replicate; VP: viral passage; FG: fly genotype, GR: genomic region.  Table 2 Mutations for which significant estimates of fitness effects have been obtained.
For each mutation, we indicate whether it already existed in the DCV starting stocks (and at which frequency) or arose during the evolution experiment. We also provide the estimated selection coefficient, its SEM and statistical significance. Cases significant after FDR correction are marked with an asterisk.