Release from sexual selection leads to rapid genome-wide evolution in Aedes aegypti

The yellow fever mosquito, Aedes aegypti, mates in flight as part of ephemeral aggregations termed swarms. Swarms contain many more males than females, and males are thought to be subject to intense sexual selection.1,2 However, which male traits are involved in mating success and the genetic basis of these traits remains unclear. We used an experimental evolution approach to measure genome-wide responses of Ae. aegypti evolved in the presence and absence of sexual selection. These data revealed for the first time how sexual selection shapes the genome of this important species. We found that populations evolved under sexual selection retained greater genetic similarity to the ancestral population and a higher effective population size than populations evolving without sexual selection. When we compared evolutionary regimes, we found that genes associated with chemosensation responded rapidly to the elimination of sexual selection. Knockdown of one high-confidence candidate gene identified in our analysis significantly decreased male insemination success, further suggesting that genes related to male sensory perception are under sexual selection. Several mosquito control technologies involve the release of males from captive populations into the wild. For these interventions to work, a released male must compete against wild males to successfully inseminate a female. Our results suggest that maintaining the intensity of sexual selection in captive populations used in mass-releases is important for sustaining both male competitive ability and overall genetic similarity to field populations.


RESULTS AND DISCUSSION
Previously, we experimentally evolved replicate populations of Aedes aegypti in the presence or absence of sexual selection. The populations evolved in the presence of sexual selection experienced male-male competition and female choice, both of which are known to occur in the natural mating system. 2 In these populations, 5 males were held with a single female during a fixed mating period. In contrast, populations evolved in the absence of sexual selection were subjected to enforced randomized monogamy during the same period. We found that only 5 generations of evolution in the absence of sexual selection resulted in consistently decreased male competitive mating success. 3 Because the observed phenotypic changes indicated rapid evolution in response to the removal of sexual selection, we aimed to identify the signature of sexual selection in the genome using whole-genome pooled sequencing of these same populations to examine allele frequencies. We analyzed both patterns of genetic diversity in the two mating systems and the extent to which evolutionary change occurred in parallel across replicate experimental mosquito populations, allowing the identification, annotation, and functional validation of candidate SNPs underlying adaptative divergence.
Genomic variation is shaped by the presence or absence of sexual selection We sampled adult mosquitoes from replicate populations, which had evolved either with sexual selection (+SS: n = 3) or without sexual selection (ÀSS: n = 3), and the shared ancestral population, which was recently derived from the field (ANC: n = 1). We sequenced pools of 80 males from each population to a coverage of between 503 and 5003 using Illumina 150 bp paired-end reads. Sequences were aligned to the Ae. aegypti reference genome. 4 After stringent filtering, we generated a dataset of $8.7 million SNPs from the seven populations, a large but unsurprising number given the 1.25 Gb genome size (for full details of the sample preparation and bioinformatic pipeline, see the STAR Methods).
We found that sexual selection had a major effect on patterns of allele frequencies in evolved Ae. aegypti populations (Figure 1). Principal component analysis (PCA) of all SNPs indicated that populations evolved in the presence of sexual selection remained more similar on the first two principal component axes to the ancestral population than populations where sexual selection was eliminated ( Figure 1A; independent-samples t test on Euclidean distance from ancestral; t = À5.63, df = 4, p = 0.005). Although these first two axes revealed movement of the populations released from sexual selection away from the ancestral population in multivariate space, the third principal component axis consistently separated populations evolving with sexual selection from those where sexual selection was eliminated ( Figure 1B; independent-samples t test on PC3; t = 4.20, df = 4, p = 0.01). This clear separation by regime on PC3, which explains 18.6% of the total variance, indicates some degree of parallel adaptation in response to the mating system manipulation after only five generations of evolution.
We compared genetic differentiation between populations by calculating pairwise F ST between all evolved populations, as well as the ancestral population and the evolved populations, for non-overlapping windows of 500 kb across the genome ( Figure 1C). This analysis recapitulated the results of the PCA, where populations experiencing sexual selection had relatively low differentiation from one another (mean F ST = 0.023 ± 0.0005 SE), while populations evolving without sexual selection showed higher divergence from one another (mean F ST = 0.043 ± 0.001 SE, independent-samples t test on +SS versus +-SS and ÀSS versus ÀSS; t = À7.44, df = 4, p = 0.002). Genetic differentiation arises due to the cumulative action of both selection and drift. 5 Individual populations not experiencing sexual selection appeared more genetically similar to populations from the other evolutionary regime than they did to one another in both the PCA and our analysis of genome-wide F ST (independent-samples t test on +SS versus ÀSS and ÀSS versus ÀSS; t = À4.50, df = 10, p = 0.001). There was also greater differentiation between the ancestral population and populations evolved in the absence of sexual selection than there was between the ancestral population and populations evolved in the presence of sexual selection (independent-samples t test on ANC versus +SS and ANC versus ÀSS: t = À8.29, df = 4, p = 0.001) ( Figure 1C). These results suggest elevated genetic drift after the elimination of male-male competition and female choice and provide an explanation for the high levels of genetic differentiation among these populations. If the relaxation of (A) Populations exposed to sexual selection remained more similar to the ancestral population than populations not exposed to sexual selection on the first two principal component axes. (B) Evolutionary regimes are separated by PCA3. In both PCAs, all 8,704,928 SNPs were included (minimum coverage of 50 and maximum coverage of 500 in each population). (C) F ST comparisons (500 kb non-overlapping windows) between all replicate populations, grouped by within-and between-regime comparisons. Populations retaining sexual selection (+SS) showed the least divergence from one another, while across evolved regime F ST values were significantly higher. Populations evolving without sexual selection (ÀSS) showed the highest divergence from one another, indicating a prominent role of genetic drift over the course of experimental evolution. Between the evolved and ancestral populations, populations evolving with sexual selection showed the least divergence from the ancestral population. Colored points indicate pairwise means for individual populations, black points indicate overall mean, error bars indicate ±1 SE. Lower case letters a-e indicate significant differences in overall mean with p < 0.05. (D) Populations with sexual selection (+SS) maintained a greater estimated effective population size (N e ) than populations without sexual selection (ÀSS). Colored points indicate N e for individual populations, black points indicate overall mean and error bars indicate ±1 SE. Lower case letters indicate significant differences in overall mean with p < 0.05. (E) Estimated marginal means of chromosome-wide nucleotide diversity (p) are lower in the ÀSS populations for each chromosome. p was calculated using sliding windows of 10 kb on sites where 40% of coverage was greater than 20 and less than 500. Error bars indicate ±1 SE; *p < 0.05, **p < 0.01, and ***p < 0.001. See also Figure S1 and Data S1. ll OPEN ACCESS sexual selection freed populations to evolve in unique or arbitrary directions, then we would expect a reduction in variance effective population size (N e ) because these estimates are based on the magnitude of allele frequency changes over time. We used temporal allele frequency data from the ancestral population and our six evolved populations to estimate N e . 6 We found that populations evolved in the absence of sexual selection indeed showed reduced N e (mean = 23.97 ± 1.29 SE) relative to those for which sexual selection was present (mean = 52.84 ± 2.35 SE; Figure 1D; t = 10.77, df = 4, p < 0.001). This observed reduction in N e cannot be attributed to variation in offspring production or fewer breeding individuals of one sex 7 arising from our experimental design, because in both previous work 3 and in supporting experiments for this study (Data S1), there was no evidence that female mating probability or fecundity differed between selective regimes. Similarly, because Ae. aegypti females are largely monandrous, 2 and the time spent with males was limited for each generation, it is unlikely that evolutionary regimes differed in the number of contributing males. Our estimates of effective population sizes, combined with patterns of genetic differentiation between the evolved populations, strongly support the idea that evolutionary change occurred via elevated genetic drift after the relaxation of sexual selection. A less likely alternative explanation is that the removal of sexual selection caused populations to rapidly evolve toward new allele frequency optima, which would also result in a lower estimated N e . However, we do not see any evidence of large, consistent changes in these populations to support this alternative ( Figure 1A).
Consistent with the observed differences in N e between regimes, nucleotide diversity (p), calculated in non-overlapping 10 kb windows for each replicate population, also differed between evolutionary regimes. There was spatial variation across the genome in nucleotide diversity (chromosome effect: c 2 = 324.91, p < 0.001), but a consistent effect of reduced nucleotide variation in populations that evolved in the absence of sexual selection (regime effect: c 2 = 11.47, p < 0.001). Although significant on all three chromosomes, the reduction in p in ÀSS compared with +SS is most pronounced on the 2 nd and 3 rd chromosomes and smaller on the 1 st chromosome ( Figure 1E). Interestingly, p on chromosome 1 was higher than on chromosome 2 and 3 for both regimes. We found that this elevated p was unlikely to be attributable to the presence of the sex-determining M locus ( Figure S1), given that we observed locally reduced p in this region. Other studies reporting a relatively elevated p on chromosome 1 in Ae. aegypti 8 suggest that additional features of the chromosome, such as relatively low gene and transposable element density, or a high number of satellites, 9 might be creating this pattern.
Altogether, our analyses of genome-wide allele frequency data show that the removal of sexual selection both leads to the loss of genetic diversity and drives rapid evolutionary change away from the ancestor.

Consistent allele frequency differences between selection regimes
In addition to evidence for an important role of genetic drift in the evolutionary changes observed in our experimental populations, we detected a clear signal of differentiation between the regimes in our multivariate analyses ( Figure 1B). We identified candidate SNPs underlying differentiation between evolutionary regimes using generalized linear mixed models (GLMMs) with a false discovery rate (FDR) of 10%. This approach identified 50,988 SNPs spanning all three chromosomes (Figure 2A), with large and consistent differences in allele frequency between evolutionary regimes that represent either the targets of selection or linked genetic variation.
In order to determine to what extent reduced opportunities for recombination had inflated our number of candidate SNPs, we calculated Pearson's correlation coefficient (r w ) between residuals of all pairs of candidate SNPs on the same chromosome (STAR Methods). Because the distribution of r w for pairs of candidate SNPs on the same chromosome and the distribution of r w for pairs of candidate SNPs on different chromosomes were similar (Figure S2B), we assessed the level of linkage to be low ( Figures S2B  and S2C). Under stringent clustering criteria (candidate SNPs were grouped into clusters if they were located within 100 kb of the midpoint between one another and had highly correlated allele frequencies of r w > 0.8), we found 2,161 clusters ( Figure S2A). Candidate SNPs assigned to different clusters appear to be no more correlated to one another than they are to candidate SNPs on different chromosomes, appearing effectively independent ( Figures S2B and S2C). Further, all gene-level candidates in our functional experiments (see below) were in independent clusters, eliminating the likelihood that linkage disequilibrium (LD) caused these genes to evolve together, and consistent with the likely polygenic architecture of quantitative traits. This pattern is consistent with typical ''islands'' of divergence seen in comparable experimental evolution studies in Drosophila pseudoobscura, 10 though with predictably lower resolution that might be expected given fewer recombination events in five generations.
A functional annotation of the 50,988 candidate SNPs identified a subset with a moderate (n = 713) or high (n = 12) impact on amino acid sequence. Allele frequency trajectories of high-impact SNPs ( Figure 2B) show a large (mean 19.4%) and consistent average allele frequency difference between evolutionary regimes. A gene ontology (GO) enrichment analysis on high-and moderate-impact SNPs (n = 725) indicated that they colocalized with genes enriched for GO terms related to olfactory receptor activity, sensory perception of smell, and other chemosensory terms (Data S2A and S2B), suggesting that odor perception may be a target of sexual selection. We additionally ran this analysis to include non-coding regions 2,000 bp either side of genes and recovered the same GO terms (Data S2C). Among our list of moderate-and high-impact genes were five odorant receptor genes (or53 [aael015286], or24 [aael017557], or44 [aael006465], or116 [aael025139], and or63 [aael000628]). 4 Sensory systems underpin key components of mosquito mating behavior 11,12 and recent work suggests that chemical cues may play a role in mating swarm formation and mating success. 13,14 Alongside this previous experimental work, our results implicating sensory pathways as rapidly evolving targets of sexual selection make a compelling case for further investigation into the role of chemosensation in mating behavior.
Depleting ppk315 expression reduced male insemination capacity Eleven genes colocalized with candidate high-impact (null) variants. Using previously published expression data from Ae. aegypti, 15 we determined that seven of the 11 genes showed high expression in adult male tissues. We targeted these seven genes for RNAi-mediated gene silencing in order to test their role in male mating-related phenotypes. Four of the seven genes were effectively silenced 2 days following intrathoracic injection with double-stranded RNA (dsRNA) (gustatory receptor11, pick-pocket315, AAEL022941, and AAEL005908) and were available to screen for effects on male mating success ( Figure S3A).
To measure insemination capacity of males with reduced gene expression, 2 days following dsRNA injection, virgin male mosquitoes were held singly with six virgin females and allowed to mate for 24 h. The reproductive tracts of females were dissected to determine the number of females that males were able to successfully inseminate in this period. Males with reduced pickpocket315 (ppk315) expression inseminated significantly fewer females relative to the dsGFP control ( Figure 3A; GLM: p < 0.001). Strikingly, while 89% of control males were able to mate in the 24-h period, only 55% of the males with depleted ppk315 expression were able to achieve a single successful mating ( Figure 3B; GLM: p < 0.001). This was not a result of a general deleterious effect of ppk315 knockdown in these males, as we did not observe any significant decrease in survival associated with the knockdown of ppk315 expression in the first 7 days ( Figure S3B). Knockdown of gustatory receptor11, AAEL022941, and AAEL005908 did not affect insemination capacity (GLM with replicate as factor: p = 0.38, p = 0.70, and p = 0.13, respectively).
The insemination capacity assay used to assess the effects of knockdown here required that males only be able to locate and inseminate females in small cups in the absence of competing males. Thus, these results do not necessarily indicate that the remaining candidate genes have no effect on male mating success (and were therefore false positives in our screen). It may be that these genes instead have effects that were not captured by our assay. Future work will need to use more sensitive assays, incorporating male competition and behavioral observation to detect more subtle effects on male mating phenotypes. Conversely, it is possible that off-target effects may have altered mRNA levels of genes elsewhere in the genome. Subsequent work using CRISPR-based methods could be used to silence our high-confidence candidates. Finally, although our focus was on male mating phenotypes, both sexes were under selection in these regimes. Although previous work did not identify changes in female mating responses, further work on how these candidates affect female biology are warranted.
There is evidence in the literature to support the role of pickpocket (ppk) genes in mating behaviors. Indeed, a ppk gene was identified as a potential candidate in an evolve and resequence experiment manipulating the intensity of sexual selection in Drosophila. 10 ppk genes belong to a superfamily of degenerin/ ENaC channel genes that encode cation channel subunit receptors of a wide variety of stimuli and have been shown to be involved in the perception of female pheromones. 16,17 Ae.

OPEN ACCESS
aegypti has 31 ppk genes, and ppk301 has been shown to be important for the detection of salt in Ae. aegypti. 4,18 Transcriptomic data show expression of ppk315 in the rostrum (composed of the maxillary palps and the proboscis) and the terminal abdominal segment of male Ae. aegypti mosquitoes. 15 Our functional analysis ( Figure 3) and genetic analysis ( Figure 2) confirm that ppk315 is indeed under sexual selection and suggest a full behavioral and molecular characterization of this gene is needed.

Conclusions
It is increasingly recognized that male mosquitoes experience intense sexual selection due to the presence of male mate competition and female choice. 2 The majority of research on traits involved in male mosquito mating success has focused on the contribution of environmental factors (reviewed in Cator et al. 2 ). Here, we provide evidence, as suggested by previous phenotypic observations, 3 that mating success also has genetic determinants which can rapidly evolve.
Many investigations of sexual selection in mosquitoes have also been observational and have focused on correlations between traits and mating success. 2 In addition to the limits on inference that are imposed by work relying on correlations, observational studies are also challenged by the fact that they rely on researchers to have some idea a priori which traits might be targeted by sexual selection and should therefore be measured. The use of evolve and resequence to investigate sexual selection in mosquitoes is a powerful complement to this kind of work because the genetic basis of mating success can instead be identified by the cumulative and parallel action of selection. In addition, this approach allows estimation of the net effect of sexual selection on patterns of genetic variation.
This technique shows great potential for providing novel insights into control relevant aspects of mosquito behavior and evolution. As in other studies manipulating sexual selection, 10,19-21 there was a limited period during the life cycle where social dynamics differed between our opposing sexual selection regimes. There is evidence from Drosophila suggesting that altered social environments might favor different behaviors (e.g., longer mating duration in response to the perception of more intense mating competition 22,23 ), and the potential effects of altered social experience should be considered in the design of future experiments.
In our design, populations evolving with sexual selection that retained competitive sexual performance 3 were more similar genetically to the wild populations from which they originated and harbored more genetic variation. This suggests that the selective environment, including intensity of sexual selection, experienced by lab-reared populations destined for use in reproductive control could have significant impacts on male mating success and control efficacy.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:

Lead contact
Further information and requests for resources should be directed to and will be fulfilled by the lead contact, Lauren Cator (l.cator@ imperial.ac.uk).

Materials availability
This study did not generate new unique reagents.

Data and code availability
Whole-genome sequence data have been deposited at the European Nucleotide Archive (ENA) and are publicly available as of the date of publication. Accession numbers are listed in the key resources table. All data analyzed in this study are available as supplemental files. This paper does not report original code. Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

EXPERIMENTAL MODEL AND SUBJECT DETAILS
Experimental populations and selection regime Mosquito populations were established as described in Qureshi et al. 3 Briefly, around 4500 Aedes aegypti eggs were collected from 17 water storage containers from two villages in Muang District, Kamphaeng Phet Province (KPP), Thailand between February and April 2016. Eggs were hatched under a vacuum for 20 minutes. Newly emerged larvae were provided with 0.1 mg of ground fish food (Cichlid Gold [#04328], Hikari, Himeji City, Japan) and held in a 27 C incubator overnight. First instar larvae were sorted into trays of 500 larvae/1L water and provided fish food pellets ad libitum. The initial egg collection was divided into six experimental populations, with three populations assigned to a sexual selection mating regime (+SS) and three populations assigned to a no sexual selection mating regime (-SS). In each of the +SS cages, five males competed to inseminate a single female (allowing both male-male competition and female choice). In -SS cages, one male was placed with one female, eliminating any opportunity for male-male competition or mate choice. One hundred mated females were pooled for each population (300 females per mating regime), offered a bloodmeal (FirstLink, UK) using a Hemotek feeder, and allowed to lay eggs on damp filter paper. Eggs were again hatched under a vacuum, larvae fed ground fish food, then sorted at first larval instar and fed ad libitum until adulthood. This process was repeated for 5 generations, after which the populations underwent a single generation of common garden rearing. For more details on the experimental evolution set-up see Qureshi et al. 3

METHOD DETAILS
Genomic DNA extraction Seven different samples were collected from the experimental evolution set-up: a single sample from the ancestral population, a sample from each of the three replicate populations of the evolved sexual selection regime, and a sample from each of the three replicate populations of the evolved no sexual selection regime. Samples consisted of pools of 80 male mosquitoes. Genomic DNA (gDNA) was extracted from each of the seven samples using the QIAGEN DNeasy Blood & Tissue kit according to the QIAGEN protocol for purification of total DNA from insects. 39 In order to ensure approximately equal amounts of gDNA were obtained from each individual male, pools for each replicate population were subsampled in groups of 20, gDNA extracted and then equal quantities of gDNA combined from each subset.
DNA sequencing, aligning, and filtering Sequencing libraries were prepared with TruSeq DNA Nano kit (Illumina) and sequenced on the NovaSeq 6000 platform (CeGaT GmbH, Germany) with 150-base pair paired-end reads. FastQ files were quality checked using FastQC version 0.11.9 with default settings. 40 Reads were mapped to the Ae. aegypti genome (Liverpool AGWG-AaegL5.0) 4 using BWA-mem (version 0.7.17). 31 PCR duplicates were removed using Samtools version 1.3.1 32 and reads realigned around intervals using PicardTools (version 2.6.0) and GATK version (3.4). Realigned sequence files were compiled using the 'mpileup' command in Samtools and then converted to sync file format using the 'mpileup2sync' command in Java multi-threading from the package Popoolation2. 35 Allele frequency differences (_rc file) were called using the 'snp-frequency-diff' command in Popoolation2 based on a minimum coverage of 50, maximum coverage of 500 and minimum minor allele count of 14. 35 A custom R script removed monoallelic and triallelic loci. The mitochondrial genome was excluded from the analysis due to excessively high coverage. We also removed sites with a minor allele frequency < 5 %. After these filtering steps, the remaining 8,704,928 SNPs were used for downstream analysis.

Genome-wide patterns of variation and differentiation
Principal component analysis was run on all filtered SNPs (>8.7 million) using the R function 'prcomp'. Next, we calculated the fixation index (F ST ) to determine genetic differentiation between pairs of evolved populations. F ST is the proportion of genetic variance in a subpopulation relative to the total variance. High F ST values are indicative of a large amount of genetic differentiation between subpopulations. We first converted filtered BAM files used in earlier analyses to 'mpileup' format using Samtools, 32 then used a script from Popoolation2 to convert the 'mpileup' to a 'sync' format file. We used PoPoolation2 34 to average F ST values for non-overlapping windows of 500kb across the three chromosomes, using parameters of a minimum allele count of 14, minimum coverage of 50 and maximum coverage 500 reads. Population size (N e ), which greatly affects changes in allele frequency, was estimated for each evolved population using a method designed for pool-seq data which corrects for sampling at both the pooling and sequencing stages. 6 Nucleotide diversity (p) is the average pairwise difference between all samples, and a measure of overall genetic variation between populations. Low nucleotide diversity may indicate selection at a locus. p was calculated for each of the three replicates of the two evolved populations. BAM files for each sample were combined using the Samtools 'mpileup' command, 32 subsampled to a uniform coverage of minimum 50 and maximum 500 reads, then p was calculated in 10kb non-overlapping windows using PoPoolation. For any given window to be assigned a p value, at least 40 % of the window had to comply with a minimum and maximum coverage of 20 and 500, respectively. 34 We fit a linear mixed effect model (lmer) using the R package lme4 27 to determine if there was a difference in p between evolutionary regimes genome-wide and for individual chromosomes, incorporating population and the position of the 10 kb window in the genome into the model as random effects.
To identify SNPs that had diverged due to the selective regime, we employed a generalized linear mixed model (GLMM) to identify SNPs with consistent allele frequency differences between all 3 replicate populations of each evolutionary regime. The GLMM compared major allele frequencies of the 8,704,928 SNPs between the two evolutionary regimes, with evolutionary regime as the fixed effect and population as a random effect. The alternative and null models were then compared with a likelihood ratio test to obtain variant-level significance tests. We then calculated q-values (false discovery rate) using the method of Storey & Tibshirani and applied a threshold of 0.1 or 10 % using the qvalue package in R. 41 ''Significant SNPs'' were those which consistently diverged between treatments with a q value < 0.1.

Linkage disequilibrium estimations
Linkage disequilibrium (LD) is the non-random association of neighboring polymorphisms due to limited recombination. 42 Experimental populations were evolved for five generations and so there were limited opportunities for recombination to occur. We therefore might expect high levels of linkage due to hitchhiking of neutral variants along with any selected variants, and so an inflated number of candidate SNPs, as is fairly common with pooled sequencing data. 43 To determine the extent to which linkage disequilibrium caused our 50,988 diverging candidate SNPs to evolve non-randomly we calculated the number of independently evolving clusters on each chromosome, based on a method described in Kawecki et al. 44 using a custom R script and the packages dplyr 29 and data.table. 30 First, we calculated residual allele frequencies for each evolved replicate population by subtracting allele frequencies at each SNP locus from the regime mean. We then calculated the Pearson's correlation coefficient (r w ) between residuals of all pairs of SNPs on the chromosome. We expected that if linkage disequilibrium was high the distribution of r w between SNPs on the same chromosome would be more positively skewed than the distribution of r w between SNPs on different chromosomes. Conversely if the linkage was low these two distributions would appear similar. To group SNPs into clusters, we needed to determine a threshold correlation coefficient and maximum distance between SNPs that could be clustered together. To do this we calculated the number of clusters for a range of threshold correlation coefficient values (r w = 0.7, 0.8, 0.9) and distances between SNPs (100kb, 200kb, 500kb, 1000kb), and selected values for which the correlation between clusters on a single chromosomal arm was similar to the correlation between SNPs on different chromosomes. As such, SNPs were assigned to the same cluster if they were located within 100kb left or right of the midpoint between each other, and if they were found to be highly correlated (r w > 0.8).

Genetic variant effect prediction
To determine the functional effects of the significantly diverged SNPs we used the program SNPEff (version 4.3). 36 First a custom database was created with the most recently annotated Ae. aegypti genome. 4 A Variant Call Format (VCF) file was compiled for all 8.7 million SNPs (biallelic variants with maximum 500x coverage, minimum count 14 and global frequency > 5 %) and a subset of 50,988 significant SNPs (< 0.1 FDR) from the GLMM analysis. SNPeff requires input VCF files to contain a reference allele and an alternative allele for each variant. The reference was extracted from the reference genome using BEDTools 37 and the alternative allele as both alleles present at the SNP site. SNPeff generated a list of variants with either high or moderate effect on gene function; moderate effect substitutions resulted in a change of amino acid, and high effect substitutions resulted in a stop codon. To test whether moderate and high impact variants were in genes enriched for any functional groups, we performed a gene ontology analysis with Gowinda. 38 The program performed 100,000 simulations and counted every gene with a SNP in it only once.

Functional validation of high-impact variants
The 12 high impact variants were found in 11 candidate genes identified by SNPeff. The seven of these candidates that were differentially expressed in adult male Ae. aegypti tissues based on RNA-seq data 15 were selected for functional validation. Expression of candidate genes was disrupted using RNAi-mediated gene silencing by intrathoracic injection of dsRNA. The reference strain, Ae. aegypti Liverpool was used for functional validation of candidate genes as sequence fidelity was assumed to be higher.

dsRNA synthesis
To achieve gene knockdown total RNA was extracted from adult Ae. aegypti Liverpool mosquitoes using Trizol (ThermoFisher), and residual genomic DNA removed using TURBO DNA-free Kit (Life Technologies). Complementary DNA (cDNA) was synthesized using the LunaScript RT SuperMix (NEB). Target regions (200-600 bp) were amplified from cDNA by PCR using CloneAmp HiFi PCR Premix (Takara) and custom primers flanked with the T7 promotor sequence TAATACGACTCACTATAGGG (see key resources table and Data S3A). Ribosomal protein S7 primers (custom design) were used as an internal control to confirm cDNA presence. GFP served as a negative control (primers from Eggleston and Adelman 24 ) and was synthesized using a plasmid template containing GFP (gift from D. Ellis).
To generate sufficient template, three 50 mL PCR reactions were performed for each gene. Each 50 mL reaction consisted of 25 mL Hifi polymerase, 5 mL of each forward and reverse primers, 2.5 mL of cDNA and 12.5 mL of dH 2 0. Amplification was performed on a Veriti 96-well thermal cycler at the following conditions: 1 m at 98 C, followed by 40 cycles of 10 s at 98 C, 10 s at 56 C, and 30 s at 72 C, and lastly a 2 m extension time at 72 C. To verify the amplicon size, 5 mL of PCR products were run on a 2 % agarose gel at 90 V for 40 mins. The remaining 145 mL of the pooled PCR reactions were then purified using the QiaQuick PCR purification kit (QIAGEN) to remove any off-target amplification products. To generate dsRNA 1 mg of purified PCR product was used as the template for in vitro transcription using T7 high yield transcription kit (ThermoFisher). Two 20 mL reactions were carried out for each gene and incubated at 37 C overnight to achieve maximum dsRNA yield. Both in vitro transcription reactions were then pooled and purified using the RNeasy kit (QIAGEN), eluting in 20 mL of RNase-free water (Qiagen). The purified dsRNA was then concentrated using 3 M sodium acetate (pH 5.2) and 100 % ethanol then stored at -20 C prior to use. dsRNA concentration was determined by measuring a 10x diluted sample on a spectrophotometer and 2 mL run on a 1 % agarose gel at 90 V for 40 mins to check integrity. All samples were adjusted to 7.25 mg/mL to facilitate delivery of 500 ng dsRNA in a 69 nL volume injection (max capacity of injector). Immediately prior to microinjection the dsRNA was heated to 95 C for 5 mins and cooled back to 25 C at a rate of 0.1 C /s to anneal any single stranded RNAs (ssRNA) (G. P. League, personal communication).

Microinjection of dsRNA
Ae. aegypti Liverpool eggs were hatched under a vacuum for 20 minutes. Newly emerged larvae were provided with 0.1 mg of ground fish food (Cichlid Gold [#04328], Hikari, Himeji City, Japan) and held in a 27 C incubator overnight. First instar larvae were sorted into trays of 500 larvae / 1 L water and provided fish food pellets ad libitum. Pupae were sorted by sex and allowed to emerge in separate cages. dsRNA was injected into the thorax of cold-anesthetized 0-2 day old virgin male Ae. aegypti mosquitoes using the Nanoject II microinjector (Drummond Scientific) and borosilicate glass capillary needles (500 ng per mosquito). Mosquitoes were allowed two days to recover post injection and provided with 10 % sucrose ad libitum.
Quantitative real-time PCR (q-PCR) qPCR was performed to detect expression of target genes and evaluate knockdown efficiency. Primers were designed for each of the seven target genes using Benchling and a previously published S7 internal control 25 (see key resources table and Data S3B). Primer efficiencies were first calculated to compare expression of target and housekeeping genes. Primer efficiencies were calculated in Excel using the average Cq values and log 10 values for the cDNA dilution series. Primers all had an efficiency of 80-120 %. Meltcurves were also performed for all primer pairs to confirm the presence of just one amplicon.
cDNA template was synthesized from RNA extracted from a pool of 20 0-2 day old male mosquitoes as described above. Reaction mixtures totaling 10 mL were composed of 1 mL of male cDNA, 5 mL of PowerUp SYBR Green Master Mix (Life Technologies), 0.5 mL of each the forward and reverse primer (10 mM) and 3 mL of dH 2 0. Amplification was performed on the LightCycler 96 (Roche) using the standard cycling mode recommended for PowerUp SYBR Green Master Mix and primers with a T m R 60 C (50 C for 2 min, 95 C for 2 min, 40 cycles at 95 C for 15 sec and 60 C for 1 min, followed by a single cycle of melt curve 95 C for 10 sec, 65 C for 1 min and 97 C for 1 sec). Relative gene expression was calculated for each treatment relative to the GFP control based on 3 biological replicates and 2 technical replicates using the Pfaffl method. 45 A statistically significant reduction in gene expression was determined by comparing the gene expression ratios of the control and the treatment gene with and independent t-test.

Survival assay
To test whether knockdown of any of the genes of interest increased mortality, survival of injected male mosquitoes was recorded. Three replicates of 20 0-2 day old Ae. aegypti males were injected with dsRNA and held in 16 oz paper cups and provided with 10 % sucrose ad libitum. Mortality was recorded every other day from one day post-injection (dpi) until 30 days post-injection. A log-rank test was used to determine if there was a difference in survival between each test gene and the control pre-and post-seven days following dsRNA injection.

Insemination capacity
To determine if dsRNA injection influenced the ability of a male to mate a female or multiple females, insemination capacity was measured. dsRNA-injected males (2dpi) were held singly in 16 oz cups with six virgin females and provided 10 % sucrose solution ad libitum. After 24 hours all females were collected and their spermathecae dissected in 1 x phosphate buffered saline (PBS) and examined under a brightfield microscope to determine insemination status. The insemination capacity of around 60 males per gene were assessed over 3 replicates. Whether males were able to achieve a mating in 24 hours (Y/N) was analyzed with a binomial Generalized Linear Model (GLM) with treatment and replicate as predictors. Number of females mated out of a maximum of six in 24 hours was analyzed with GLM (Poisson distribution) with treatment and replicate as predictors.

Mating and fecundity measurements
To determine if the selective regimes themselves lead to changes in fecundity and insemination capacity of evolved populations, the proportion of females mating, and the number of eggs laid after mating, were measured for the -SS and the +SS treatments. We set up 50 pairings per treatment (-SS-1_:1\, +SS-5_:1\) in the same containers used in our experimental evolution design. 3 As in the selective regime, males were transferred into containers the night before the experiment and provided with 10 % sucrose solution. Females were introduced into containers at 0700 and moved to pooled cages 8 hours later. The following day females were provided with a bloodmeal and transferred into individual tubes for oviposition and provided with a 10 % sucrose solution. After six days, eggs were removed and counted. Females that failed to lay eggs were dissected to confirm insemination. This experiment was conducted twice. Insemination data was analyzed with a binomial with replicate and treatment as predictors and inseminated (Y/N) as the response variable, and number of eggs produced per females was analyzed with a GLM with replicate, wing length, and treatment as the predictors and eggs as response.

QUANTIFICATION AND STATISTICAL ANALYSIS
Quantification was performed as described in the relevant method details sections above. Statistical analysis was performed with t-tests for comparisons of Euclidean distances in the PCA, comparison of F ST , and N e . A linear mixed effect model was used to detect differences in p between evolutionary regimes both genome-wide and within individual chromosomes. GLMMs were used for comparison of allele frequencies between evolutionary regimes, with treatment as a fixed effect and population as a random effect, with significance determined using Storey & Tibshirani q values < 0.1. Comparison of gene expression ratio for control and treatment genes following dsRNA injection was also performed with t-tests. Log-rank tests were used to assess survival following dsRNA injection. Binomial GLMs were used to analyze insemination capacity measurements, as well as mating and fecundity measurements.