Rapid evolution of knockdown resistance haplotypes in response to pyrethroid selection in Aedes aegypti

Abstract This study describes the evolution of knockdown resistance (kdr) haplotypes in Aedes aegypti in response to pyrethroid insecticide use over the course of 18 years in Iquitos, Peru. Based on the duration and intensiveness of sampling (~10,000 samples), this is the most thorough study of kdr population genetics in Ae. aegypti to date within a city. We provide evidence for the direct connection between programmatic citywide pyrethroid spraying and the increase in frequency of specific kdr haplotypes by identifying two evolutionary events in the population. The relatively high selection coefficients, even under infrequent insecticide pressure, emphasize how quickly Ae. aegypti populations can evolve. In our examination of the literature on mosquitoes and other insect pests, we could find no cases where a pest evolved so quickly to so few exposures to low or nonresidual insecticide applications. The observed rapid increase in frequency of resistance alleles might have been aided by the incomplete dominance of resistance‐conferring alleles over corresponding susceptibility alleles. In addition to dramatic temporal shifts, spatial suppression experiments reveal that genetic heterogeneity existed not only at the citywide scale, but also on a very fine scale within the city.

(F1534C) (numbered according to homology in the house fly, Musca domestica) in the VGSC protein (Table 1) (Deming et al., 2016;Linss et al., 2014). Another common single nucleotide mutation in Brazil and Mexico that is associated with pyrethroid resistance is the valine (V) to leucine (L) amino acid substitution at locus 410 (Haddi et al., 2017;Saavedra-Rodriguez et al., 2018). Different nonsynonymous mutations in the VGSC gene are more common in other parts of the world. For example, the single nucleotide mutation resulting in the V1016G substitution is common in Asia (Li et al., 2015) and has only recently been reported in the Western Hemisphere (Murcia et al., 2019).
Because the mutations causing V1016I and F1534C are within the VGSC gene, they are expected to be in linkage disequilibrium when polymorphic (Hartl, 2000). Haplotypes that contain both the Cys1534 allele and the wild-type Val1016 allele are common in wild populations (Deming et al., 2016;Linss et al., 2014). In contrast, the Ile1016 allele is rarely found on the same chromosome as the wildtype Phe1534 allele (Linss et al., 2014). Vera-Maloof et al. (2015) suggest that the Cys1534 allele may be required to support the presence of the Ile1016 allele, yet these resistance alleles do not always rise in frequency together. In Brazil, the frequency of both resistance alleles, Ile1016 and Cys1534, increased over time at all sites sampled, but the ratio of different haplotypes of the gene with one or both nucleotide changes varied geographically (Linss et al., 2014). More recently, Ae. aegypti kdr haplotype analysis has confirmed that the Val1016/Cys1534 haplotype is found worldwide, while the Ile1016/Cys1534 haplotype is found only in the Americas (Cosme et al., 2020;Fan et al., 2020). This implies that both widely dispersed and geographically constrained haplotypes contribute to the widespread pyrethroid resistance observed in Ae. aegypti.
Although specific resistance haplotypes are commonly observed in Ae. aegypti in the Western Hemisphere, results from two studies in Merida, Mexico, indicate that heterogeneity in frequencies of common resistance haplotypes can exist at the scale of a city block (Deming et al., 2016). High spatial heterogeneity means that fine-scale spatial variability is not always observable when assessing insecticide resistance at city or region-wide scales (Grossman et al., 2019). These two studies are important for understanding the scale of geographic variability in Ae. aegypti insecticide resistance and the role that limited mobility of this species plays in this evolutionary process. These studies were not, however, able to determine whether the heterogeneity in resistance is due to spatial variation in insecticide pressure or to variation in initial resistance allele frequencies. A better understanding of fine-scale spatial structure of resistance genetics would help public health officials and scientist preserve the efficacy of important insecticides and aid in the mathematical modeling of new insect control technologies.
Iquitos, Peru, is a well-established study site for Ae. aegypti and dengue in the Western Hemisphere with a long history of entomological research and sampling that can help elucidate the factors contributing to pyrethroid resistance evolution in this mosquito (Cromwell TA B L E 1 The number of mosquitoes that were genotyped and used for analysis in this study. Total males analyzed are the number of males that were genotyped at both V1016I and F1534C loci and used in the haplotype analysis. Buffer and spray zones correspond to the subset of the total males analyzed that were collected in each area during the 2013 and 2014 suppression experiments. The discrepancy between the totals is due to mosquitoes collected outside of the study area in that year. The V410L genotyped column indicates the subset of total males that were genotyped at the 410 locus in each year  Gunning et al., 2018;LaCon et al., 2014;Morrison et al., 2010). Iquitos is located in the northeastern Peruvian Amazon and is only accessible by boat or plane, and there are no roads into Iquitos.
The population was estimated to be ~437,300 in 2015 (Nacional & de Estadistica e Informatica, 2012). Ae. aegypti were initially eradicated from Peru in 1958, but reinvaded and were detected again in 1984 (Phillips et al., 1992). In anticipation of dengue moving into areas with Ae. aegypti populations, the U.S. Naval Medical Research Unit No. 6 (NAMRU-6) established a field office in Iquitos in the late 1980s to monitor for Aedes-transmitted viruses. Because local DENV transmission has occurred since the early 1990s, repeated citywide insecticide spraying and long-term epidemiological monitoring efforts have been carried out to control Ae. aegypti populations in an effort to reduce disease (Morrison et al., 2010;Stoddard et al., 2014). Dengue-1 was detected in Iquitos in 1990 (Phillips et al., 1992), and dengue-2 was confirmed in 1995 .

Continuous hospital and outpatient clinical surveillance began in
1993 , and household vector surveillance began in 1998 . Household vector surveillance collection spans the years before, during, and after pyrethroid use in the city, making Iquitos an excellent location to explore the natural evolution of pyrethroid resistance in Ae. aegypti.
A priori, we anticipated that the two most commonly observed kdr mutations in Central and South America, V1016I and F1534C, would rise in frequency in response to the application of pyrethroids.
We also expected that existing heterogeneity in the population would impact the patterns of resistance haplotype frequencies both temporally throughout the city and spatially on a fine scale among city blocks. In addition, the dominance and selection coefficients of the resistance alleles were anticipated to affect the rate that resistance haplotypes increased in frequency in the population.

| Temporal collections
Aedes aegypti were collected and stored at −80°C by NAMRU-6 and University of California at Davis personnel since the late 1990s.
Specimens dating back to the year 2000 were available for study.
Mosquitoes were collected by backpack aspirator (Clark et al., 1994) prior to June 2009 and by Prokopack Aspirator following June 2009 (Reiner et al., 2019;Vazquez-Prokopec et al., 2009). Each mosquito in the repository was identified to species, sex, collection date, and collection site. Each collection site, typically an individual household, was associated with GPS coordinates (Figure 1).

| Spatial collections
Intense suppression experiments based on pyrethroid spraying were conducted in 2013 and 2014 (Gunning et al., 2018) to test the predictions of a detailed Ae. aegypti population dynamics model (Magori et al., 2009). In brief, two areas of the city were identified as having relatively high densities of Ae. aegypti and were configured spatially in a way that allowed for a central spray sector with an outer buffer sector to act as an experimental control region ( Figure 2). To limit the impact of migration on resistance allele frequency, site dimensions were selected to be 3-5 times larger than the expected Ae. aegypti lifetime flight distance of approximately 150 m (Harrington et al., 2005).  were collected and stored as described above.

| DNA extractions and quantification
Male mosquitoes were chosen for genetic analysis throughout this study because female mosquitoes were typically the focus of virology and epidemiological studies, and therefore, more males were available in the repository. Using females would also have brought the risk of genomic contamination from male mosquitoes (via insemination) and from humans (via human blood feeding

| Genotyping
Allele-specific quantitative PCR and melting curve analysis (AS-PCR) was used to genotype all mosquitoes in duplicate for each of the mutations most commonly found in Central and South America (V1016I and F1534C). If the two reactions were not scored identically, the sample was discarded from further analysis. Mismatches were rare and typically due to nonamplification of a sample or because certain criteria for scoring were not met; that is, melting peak did not cross threshold. A smaller number (n = 92) of individuals were additionally genotyped at the V410L locus to verify the strong linkage disequilibrium that has been previously reported between it and locus V1016I (Saavedra-Rodriguez et al., 2018). Each mosquito genotyped at the V410L locus was also genotyped twice to ensure accuracy.

| Genotyping of V1016I
AS-PCR for the V1016I locus was based on the method reported by Real-Time PCR machine in the GSL, with the following thermal conditions: step 1-95°C for 3 min, step 2-95°C for 10 s, step 3-60°C for 10 s, step 4-72°C for 10 s, step 5-go to step 2, 39 times, step 6-95°C for 10 s, and step 7-melting curve 65-95°C, increment 0.2°C per 10 s plus a plate read.

| Genotyping of F1534C
AS-PCR for the F1534C locus was performed following the method reported by Yanola et al. (2011) with the following modifications.
The PCR volume was reduced to a total of 10 µl per reaction and con-

| Genotyping of V410L
The AS-PCR for the V410L locus was based on a protocol developed

| Statistical analyses
For each year, we conducted a logistic regression to compute confidence intervals on monthly resistance allele frequencies in the spray and buffer zones. We then used these regression models to compute contrast ratios between zones and their corresponding p-values in R via the "emmeans" package (Lenth, 2020). For specific months estimated a wide window for N e (25 -3000), but found that most populations averaged between 400 and 600 for this parameter. We set min_s = −1, max_s = 0, min_h = 0, and max_h = 1. In the WFABC model, genotype fitness assumes a selection benefit to the susceptible allele. Because the susceptible allele was decreasing during the sampled period, we expect a negative selection coefficient. The additional modeling that we implemented, described below, used relative genotype fitness formulas that assumed a selection cost to the susceptible allele. Therefore, we report the absolute value of the WFABC selection parameters to present them in terms of cost to the susceptible allele and to more easily compare results between the WFABC and genotype-frequency models. In our results, when dominance is zero, the heterozygote behaves exactly like the resistant homozygote and the resistance allele is considered dominant, whereas, when dominance is one, the heterozygote behaves exactly like the susceptible homozygote and the resistance allele is considered recessive.

| Parameter estimation using a genotypefrequency model
We use a genotype-frequency model to evaluate the spread of the resistance allele at the 1534 locus of the VGSC. We assumed the population was well mixed with random mating and that there were discrete (nonoverlapping) generations of mosquitoes, each lasting one month. We also assumed the population was isolated, without immigration or emigration of mosquitoes. These assumptions were used to maintain the simplicity and interpretability of the model.
The full population genetics model (Appendix S1) and sampling distribution form a hidden Markov model (HMM

| Linkage disequilibrium
Linkage disequilibrium between alleles at loci coding for V1016I and V410L was calculated for each year from 2010 to 2017 (  Notably, the Ile1016/Cys1534 haplotype was not detected until 2010, but then increased in frequency until 2014, when the city stopped spraying pyrethroids and switched to malathion, which has a different mode of action than pyrethroids (Figure 3).

| Temporal kdr haplotype frequencies
We observed a temporary increase in the wild-type haplotype spectively), as well as between h and R 0 (0.536) ( Figure S2). For example, without large samples in the first several generations, it may not be possible to distinguish a low initial frequency and large fitness cost from a higher initial frequency but lower cost. We also used the genotype-frequency time-series samples from pMCMC to construct mean genotype estimates and 95% credible intervals ( Figure S3). In later years, these genotype estimates suggest that the susceptible allele is maintained in the population in heterozygotes, which is a result of the low estimates for h.
We noticed that our pMCMC model appeared to predict a greater abundance of heterozygotes than were observed in the empirical data during the period when the Val1016/Cys1534 haplotype was increasing ( Figure S3). However, the HWE analysis shows only three months (January 2003, March 2003, and December 2003 where our observed genotype frequencies were not in HWE (Table S2).

| Spatial kdr haplotype frequencies
During 2013 and 2014, the haplotype Ile1016/Cys1534 was increasing in frequency in the overall population. We had hypothesized that within each year, the frequency of Ile1016/Cys1534 would increase faster in the experimental zone receiving insecticide treatment than in the buffer zone even though there was likely to be movement between the two areas. The sample size from the spray zone for the 2013 suppression experiment was relatively small due to the efficacy of the spraying regime, with 252 mosquitoes analyzed from the buffer zone but only 120 analyzed from the spray zone. After spraying, the resistance allele frequency increased in the spray zone but

| D ISCUSS I ON
We examined the temporal and spatial patterns of kdr alleles in We demonstrated that both resistance alleles, Ile1016 and Cys1534, experienced strong selective pressure and that Cys1534 appears to be partially dominant over the wild-type allele. In addition, we detected significant spatial heterogeneity on the scale of a few city blocks, which supports findings reported elsewhere (Deming et al., 2016;Estep et al., 2018;Grossman et al., 2019). Based on the extensive duration and area of sampling, this is the most thorough study of kdr population genetics in Ae. aegypti to date and it confirms the remarkable capacity for this species to adapt rapidly to insecticidal pressures. Temporal results indicate that, prior to the municipal citywide spraying of pyrethroids that commenced in 2002, the majority (>99%) of haplotypes in the population were wild type and <1%

| Haplotype dynamics
of the haplotypes were Val1016/Cys1534 (Figure 3). Individuals carrying the Val1016/Cys1534 haplotype could either represent standing genetic variation in the population due to mutation/selection balance or to historical selection of the population exposed to DDT during the hemisphere-wide Ae. aegypti eradication campaign.

| Selection and dominance
Our data provide field-based evidence that pyrethroid spraying selected for multiple kdr haplotypes. This is consistent with previous studies that demonstrated functional pyrethroid resistance caused F I G U R E 7 Frequency of Ile1016/Cys1534 haplotype by treatment group across the year 2013 (left) and 2014 (right). A 6-week experimental cypermethrin application, highlighted with yellow blocks, occurred in the spray zone from April 29 to June 3, 2013, and from April 28 to June 2, 2014. A citywide cypermethrin spray in response to a dengue outbreak occurred in February 2014 and is also highlighted with a yellow block. p-values for a test of equality between treatment groups by month were computed from generalized linear models (shown for only for p < 0.1) by kdr mutations in the laboratory (Du et al., 2013;Hirata et al., 2014) and other field-based studies that reported high frequencies of kdr mutations in resistant populations (e.g. Brito et al., 2018;Deming et al., 2016;Li et al., 2015;Linss et al., 2014).  (Linss et al., 2014). While Iquitos is a fairly isolated city with no roads coming into it from other cities, planes and boats arrive daily, which may facilitate mosquito immigration (Guagliardo et al., 2015).
We found evidence of selection pressure in favor of kdr mutations.
The true selection pressure, however, was likely underestimated because several assumptions of both WFABC and genotype-frequency models were violated. The selective environment was not constant because the insecticide used, frequency of sprays, and the number of houses sprayed in any given pyrethroid application varied across time and space (Table S1). The model also assumed that the loci are independent, but V1016I and F1534C are tightly linked in the genome. Despite the possible underestimation of selection strength, the estimated mean values indicate strong selection pressure for both resistance alleles.
The Cys1534 resistance allele was estimated to be partially dominant to the Phe1534 susceptible allele by both the WFABC model and our genotype-frequency model using pMCMC. While some previous investigators found evidence that kdr alleles are recessive, they used laboratory, phenotypic exposure assays to determine inheritance (Brito et al., 2018;Saavedra-Rodriguez et al., 2007). Others suggested an additive effect with a certain combination of kdr alleles (Ishak et al., 2015;Plernsub et al., 2016). The WFABC method employed in this paper determined inheritance based on the number of generations that an allele existed at low frequencies and did not include phenotypic data (Foll et al., 2014). Phenotypic resistance can be caused by multiple mechanisms, and while the kdr mutations can be useful genetic markers for populations under selective pressure, resistance itself is more complex and involves additional genes and metabolic mechanisms (Smith et al., 2019). In addition, the dominance of an allele can vary depending on environment and genetic background (Bourguet et al., 2000). All of this taken together may explain the estimated differences in dominance between results from this study and those from previous investigations.
The incomplete dominance estimated here by both the WFABC model and our own population genetics model predicted increases in the rate at which the allele frequencies would rise in the population from small frequencies. This is because most resistance alleles were present as heterozygotes, and higher dominance would give greater selective advantage to the heterozygote compared to the homozygous susceptible under spray conditions. Similarly, when the resistance allele was detected at high frequencies, the susceptible allele was present mostly as heterozygotes and quickly fell out of the population (Conner & Hartl, 2004).
While neither the WFABC model nor our own population genetics model accounted for spatial heterogeneity or migration, we know that these phenomena play important roles in Ae. aegypti population dynamics. Due in part to their limited flight dispersal, this species is thought to have strong spatial structuring. We also know that in the early years of spraying pyrethroids in Iquitos, that the population sizes were significantly reduced (Harrington et al., 2005;Hlaing et al., 2010;Morrison et al., 2008). It is possible that these factors are contributing to the rate of resistance evolution rather than the dominance of the Cys1534 resistance allele. Because the dominance estimate for Ile1016 indicated partial recessivity, future investigations may be able to clarify the roles that spatial heterogeneity, migration, and dominance play in the rate of evolution at these two loci.

| Spatial heterogeneity
As has been found in previous studies, we were able to detect heterogeneity in kdr frequencies at the citywide scale in Iquitos. This is best illustrated by the haplotype frequency data from 2007 and 2008. There appeared to be a decrease in the frequency of resistance alleles in those years (Figure 3), but this result was likely an artifact of an expanded collection regime and not due to an actual decrease in allele frequency across the city (Figure 4) and July 2004 ( Figure S3). For simplicity, our model assumes a homogeneous, randomly mating population. However, we know that Ae. aegypti tend not to fly far and migrate on average only 150 m in their lifetime (Harrington et al., 2005). In addition, the mosquitoes included in this analysis were collected over the Northern part of Iquitos, but collections from any given month were typically spatially limited to just a few blocks. Therefore, while the genotype frequencies for this analysis are in Hardy-Weinberg equilibrium for 19 out of 22 months (Table S2), our sampling design was likely to detect existing genetic structure. This is highlighted by the difference between the predicted and observed heterozygote frequencies.
Previous work indicated that the frequency of resistance alleles within an area can be highly heterogeneous; however, prior studies were unable to determine whether the heterogeneity was due to heterogeneity in insecticidal pressures or strong population structure and variation in initial resistance allele frequencies (Deming et al., 2016;Linss et al., 2014). We examined the outcomes from two suppression experiments over 16 and 44 calendar weeks in years 2013 and 2014, respectively. Significant differences were detected between the spray and buffer areas in both years. While initial resistance allele frequencies were higher in 2014 than in 2013, a large divergence in resistance allele frequency was detected following the experimental 6-week pyrethroid spray treatment, providing further evidence that genetic heterogeneity in Ae. aegypti can exist at a very fine spatial scale. In the 2013 experiment, trends indicated that resistance allele frequencies increased in the spray zone compared to the buffer zone. The larger confidence intervals were likely due to the small sample size available for analysis during that year. Fewer houses were sampled and the spraying regime was much more effective at reducing population size in 2013, likely due to lower levels of resistance compared to the 2014 sampling period.

| CON CLUS ION
This study demonstrated the evolution of two important kdr haplotypes in Ae. aegypti over the course of 18 years in Iquitos, Peru. We document a strong association between citywide pyrethroid spraying and dramatic increases in the frequency of specific kdr haplotypes, including two periods when dramatic increases in resistance allele frequency were detected. In our examination of the literature on mosquitoes and other insect pests, we could find no cases where a pest evolved so quickly to so few exposures to low or nonresidual insecticide applications. A combination of high estimated selection coefficients for both kdr haplotypes and the partial dominance that the resistance allele, Cys1534, exerted over the susceptible allele when under the selection pressure of pyrethroid insecticides likely contributed to the rapid rate of evolution. Results from spatial suppression experiments showed that kdr allele frequencies varied not only at the citywide scale, but also on a fine scale within the city.
Taken together, our results highlight the importance of widespread monitoring of Ae. aegypti for insecticide resistance from the beginning of first applications. Similarly, it will be especially important to adjust the spatial scale of monitoring efforts to better capture the potential heterogeneity in allele frequencies within an area. The future use of pyrethroids for vector control in Iquitos remains uncertain. Even if the kdr resistance haplotype frequencies fall in the absence of pyrethroid selection pressure, resistance is likely to rebound from low frequencies if previously used pyrethroids become used widely once again. Currently, the Ministry of Health is using the organophosphate, malathion, for vector control; however, the results of this study should caution against relying on the use of a single class of insecticides.
In fact, a wide range of traditional and novel techniques are being considered for vector control in the city. Possible strategies include placing more emphasis on larval habit management, using spatial repellents instead of ULV spraying, and targeted residual spraying with novel nonpyrethroids. Genetic-based control strategies such as Wolbachia for population suppression or replacement and the sterile insect technique are also being considered. It is likely that a combination of methods may need to be utilized to achieve optimal vector control. No matter which solutions that are implemented, resistance management should be a key goal from the beginning.