The impact of artificial selection for Wolbachia-mediated dengue virus blocking on phage WO

Wolbachia is currently at the forefront of global efforts to control arbovirus transmission from the vector Aedes aegypti. The use of Wolbachia relies on two phenotypes—cytoplasmic incompatibility (CI), conferred by cifA and cifB genes in prophage WO, and Wolbachia-mediated pathogen blocking (WMPB). These traits allow for local, self-sustaining reductions in transmission of dengue (DENV) following release of Wolbachia-infected A. aegypti. Here, aided by previous artificial selection experiment that generated Low and High pathogen blocking lines, we examined the potential link between WMPB and phage WO. We found no evidence that Wolbachia or phage WO relative densities predict DENV blocking strength across selected lines. However, selection resulted in reduced phage WO relative density for the Low WMPB line. The Low blocking line was previously shown to have reduced fitness as a result of selection. Through subsequent genomic analyses, we demonstrate that SNP variation underpinning selection for low blocking led to elevated frequency of potential deleterious SNPs on chromosome 1. The key region on chromosome 1 contains genes relating to cell cycle regulation, oxidative stress, transcriptional pausing, among others, that may have cascading effects on Wolbachia intracellular environment. We hypothesize that reduction in phage WO may be driven by changes in the loci directly under selection for blocking, or by the accumulation of predicted deleterious alleles in linkage disequilibrium with blocking loci resulting from hitchhiking. For the Low line with fewer phage WO, we also detected reduced expression of cifA and cifB CI genes, with patterns of expression varying between somatic and reproductive tissues. In conclusion, we propose that artificial selection for WMPB trait had corresponding impacts on phage WO densities, and also the transcription of CI-causing genes. Future studies may include a more detailed analysis of the regions the A. aegypti chromosome 1’s ability to affect WMPB and other Wolbachia-associated intrinsic factors such as phage WO.


Introduction
Mosquitoes transmit many medically relevant pathogens [1]. Aedes aegypti is a relevant vector of arboviruses as it transmits dengue (DENV), Zika, chikungunya and yellow fever viruses. DENV infection causes approximately 100 million symptomatic infections per year [2], with an estimated global annual cost to health care systems of US $8-9 billion [3]. Currently, half the world's population is at risk of dengue infection. The risk will likely intensify, as the vector's geographical range is predicted to increase due to climate change and socioeconomic factors [2,4]. Without an effective antiviral drug or vaccine [5], vector control remains the best strategy to curb vector-borne disease transmission.
Classical vector control initiatives such as active removal of mosquito breeding sites and insecticide spraying have limitations such as poor community compliance and development of insecticide resistance [6][7][8]. As such, the development of alternative approaches like gene drives [9,10] and the use of bacteria and fungi for biocontrol [11] is critical. Amongst these, the use of Wolbachia is one of the most promising biological strategies for controlling arbovirus transmission in the field [12,13]. Wolbachia is a genus of maternally-inherited, endosymbiotic bacteria present in many arthropod species [14]. The bacteria can induce two main effects in mosquitoes that together make them useful for biological control. First, they reduce the replication of co-infecting arboviruses, a phenotype known as "Wolbachia-mediated pathogen blocking" (WMPB) [15,16]. Second, they manipulate host reproduction in a manner that enhances the symbiont's spread in arthropod populations through induction of a phenotype known as cytoplasmic incompatibility (CI) [17]. CI is typically characterized by reduced viability of embryos from matings specifically between Wolbachia-free females and Wolbachiainfected males [17,18]. Wolbachia-infected females, in contrast, reproduce successfully regardless of their mate's infection status and produce more infected offspring; hence, Wolbachia spreads. Numerous field releases of A. aegypti stably transinfected with Wolbachia, have demonstrated the symbiont's ability to spread, and more recently, to reduce the transmission of DENV to humans [12,13,19].
The genetic basis of wMel-induced CI was recently resolved to two genes, here referred as cytoplasmic incompatibility factor A (cifA) and cytoplasmic incompatibility factor B (cifB), in the eukaryotic association module of phage WO (named after Wolbachia) [20][21][22] (please refer to [23][24][25] for an in-depth discussion on the proposed nomenclature of CI genes according to phylogenetic and enzymatic data on the matter). Considering all insect hosts surveyed, phage WO is temperate and may exist as a resident in Wolbachia genomes containing prophage WO (integrated lysogenic form in Wolbachia chromosome) or as a lytic bacteriophage (infectious capsids capable of exiting Wolbachia) [26][27][28]. Prophage WO and/or its genetic remnants are present in at least 89% of surveyed Wolbachia strains [29]. Most Wolbachia genomes are known to contain multiple phage WO haplotypes, or genetic variants of the phage WO family. In most cases, defective forms of the prophage WO (unable to form a lytic particle) are generally accompanied by fully-intact phage [22,30]. For instance, the genome of wMel contains two prophage elements: a small pyocin-like element named WOMelA and a larger intact prophage named WOMelB. It is still unknown if WOMelB is capable of forming active lytic particles, and without complete reference genomes for both the integrated prophage and active phage, there is no way of precisely distinguishing between lysogenic and lytic particles. WOMelB modular structure is intact, however, and composed of two prophage fragments that complement each other, WOMelB1 (containing the head and tail modules), and WOMelB2 (harboring the baseplate module and recombinase) [30][31][32].
Although phage WO commonly exists in multiple copies [31,33] and comprises a significant percentage of Wolbachia's genome, only a single pair of cifA and cifB is required to recapitulate CI in transgenic D. melanogaster [21,34]. Furthermore, rescue of CI can occur by expression of cifA alone in ovaries [20,35]. While the expression of CI assists Wolbachia spread, WO is not entirely beneficial as described in the "phage density model". When in the lytic phase, mature phage WO particles lyse Wolbachia, leading to corresponding decreases in bacterial density and consequentially CI expression because the cif genes reside in the prophage region of Wolbachia.
The mechanism of WMPB in mosquitoes remains unsolved. Several factors are associated with WMPB [36][37][38][39]. Recently, we employed experimental evolution in an attempt to solve the mechanism of WMPB, as well as to predict its likely long-term evolutionary stability [40]. We artificially selected for both stronger (referred throughout the manuscript as High) and weaker (referred throughout the manuscript as Low) DENV blocking in mosquitoes infected with wMel Wolbachia. In the resulting selected lines, our group tested for differences in mosquito fitness between evolved populations, and the associated differences in SNPs in the mosquito genome [40,41]. The study revealed two novel major gene candidates, α-mannosidase, and cadherin, along with another 58 minor genes in A. aegypti whose SNP variation associated with changes in WMPB phenotype. It also revealed that mosquito lines selected for higher blocking also exhibited faster population growth rates, an association that may be due to pleiotropy between the two traits. This association may help natural selection maintain or drive strong blocking in populations.
While the mechanistic basis of both CI and pathogen blocking are two 'holy grails' of Wolbachia biology, it has not been explored whether there may be interactions between these two traits. Here, we took advantage of the above artificial selection experiment, where we generated lines exhibiting phenotypic extremes in blocking. Using DNA and RNA samples from these lines, we tested for correlations between blocking strength and SNP variation, and the loads of both Wolbachia and phage WO. Last, we examined cifA and cifB gene expression in tissues of the selected lines to explore the potential consequences of selection on blocking for the expression of genes known to impact the strength of CI [20].

Artificial selection on DENV blocking strength alters phage-Wolbachia associations
We quantified the parameters related to both Wolbachia and phage WO densities utilizing a set of samples previously collected as part of the original selection experiment by Ford et al. [40]. We provide this quantification for samples collected from two lines (each in triplicate) post four generations of selection for High and Low DENV blocking as per Ford et al. [40]. In this study, both the Ancestral base population and the Random selected lines were excluded from our analyses since neither of the comparisons act as true controls to the evolved lines (see Materials and methods sub-section "mosquitoes" for details).
We first examined the absolute abundances of both Wolbachia and phage WO in the evolved lines (Fig 1A and 1B) (High and Low blocking) and saw a strong positive correlation between the two variables, as expected. Next, we utilized a relative density measurement to explore the interaction between phage WO and Wolbachia (see materials and methods subsection "Wolbachia, phage and mosquito gene quantification" for details on quantification strategies). Relative phage WO densities can estimate lytic phage development [26] when the number of phage genomes relative to Wolbachia genomes exceeds the baseline expectation of one prophage WO copy for every wMel genome. Alternatively, molecular makers for phage WO genome rearrangements can be used to distinguish lytic versus lysogenic genomes [42], but that is not currently possible as the complete genome sequence for WOMelB phage particles has not yet been resolved. Whole-body analyses (minus ovaries) on individuals revealed that while Wolbachia relative density values were similar between the two evolved populations (Mann-Whitney U = 364, d.f. = 1; P = 0.3) (Fig 1C), phage WO relative densities were significantly lower in the Low blocking group (1.38 ± 0.21 -mean ± standard error) compared to High blocking (1.89 ± 0.62) (Mann-Whitney U = 63, d.f. = 1; P = 0.01) (Fig 1D). This finding indicates alterations on phage relative density between the evolved populations as a consequence of selection for WMPB.
We then tested for correlations between the relative densities of both Wolbachia and phage WO in selected populations using Spearman's correlation coefficient analysis, as to test for the phage density model [26]. We observed marginally non-significant, negative correlations for both the High (Spearman's ρ = -0.44, P = 0.06) and Low blocking groups (Spearman's ρ = -0.39, P = 0.12). The significance values of these correlations are very likely due to sample sizes, as combining both sets of selected lines resulted in a significant and expected negative correlation of relative phage WO density against Wolbachia density (Spearman's ρ = -0.37, P = 0.02) (Fig 2). Under the phage WO density model, high Wolbachia densities associate with baseline prophage WO densities of about 1 in this case, and low Wolbachia densities associate with higher phage WO densities likely due to phage WO replication and lysis of Wolbachia cells.

Low DENV blocking mosquitoes have reduced overall fitness and higher predicted deleterious mutation load
While we did not observe a linear relationship between phage/Wolbachia relative densities and DENV loads, there were lower phage densities in the Low lines alone. From the original study [40], the Low blocking line also exhibited a reduced population growth rate (r) compared to Whole body (minus ovaries) Wolbachia and phage relative densities for wMel-infected A. aegypti mosquitoes from the evolved lines after 4 generations of selection (High: dark blue-Low: yellow). Correlation between Wolbachia (absolute copy number of groEL) and phage WO (absolute copy number of orf7) abundances in whole body (minus ovaries) of wMel-infected A. aegypti mosquitoes selected for (A) High and (B) Low WMPB. Red, orange and blue dots represent replicated populations 1, 2 and 3, respectively. (C) Wolbachia densities are reported as density per host cell (absolute copy number of groEL relative to the absolute copy number of rps17) and (D) phage density as per Wolbachia cell (absolute copy number of orf7 relative to the absolute copy number of groEL). Each dot represents a single mosquito and Boxplot depict medians. Sample sizes for each group are depicted on each graph. Numbers "1", "2" and "3" represents distinct replicated populations of the same selected group. Significant comparisons expressed via P-values were determined via Spearman's rank correlation coefficient across evolved groups for graphs (A) and (B). Graphs (C) and (D) were analyzed via Mann-Whitney U test.
https://doi.org/10.1371/journal.pntd.0009637.g001 Correlation between Wolbachia and phage WO relative densities in whole body (minus ovaries) wMel-infected A. aegypti mosquitoes selected for High and Low WMPB. Each dot represents a single mosquito. Sample sizes are depicted within the figure for each group. Red, orange and blue dots represent replicated populations 1, 2 and 3, respectively. Spearman's rank correlation coefficient analyses were performed to describe the relationship between Wolbachia density (absolute copy number of groEL relative to the absolute copy number of rps17) and phage WO density (absolute copy number of orf7 relative to the absolute copy number of groEL).
https://doi.org/10.1371/journal.pntd.0009637.g002 the High blocking line. Significant SNP variation from the Poolseq approach (100 individuals per line) differentiated the High and Low Blocking lines, which may explain both differences in blocking as well as fitness, with potential consequences for phage WO. As such, we explored if there were elevated patterns of mildly predicted deleterious alleles in the selected lines and across the different chromosomes. While there is no evidence of A. aegypti genes that affect phage WO density, stress in Wolbachia due to external influences, such as temperature exposure in the host, leads to activation of phage WO [43]. Similarly, physiological changes in A. aegypti cells brought on by deleterious genomic changes could theoretically have cascading impacts on phage WO activity, though cautiously there are other potential reasons for these associations as well.
We quantified the potential deleterious effect of all protein coding SNPs and associated this with the population frequency to measure predicted deleterious load. Across all groups (High and Low blocking), we observed a significantly greater predicted mean deleterious load on chromosome 1 relative to chromosomes 2 and 3 (ANOVA, F = 20.62, d.f. = 2, P = 6.7x10 -5 ) (Fig 4A-left panel). Additionally, Low DENV blocking lines on average displayed a significantly greater predicted mean deleterious mutation load than High DENV blockers (ANOVA, F = 8, d.f. = 1, P = 0.01) (Fig 4A-right panel). Further analysis revealed that this increase in the predicted mean deleterious load is primarily driven by higher mean allele frequency of pre-  4C-right panel). From the original study, chromosome 1 is where the bulk of SNPs associating with differences between the High and Low blocking lines occurred [40].

cifA and cifB expression pattern differs between somatic and reproductive tissues in wMel-infected A. aegypti regardless of selection regime
Two key genes in the phage WO genome responsible for induction and rescue of CI are cifA and cifB. To test if changes in phage WO densities in the selected lines associate with phageinduced CI functions, and to examine a potential link between CI and WMPB, we performed whole-body gene expression analyses of the cifA and cifB genes in the selected lines. We To explore if these findings extended to specific tissues, we analyzed ovaries because they are a key site for expression of cifA-mediated rescue of CI [20] and transmission of Wolbachia [44]. There were several significant differences (Kruskal-Wallis H = 33.8, d.f. = 3, P < 0.0001). cifB expression levels were higher than cifA in both selection groups (High: Mann-Whitney U = 138.5, d.f. = 1, P = 0.0016; [1.9 fold], and Low: Mann-Whitney U = 85.5, d.f. = 1, P < 0.0001; [1.9 fold]), and the expression levels for both genes were significantly higher in the High blocking group (Fig 5, cifA, Mann-Whitney U = 144, d.f. = 1, P< 0.004; cifB, Mann-Whitney U = 151, d.f. = 1, P = 0.007). Taken together across the two body samplings, these results demonstrate that the transcriptional regulation of the CI genes likely responded to selection for WMPB in a tissue-specific manner.

Discussion
Wolbachia are currently at the forefront of vector control efforts to reduce transmission of arboviruses like dengue, Zika and Yellow fever to humans [12,13,45]. Successful field deployment of Wolbachia-infected mosquitoes in a population replacement strategy heavily relies on the long-term stability of both WMPB and CI. Previously, Ford et al., [40] demonstrated the potential for natural selection to maintain WMPB given a positive association with faster population growth. Here, aided by samples from this same experimental evolution study, we explored whether there is any evidence that temperate phage WO [26,46,47] play a role in Wolbachia-mediated DENV blocking.
We observed no evidence that phage WO relative densities in the somatic tissues predicted DENV load across selected lines, suggesting no direct link between phage load and WMPB. Our results did show both (i) a negative correlation between Wolbachia and phage WO relative densities when evolved populations were combined (Fig 2) and (ii) a positive correlation between Wolbachia and phage WO absolute abundance (Fig 1A and 1B), which is in line with Low DENV blocker mosquitoes have a higher predicted deleterious mutation load than High DENV blockers; an effect primarily driven by higher allele frequency of potential deleterious SNPs on chromosome 1 in the Low blocking line. Boxplots represent screened wMel-infected A. aegypti mosquitoes from WMPB selected lines for differences in (A) predicted mean deleterious load, (B) mean allele frequency of potential deleterious SNPs and (C) mean inverted SIFT score, as a proxy for accumulation of more severe SNPs across groups (High: dark blue-Low: yellow) and chromosomes (red). Significant comparisons expressed via P-values were determined via Type 2 and type 3 ANOVA and are listed in the results section of the main text. Each dot corresponds to a treatment group and its corresponded independent replicate line. N = 6 samples per chromosome and N = 3 samples per group for treatments in all graphs. https://doi.org/10.1371/journal.pntd.0009637.g004

Fig 5. Expression of CI-inducing genes differs in reproductive but not somatic tissues of wMel-infected A. aegypti from the 4 th generation of selection. Gene expression in whole body (minus ovaries) and ovaries of wMel-infected A.
aegypti mosquitoes selected for WMPB (High: dark blue-Low: yellow). Expression levels of cifA and cifB were normalized to Wolbachia groEL gene. Scatter plot depicts median levels of both genes for both the whole body and ovaries. Each dot represents a single mosquito. Sample sizes are depicted within the figure, in parentheses above each group, for the whole body minus ovaries and ovaries datasets. Significant comparisons expressed via P-values were determined via Kruskal-Wallis followed by post-hoc Dunn's multiple comparison test.
As for Wolbachia, the lack of a significant correlation between its densities and pathogen blocking (Fig 3) is inconsistent with what has been described in Drosophila [48,49]. However, other experiments involving A. aegypti have also failed to detect a significant correlation between WMPB and bacterial density [15]. This lack of correlation is part of an increasing body of evidence suggesting that WMPB may arise as a collective result of many distinct mechanisms [37,38], not simply density.
Surprisingly, the selection regime drove differences in phage WO density between High and Low blocking lines, but not in Wolbachia densities. As our previous analysis of the selected lines [40] also did not detect significant SNPs in the Wolbachia genome, the change in phage WO relative density is not likely due to Wolbachia. We propose that the change in WO relative density between High and Low lines may result from SNPs in mosquito loci, resulting directly from selection from linkage with those under positive selection. Since our list of SNP candidates associated with changes in blocking come from Low vs High comparisons, it cannot be determined in which lineages the changes occurred. Any of the selected SNPs may also have pleiotropic effects on Wolbachia physiology and/or phage stability. Similarly, the selection regime may also have resulted in hitchhiking of mildly predicted deleterious alleles in linkage to key blocking SNPs. Co-authors of this study have previously reported that Low DENV blocking mosquitoes displayed a significantly lower average fitness when compared to High DENV blockers [40]. Our bioinformatic analyses here identified a higher predicted mean deleterious load in Low DENV blocking mosquitoes versus the High blocking line (Fig 4A). A significant number of genes on chromosome 1 [40], near the SNPs selected upon (Fig 4B), are involved in processes such as cell cycle, transcriptional pausing, oxidative stress and others [41]. The increase in genetic changes in chromosome 1 could have disrupted Wolbachia's homeostasis, interfering with the balance between lytic and lysogenic cycles of the phage WO population. Phages are known to "sense" and respond to disturbances in the intracellular environment of the host where they reside [43]. Phage WO has numerous ankyrin repeat containing genes as well as other candidates, that appear to allow for direct protein-protein interaction with the host and that may provide a host sensing function [47,[50][51][52]. For instance, The cell cycle has been directly associated with bacteriophage productivity in Escherichia coli, in which infections of larger cells late in the cell cycle resulted in a significant increase in phage numbers [53], with consequences for bacterial populations and host symbiosis.
The reduced phage density in the Low blocking line also corresponded, not surprisingly, to lower expression of both CI-inducing genes cifA and cifB in the ovaries, relative to the High blocking line. This is, to the best of our knowledge, the first time in which the expression of the CI genes has been quantitated in this tissue for mosquitoes. The observation of higher ovarian expression of cifB relative to cifA is surprising but plausible. Lindsey et al., observed higher cifB mRNA levels relative to cifA during the first two-thirds stages of embryogenesis, as well as during larval stages [54]. Work in D. melanogaster has shown that cifA expression levels are dramatically higher than cifB when analyzing 1, 3 or 7 days old adult individuals [21,54]. This observation also seems to be consistent with reports in 11-14 day old ovaries of Drosophila pandora infected with wPanCI and wPanMK [55], and with the observation here for whole body samples. Lindsey et al., detected the existence of an intergenic region between cifA and cifB in wMel [54] that contains a hairpin termination signal that may impair the co-transcription of cifA and cifB. As such, the contrasting shift in expression pattern between these two genes could be attributed to yet unknown factors present in the ovaries that would regulate or interfere with such an intergenic regulatory region.
In conclusion, we propose that while WMPB and Phage WO load/CI-causing gene expression seem largely independent, there are potential scenarios involving pleiotropy or chromosomal linkage where selection on one trait may drive change in the other. Such relationships could lead to unforeseen consequences for Wolbachia release strategies. Future work should aim to understand the mechanistic nature of the blocking phenotype and how natural selection in the field may act upon the phenotype. Furthermore, the genomic region on chromosome 1 that is key to selected differences in blocking in our studies, should also be examined for potential roles in mosquito fitness and Wolbachia/phage WO homeostasis. Finally, without live male mosquitoes from the Ford et al. study we were unable to measure the strength of CI. Lower cif gene expression, however, would predict reduced CI in the Low lines, so future studies should assess the relationship between phenotypic shifts in blocking under selection and the expression of CI.

Mosquitoes
We used 5-7 day old, mated and blood-fed female A. aegypti mosquito samples (whole body minus ovaries; ovaries collected separately) from a previous study [40]. Briefly, in that study, mosquitoes from an initial population (Ancestral) were exposed to an artificial selection regime for Low (L) and High (H) levels of Wolbachia-mediated dengue virus (serotype 3) blocking (WMPB), with a randomly passaged control population (R) that did not undergo a selection regime. Each treatment group (L, R and H), with the exception of the original Ancestral population, was represented by three independent replicate lines. Ford et al [40] performed selection based on total body load of DENV (RNA copies) as determined by quantitative reverse transcription polymerase chain reactions (qRT-PCR) carried out over four generations. Upon completion of selection, in addition to DENV blocking, Wolbachia density and a suite of mosquito fitness traits were measured in the treatment lines. The ideal control for comparison to the H and L blocking lines is the co-reared R blocking line, as it was replicated and will have captured any drift and off-target selection resulting from the breeding regime. However, our previous study [40] demonstrated that given an association between faster population growth and higher blocking, that laboratory rearing of the R line also coselected for higher blocking. Since neither of the Ancestral or R lines are able to serve as true unselected controls, we have conservatively focused our comparisons between the evolved Low and High blocking lines, without attempting to interpret directionality of selection events.

RNA and DNA extractions
Samples were divided into two categories: 1) whole body minus ovaries and 2) ovaries. Despite such division, the ovaries used here were not matched to the same mosquito whose whole body was analyzed, as only a small fraction of paired samples were collected in the original study, and those were depleted during the progression of the experiments in the original manuscript by Ford et al. Samples were homogenized in 200μL of TRIzol reagent (ThermoFisher Scientific) containing a 2.8mm DNase/RNase-free ceramic bead (VWR) and macerated using a bead ruptor elite (Omni International). Total RNA and DNA were extracted according to the manufacturer's instructions. Nucleic acids were diluted in nuclease-free water and quantified using the NanoDrop 2000 spectrophotometer system (ThermoFisher Scientific). Each RNA sample was treated with DNase-I (Sigma) following the manufacturer's instruction to ensure no genomic DNA contamination.

Wolbachia, phage and mosquito gene quantification
Wolbachia density per mosquito host cell was obtained via quantitative polymerase chain reaction (qPCR) and quantifying the absolute copy numbers of the chaperonin groEL [see ref. [21] for the primer set] relative to the absolute copy numbers of A. aegypti ribosomal gene rps17 [see ref [40] for the primer set] [26,43]. Unfortunately, due to lack of a complete genome sequences for both integrated prophage and active phage particles, there is no current way of precisely distinguishing between lysogenic and lytic particles for WOMelB phage. As such, phage density was quantified in a similar manner to Wolbachia's relative quantification, using the absolute copy number of the phage capsid protein gene orf7 (specific to WOMelB) relative to absolute copy numbers of Wolbachia's groEL gene [43]. This quantification method potentially accounted, to a certain extent, to both integrated prophage sequences and putative lytic particles. Primers for orf7 (113bp) [orf7Fw: 5'-TGGCGAGAAAGCAGTAGAAATA-3// orf7Rv: 5'-GGTCATTTGTAGTTTTTTCATTCATAG-3'] were designed using NCBI's Primer-BLAST (www.ncbi.nlm.nih.gov/tools/primer-blast/) and both MFEPrimer 3.0 (mfeprimer.igenetech.com/) and IDT's OligoAnalyzer tool (www.idtdna.com) for quality control. Samples for quantifying the aforementioned parameters were run in technical duplicates. Prior to use in experiments, all primer pairs used in this study were examined for both specificity and amplification efficiency [56]. Specificity analysis was performed by melt curve analysis, with primer pair displaying a single peak. Efficiency analysis was achieved by examining the amplification performance under a series of sample template dilutions.
Standards for all genes, groEL, rps17 and orf7 were obtained from IDT as gBlock doublestranded DNA gene fragments resuspended in UltraPure distilled water (Invitrogen) to an equivalent of 10 12 copies of the genomic fragment. Standards were then aliquoted and frozen for one-time use. We used a total volume of 10μL per reaction, each containing: 5μL of 2x Per-feCTa SYBR Green SuperMix (Quantabio), 0.5μL of each forward and reverse primers (10μM), 2μL of nuclease-free water and 2uL of template DNA (100ng). Thermocycling conditions were as follows: initial denaturation at 95˚C for 30 sec, and 45 cycles of 95˚C for 3 sec and 60˚C for 30 sec, followed by melt curve analysis using a LightCycler 480 Real-Time PCR system (Roche).

Assessment of predicted deleterious load associated with SNPs in the Ae. aegypti genome
We have borrowed the concept of load from the field of population genetics, where load refers to the reduction in fitness of a population relative to population containing only individuals with perfect genotypes. More specifically, we have adapted the concept of mutation load to calculate what we refer to as the deleterious load [57]. The mutation load is measured as the negative selection coefficient, where larger values indicate that the mutation is more detrimental, multiplied by the allele frequency. Subsequently, the mutation load can increase either because the negative selection coefficient increases, or because the mutation increases in frequency within a population. To calculate the predicted deleterious load carried by each evolved replicate population, we substituted the negative selection coefficient used in mutation load with a measure of how deleterious a single nucleotide polymorphism (SNP) is and multiplied this by the population frequency of the SNP. The deleterious score of each SNP was acquired using the software SIFT (Sorting Intolerant From Tolerant) [58]. SIFT uses sequence homology across species to predict whether an amino acid or nucleotide substitution will affect protein function and hence, potentially alter the phenotype. The SIFT score ranges from zero to one and SNPs are by default categorized as deleterious if the score is less than 0.05.
Obtaining SIFT scores requires comparing SNPs of interest to a SIFT database. In order to build the SIFT database, we used the A. aegypti LVP_AGWG strain genome (AaegL5) [59] and the UniProtKB/Swiss-Prot database [https://www.uniprot.org/] for the comparison to other species. Overall, of the 131,171,571 protein-coding positions in the A. aegypti genome, 115,419,545 (88%) received a SIFT score, with 64,101,512 of these (56%) receiving confident scores. Low confidence scores occur when the protein-coding sequences used for prediction are closely related, and many positions appear conserved, which can lead to false positives. In the proceeding analysis, we used all positions with SIFT scores irrespective of confidence. Once the SIFT database was built, SIFT scores for SNPs segregating in each evolved treatment/ replicate were obtained via comparison of a Variant Call Format (VCF) file to the database.
In the previous analysis of this data, which identified cadherin as a candidate for resistance to DENV [40], SNPs were analyzed via a SYNC file [60,61] rather than a VCF file. Because a SYNC file cannot be readily converted to a VCF file, we had to recall SNPs. To call SNPs and produce a VCF file for this analysis, we used FreeBayes (Release 32) [62]. FreeBayes was run on each replicate with conservative settings to maximize the number of SNPs called (minimum alternate count of 1, minimum coverage of 5, maximum coverage of 200, and minimum base quality of 20). FreeBayes calls several different variant types including SNPs, multi-nucleotide polymorphisms (MNPs), insertions/deletions (INDELs), and complex haplotypes that can contain a combination of SNPs, MNPs, and INDELs. To decompose MNPs and complex haplotypes into SNPs, we used the vcfallelicprimitives and vcfbreakmulti functions from vcflib (github.com/vcflib/vcflib#vcflib). In order to produce the expected exponential distribution of alternate allele frequency in the Ancestral population, SNPs called by FreeBayes were further filtered to a minimum coverage of 20 and quality greater than 20 for the alternate call. Although each replicate's sequencing coverage was similar, the number of SNPs called after quality filtering ranged from~8M to~4M. Reproduction of a similar Manhattan plot as Fig  4A in [40] confirmed that SNP frequencies called by FreeBayes were similar to those previously called.
To confirm that there was biological signal in the SIFT scores, we compared the frequency of SNPs classified as deleterious to those classified as tolerated in the Ancestral population. Under the probable assumption that the Ancestral population was at mutation-selection-drift balance, deleterious loci should be segregating at a lower allele frequency than tolerated loci. This was indeed the case (GLM: p = 9.2X10 -25 ). To compare the predicted deleterious load of the High and Low DENV blocking treatments, load in each replicate was calculated using only those classified as potentially deleterious (SIFT score < = 0.05); load ¼ P n i¼1 ð1 À SIFT i Þ � frequency i n ; Where SIFT i is the SIFT score for a SNP, frequency i is the allele frequency of the SNP, and n is the total number of potential deleterious SNPs. Differences between selection treatments and chromosomes were then assessed via analysis of variance using the following linear model; load ¼ treatment þ chromosome þ treatment � chromosome þ error: Because the interaction was not significant (P = 0.94), the main effects of treatment and chromosome were assessed via type II sums of squares.

Expression of cytoplasmic incompatibility factor genes cifA and cifB
Detection of the number of the CI-causing transcripts cifA (locus tag: WD_RS02835 -known as WD0631) and cifB (locus tag: WD_RS06940 -known as WD0632), was performed by utilizing previously published primer sets [54]. As a reference for transcription activity of Wolbachia's core genome, we utilized the Wolbachia groEL gene. Machine-generated Cp values were used to calculate the relative expression of cifA:cifB using the ΔCt (Livak) method. We set up qRT-PCRs using a LightCycler 480 Real-Time PCR system (Roche) and a total volume of 10μL per reaction, each consisting of 5μL of 2x qScript One-Step RT-qPCR SYBR Green Master Mix (Quantabio), 0.2μL of qScript Reverse Transcriptase enzyme (RT) (50x) (Quantabio), 0.2μL of each forward and reverse primers (10μM), 2.4μL of nuclease-free water and 2uL of template RNA (100ng). Thermocycling conditions were as follows: initial cDNA synthesis step at 50˚C for 5min, RT deactivation/denaturation at 95˚C for 2 min, and 45 cycles of 95˚C for 3 sec and 60˚C for 30 sec, followed by melt curve analysis.

Statistical analyses
Normality for all Wolbachia and phage WO loads, as well as cif gene expression data were accessed using the Shapiro-Wilk test. Pairwise comparisons between non-normally distributed data were analyzed via Mann-Whitney's U test, while multiple comparisons between non-normally distributed data were analyzed using Kruskal-Wallis followed by posthoc Dunn's multiple comparison test with significance set at P < 0.05. Outliers were identified and removed prior to statistical comparisons, using 1% false discovery rate analyses, hence accounting for differences in sample size between datasets. Normally distributed data were compared with the Welch's ANOVA followed by Dunnett's T3 post hoc test with significance set at P < 0.05. Spearman's rank correlation coefficient analysis was performed in order to describe the relationship between Wolbachia and phage density across distinct generations and selection groups. The DENV-3 x phage total abundance interaction was analyzed via Spearman's rank correlation coefficient. All statistical analyses were performed using JMP Pro 15 (SAS) and Prism 8 (Graphpad).