Insecticide resistance in Aedes aegypti from Tapachula, Mexico: Spatial variation and response to historical insecticide use

Background Insecticide use continues as the main strategy to control Aedes aegypti, the vector of dengue, Zika, chikungunya, and yellow fever. In the city of Tapachula, Mexico, mosquito control programs switched from pyrethroids to organophosphates for outdoor spatial spraying in 2013. Additionally, the spraying scheme switched from total coverage to focused control, prioritizing areas with higher entomological-virological risk. Five years after this strategy had been implemented, we evaluated the status and variability of insecticide resistance among Ae. aegypti collected at 26 sites in Tapachula. Methodology/Principal findings We determined the lethal concentrations at 50% of the tested populations (LC50) using a bottle bioassay, and then, we calculated the resistance ratio (RR) relative to the susceptible New Orleans strain. Permethrin and deltamethrin (pyrethroids), chlorpyrifos and malathion (organophosphates), and bendiocarb (carbamate) were tested. The frequencies of the substitutions V1016I and F1534C, which are in the voltage-gated sodium channel and confer knockdown-resistance (kdr) to pyrethroid insecticides, were calculated. Despite 5 years having passed since the removal of pyrethroids from the control programs, Ae. aegypti remained highly resistant to permethrin and deltamethrin (RR > 10-fold). In addition, following 5 years of chlorpyrifos use, mosquitoes at 15 of 26 sites showed moderate resistance to chlorpyrifos (5- to 10-fold), and the mosquitoes from one site were highly resistant. All sites had low resistance to malathion (< 5-fold). Resistance to bendiocarb was low at 19 sites, moderate at five, and high at two. Frequencies of the V1016I ranged from 0.16–0.71, while C1534 approached fixation at 23 sites (0.8–1). Resistance profiles and kdr allele frequencies varied across Tapachula. The variability was not associated with a spatial pattern at the scale of the sampling. Conclusion/Significance Mosquito populations respond to selection pressure at a focal scale in the field. Spatial variation across sites highlights the importance of testing multiple sites within geographical regions.


Introduction
Aedes aegypti is the main vector of several arboviruses, including dengue, Zika, chikungunya, and yellow fever. The control of this mosquito species is challenging, mainly because it is highly adapted to urban and suburban areas and because it is widely dispersed in endemic regions [1]. Except for yellow fever, safe and effective treatments or vaccines for these diseases are still under study. Therefore, the suppression of Ae. aegypti remains the cornerstone to prevent transmission and control of outbreaks of these diseases [2].
Effective vector control involves several strategies, such as the elimination of potential breeding sites, application of chemical insecticides, and implementation of biological control. However, the application of chemical insecticides has become a common form of control because as a control is highly efficient and can be implemented promptly [3]. The most used insecticides in vector control have been the organophosphates temephos, used as a larvicide, and malathion, used as an adulticide by ultra-low volume application (ULV). Pyrethroids were introduced as adulticides in most Latin American countries in the 1990s [3]. In Mexico, according to the Mexican official policy for vector surveillance and control [4], the adulticide ULV formulation of permethrin, bioallethrin, and piperonyl butoxide (PBO) was used for more than 10 consecutive years (1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010). In the following 3-4 years, the pyrethroid d-phenothrin + PBO was introduced. Subsequently, the use of organophospates returned in 2013, with chlorpyrifos and malathion being used as adulticides, while carbamates were recommended for indoor residual application [5].
The prolonged use of pyrethroid insecticides resulted in the evolution of resistance to them in Ae. aegypti worldwide, including Mexico, where failures in dengue control strategies are due in part to resistance [5]. Given that resistance to insecticides has been reported in populations of Ae. aegypti globally, the World Health Organization (WHO) recommends testing to ensure an effective insecticide management program. Decisions based on evidence of the resistance and/or susceptibility of Ae. aegypti will ensure a better selection of insecticides in vector control programs [6].
Two mechanisms of resistance to insecticides have been identified: resistance due to the enhanced metabolism of insecticides and insensitivity at the target site of the insecticides. Both mechanisms are involved in resistance to pyrethroids [7]. Knockdown resistance (kdr) refers to a phenomenon in which insects are not knocked down immediately after exposure to pyrethroids. kdr is caused by specific mutations at the voltage-gated sodium channel (VGSC), which is the target site for pyrethroids and DDT [8]. The amino acid substitutions V1016I [9], F1534C [10], and V410L [11,12] frequently have been associated with resistance to pyrethroids. Once these mutations are fixed in a population, reversion to susceptibility is difficult to achieve [9]. Therefore, the detection and characterization of kdr mutations in mosquito populations before resistance fixation occurs is essential for insecticide management strategies.
In Mexico, Chiapas is one of the states with the highest rate of endemic dengue cases. In particular, the city of Tapachula reports the highest incidence of dengue in the state [13], which is attributed to the proliferation of vectors that transmit emerging and re-emerging diseases. Under the region's tropical climate conditions, Ae. aegypti maintains high densities throughout the year. Consequently, dengue and other arboviruses transmitted by this vector have been prevalent in the region for a long time [14].
In the context of insecticide resistance management, we investigated the status of insecticide resistance to five insecticides, including two pyrethroids (permethrin and deltamethrin), two organophosphates (chlorpyrifos and malathion), and one carbamate (bendiocarb) and the spatial distribution of such resistance in populations of Ae. aegypti throughout the city of Tapachula. We expect that after 5 years of heavy use of organophosphates and the removal of pyrethroids from vector control campaigns, pyrethroid resistance will be lower, whereas organophosphate resistance will appear in focal points of the city. We tested 26 collection sites located in the city of Tapachula. Each collection site consisted of nine blocks and these were selected based on vehicle access for outdoor spraying. To minimize the effects of mosquito migration by flight range (50-150 m), sites were located at least 250 m apart. The spatial correlation between resistance and geographical distance was calculated for the 26 collection sites. In addition, since Tapachula's vector control program uses a quadrant subdivision for spraying activities, we included a second analysis to test this source of variation by assigning sites to one of the four cardinal geographical quadrants (NE, NW, SE, and SW).

Collections
The study was conducted in the city of Tapachula, Chiapas, located in southern Mexico at 177 meters above sea level. The 26 collection sites located in four quadrants in the city: Northwest (NW), Northeast (NE), Southwest (SW), and Southeast (SE) are shown in Table 1. The biological material was collected from January to April 2018 using ovitraps of 1-L capacity [15]. Twelve ovitraps were installed at each collection site. Ovitraps were made by hand with transparent, inert, non-toxic polypropylene (PP) cups of a 1 L capacity, and painted black on the outside following the guidelines for Entomological Surveillance with Ovitraps of the Mexican Ministry of Health [15]. The interior of each ovitrap was lined with Whatman filter paper (No. 1) and filled with water to ¾ capacity; the paper was replaced weekly up to five times. The egg strips were transported to the insectary of the Regional Center for Research in Public Health/National Institute of Public Health (CRISP/INSP). The egg strips were submerged in 4 L of water in plastic containers (40 cm x 30 cm x 15 cm). On the third and sixth day, the hatched larvae were fed a diet of Harlan 5001 proteins, with 0.4 gr or 0.8 gr / 1.2 L for 1 st -2 nd stadium and 3 rd -4 th stadium at the 3 rd and 6 th day, respectively.
Aedes aegypti mosquitoes were identified to species and placed in cages (30 cm 3 ). Females were bloodfed from rabbit (under accepted procedures by the Ethical Commission of the Instituto Nacional de Salud Pública) to obtain the F 1 generation. Environmental conditions consisted of a temperature of 27 ± 2˚C, 70-80% humidity, and a 12:12 hour photoperiod. We used the insecticide-susceptible New Orleans reference strain of Ae. aegypti, provided by Dr. William Black and maintained in the CRISP/INSP insectary.

Bioassays
The F 1 adults were exposed to the insecticides using a modified CDC bottle bioassay (Centers for Disease Control) [16]. Sigma brand technical grade insecticides were used to determine the lethal concentrations that killed 50% (LC 50 ) at each site. The pyrethroids permethrin (Type I) and deltamethrin (Type II), the organophosphates malathion and chlorpyrifos, and the carbamate bendiocarb were used to represent the toxicological groups used by vector control programs in Mexico.
To determine the LC 50 , we tested five to six insecticide concentrations, which caused 10 to 90% mortality, in four replicates. Each insecticide LC 50 required approximately 500 mosquitoes. Table 2 shows the insecticide concentrations (μg/bottle) used to coat 250 ml Wheaton bottles using acetone as the solvent. During the bioassay, 15 to 20 (2-3 day old) females were gently aspirated into each bottle. The knockdown effect was recorded every 10 minutes for 1 hour. After 1 hour of exposure, the mosquitoes were transferred to plastic containers and maintained in the insectary to observe the mortality at 24 hours. The LC 50 of each insecticide was also determined for the susceptible New Orleans reference strain (NO) using a different set of insecticide concentrations ( Table 2). Each insecticide LC 50 was replicated at least five times during a 7-month period. As control, a bottle impregnated only with acetone was used each time a test with field or susceptible mosquitoes was run.
The LC 50 , 95% confidence intervals, slope, intercept, and p values were determined using the binary logistic regression model with QCal software [17]. The null hypothesis (Ho) assumes the observed mortality curve adjusts to a binary logistic regression model. Thus, we expected p values higher than 0.05 to accept the Ho. When the Ho was rejected, the bioassay was repeated.
To estimate the level of resistance among sites, a resistance ratio (RR) was calculated by dividing the LC 50 of the field sites by the LC 50 of the NO strain. The RR criterion according to Mazzarri and Georghiou [18] classifies high resistance as RR values greater than 10, moderate resistance as RR values between 5 and 10, and low resistance as RR values less than 5.

Genotyping kdr-associated mutations
Genomic DNA was isolated from 50 F 1 individual female mosquitoes from each collection site following the method of Black and DuTeau [19]. The DNA was resuspended in TE buffer (10 mM Tris-HCl, 1 mM EDTA pH 8) and stored at -20˚C. The V1016I and F1534C mutations were genotyped according to the protocols of Saavedra-Rodríguez et al. [9] and Yanola et al. [10], respectively. The genotype and allelic frequencies were tested for Hardy-Weinberg (HW) equilibrium. The null hypothesis is that equilibrium is present in the population, which was verified with a chi-square test (df = 1 and p value > 0.05). We tested the spatial variation of the LC 50 s between the quadrants in the city using a linear model and ANOVA in R (car package). To test the hypothesis of resistance correlation with space, we created Moran's I correlograms as implemented in PASSaGE 2.0 [20]. Mosquitoes from different collection sites were considered neighbors if the sites were within 250 meters of each other. We expected that the LC 50 s or kdr frequencies would be associated with geographical distance (i.e., that closer neighbor sites would show similar resistance levels, compared to those farther away). A second analysis to test the variation of the LC 50 s and kdr frequencies between and within quadrants using a linear regression model and ANOVA in R (car package) was conducted. Since the city is uniformly sprayed during a cycle, we did not expect variation between or within quadrants. Correlation between kdr frequencies and LC 50 s for permethrin and deltamethrin was tested using a Spearman test.

Results
The geographic distributions of the resistance ratios (RR) for each insecticide in the 26 sites in Tapachula are shown in Fig 1. The LC 50 and confidence intervals for each of the five insecticides are shown in S1 Table. For the pyrethroids, we observed high levels of resistance widespread across sites. Fig 2A shows the permethrin RRs across Tapachula. High RRs were identified at 24 sites (RR from 11.4 to 43.1-fold). Only two sites-NE-3 and NW-2-showed moderate RRs (5.3 and 5.9-fold, respectively). The variation in RRs among quadrants was not significant (F = 0.56, df = 3, p value = 0.64). For deltamethrin, high RRs were determined in all 26 sites (10.6 to 101-fold). The variation among quadrants was not significant (F = 1.08, df = 3, p value = 0.37). Except for SW, all quadrants had at least one site with RR higher than 90-fold ( Fig 2B).

Kdr-associated mutations
Genotype frequencies at the V1016I and F1534C loci in the voltage-gated sodium channel gene were determined in a sample of 45-50 individuals from each site ( Table 3). The allele frequencies of the resistant allele I1016 fluctuated from 0.16-0.71. The lowest allele frequency (0.16) was scored for NE-3, whereas the highest frequency was from NW-6 (0.71). The remaining sites ranged from 0.2 to 0.5. Except for NE-2 and SW-4, the genotype frequencies at the V1016I loci were in HW equilibrium.
High allele frequencies of the resistant C1534 allele were determined at 22 of the 26 sites, ranging from 0.85 to 1.0. Lower frequencies (0.38-0.41) were found in NE-3, NW-5, and NW-7. While NE-7 was calculated with an intermediate value of 0.61. Most sites were in HW disequilibrium due to fixation of the resistant allele. We conducted a Spearman correlation test between the pyrethroid LC 50 s and the expected frequencies of resistant homozygous genotypes. We found significant correlation coefficients among permethrin LC 50 s, I1016/I1016 homozygotes (S = 2588, rho = 0.53, p value = 0.002), and C1534/C1534 homozygotes (S = 1966, rho = 0.515, p value = 0.004). Although it is known that C1534 shows protection only against permethrin (12), for deltamethrin a significant correlation was observed between the LC 50 and I1016/I1016 homozygotes (S = 2643, rho = 0.467, p value = 0.008) and the C1534/C1534 homozygotes (S = 2945, rho = 0.507, p value = 0.002). However, the significance for both insecticides disappeared when observations from the New Orleans reference strain were removed.
To assess the correlation of LC 50 s with space, we generated Moran's I correlograms for each of the five insecticides (Fig 4). The analysis included all 26 collection sites. We did not detect a discernable pattern in any of the tested insecticides. We expected a positive correlation (Moran's I statistically > 0) between nearby sites, then as the distance increased (between the samples) the correlation would decrease, and later would turn negative (Moran's I statistically < 0). However, this was not observed. Although, few of the distance classes were statistically different from zero (eg. bendiocarb at 3250 m, 3750 m, 4750 m, and 6250 m; malathion at 1500 m; and deltamethrin at 3750 m, and 4250 m), a caveat in our analysis is the possibility that there is autocorrelation at smaller distances than the ones we selected (x < 250 m). Our experimental design was not geared towards the detection of spatial correlation at smaller distances; there were a small number of samples below 250 meters.

Discussion
Efforts to control Ae. aegypti populations are hindered by widespread insecticide resistance worldwide. Local insecticide resistance monitoring is necessary for the design of specific and successful resistance management programs [21]. In Latin America, pyrethroids have been used for adult mosquito control since the 1990s. The switch to pyrethroids was based on environmental concerns that led to the use of less toxic classes of insecticides [22]. In Mexico, vector control programs implemented the use of permethrin in 1999 and continued their use until 2010. Local selection pressure caused a rapid evolution of pyrethroid resistance in Ae. aegypti across Mexico [9,[23][24][25][26][27], resulting in policy modifications that recommended the use of insecticides with different toxicological modes of action.
In Tapachula, vector control programs replaced the use of permethrin with a different Type I pyrethroid (d-phenothrin + piperonyl butoxide) from 2010 to 2013. In 2013, pyrethroids were replaced by the organophosphate chlorpyrifos, and in 2017, by malathion. This study reveals the current status and response of local Ae. aegypti populations to these insecticide shifts. Despite the switch to organophosphates in the last 5 years, we observed that high levels of pyrethroid resistance remain widespread in Tapachula. An assumption in insecticide resistance management is that insecticide resistance has negative fitness costs. Therefore, when insecticide pressure is removed, populations are expected to reverse to susceptibility [28,29]. Currently, we are conducting a study to determine the degree of loss of resistance to pyrethroids from 2016 to 2020 in this study area, which will demonstrate whether mosquito populations in Tapachula are undergoing a process of decreasing resistance that will take several years. Another explanation is that pyrethroid resistance is maintained in Ae. aegypti populations by the domestic use of pyrethroids [30]. Surveys in Merida, Mexico, indicate that 85% of households took action to kill pests, and 89% exclusively targeted mosquitoes. Interestingly most of the aerosol spray cans contained pyrethroid insecticides [31].
Interestingly, RRs for deltamethrin-a Type II pyrethroid-were higher than permethrin RRs across sites. Deltamethrin was authorized by CENAPRECE for indoor residual use in 2009 for control of the malaria vector, but its use was restricted to rural areas. Therefore, direct selection pressure from the use of deltamethrin in public health is unlikely to be responsible for the high RRs in Ae. aegypti from Tapachula. Although all pyrethroids act at the same target site, the variability of resistance to their different types is attributed to different binding sites for Type I and Type II pyrethroids at the voltage-gated sodium channel. Additionally, the presence of enzymes that have a greater affinity to metabolize specific molecules within the same toxicological group might explain this variability [32].

Fig 2. Pyrethroids resistance ratios (RRs) of Aedes aegypti collected in 26 sites across Tapachula in 2018. A) Permethrin and B)
Deltamethrin. Dots represent the RR 50 with 95% confidence intervals for each site. Horizontal lines indicate the threshold for low resistance (< 5-fold), moderate resistance (5-to 10-fold) and high resistance (> 10-fold). https://doi.org/10.1371/journal.pntd.0009746.g002 Knockdown resistance (kdr) is a major mechanism of pyrethroid resistance in Ae. aegypti from Mexico. In this study, we measured the frequency of this mechanism by molecular tests that identify mutations that confer changes to amino acids in the VGSC. The allele frequencies of the resistant allele I1016 ranged from 0.4 to 0.7, and for the resistant allele C1534, from 0.85 to 1.0 (except for three sites that had~0.4). Historical data of kdr mutations indicated that C1534 confers low level of resistance on its own, and that resistance increased dramatically when I1016 evolved from the V1,016/C1,534 haplotype in field mosquito collected in different

PLOS NEGLECTED TROPICAL DISEASES
places from Mexico [33]. Those results demonstrated that I1016 was unlikely to have evolved independently, and that both mutations need coexist in the same mosquito in order to confer higher levels of resistance. Moderate correlations were significant between the resistant allele frequencies and the RRs for permethrin and deltamethrin only when including the New Orleans datapoints. This significance might be explained by most of the allele frequencies being distributed within a small range of variability. This study was conducted after chlorpyrifos had been used for 5 years in outdoor spraying by vector control programs. Our results provide evidence of the response of Ae. aegypti populations to chlorpyrifos pressure. Ten sites showed low RRs, 15 sites showed moderate resistance, and one site was highly resistant. Interestingly, Ae. aegypti from all 26 sites were susceptible or had low RRs to malathion, thereby indicating that resistance to chlorpyrifos does not predict the lack of effectivity of malathion. Additionally, the RRs to bendiocarb were variable: mosquitoes from 19 sites had low RRs, those from three were moderate, and those from two were highly resistant. Only a few sites showed moderate to high resistance to both chlorpyrifos and bendiocarb (NE-5, NW-6, and SE-4). The lack of cross-resistance between organophosphates and carbamates suggests that the resistance mechanisms are not due to the insensitivity of their target site (the acetylcholinesterase) [34] and, in fact, no mutations have been found in ace-1 gene in Aedes aegypti [35].
A survey in Veracruz, Mexico, identified high RRs to chlorpyrifos in Cosoleacaque (RR = 13.9), moderate RRs in Poza Rica (RR = 7.9), and low RRs in five sites in Veracruz [36]. By using a discriminating dose of 50 μg/bottle and 85 μg/bottle for 30 minutes, two additional studies were able to identify chlorpyrifos resistance in Mexico [26,37]. Since neither of these studies found a history of chlorpyrifos use in vector control programs, the resistance might be explained by indirect exposure to chlorpyrifos through the extensive use of this insecticide to control agricultural pests [36].
During vector control programs, the city of Tapachula is uniformly sprayed, using the same insecticide, frequency and intensity. More yet, we selected sites based in their accessibility to spraying-vehicles. Assuming that no spatial heterogeneity in frequency and intensity of spraying, we did not expect the high levels of variation in resistance profiles across the city. For example, significant heterogeneity in the frequency of kdr haplotypes was detected in Ae. aegypti collected between city blocks in a town of Yucatan, suggesting that selection for these haplotypes occurs at a fine spatial scale (37). However, in contrast to our study, insecticide application was highly variable in space and time, creating a mosaic of selection pressures. In our study, some sources of heterogeneity could occur from mosquito migration from untreated sites due to vehicle inaccessibility, including parks, cemeteries, steep and unpaved streets. A second source of insecticide pressure is by use of household aerosol insecticides. For example, in a previous study from Merida, Mexico approximately 87% of households used commercially available pyrethroid products to control mosquitoes in their homes (31). Future studies should include an assessment of this source of selection pressure in Tapachula.
The spatial variability in insecticide resistance observed across the 26 sites in Tapachula is likely associated with the presence or appearance of "hot spots or dengue foci," which contribute to the persistent transmission of the diseases and therefore to focal areas with greater spray intervention [38]. In addition, the spatial variability of resistance highlights the importance of evaluating resistance in multiple sites within a defined geographic area for the application of appropriate vector control decisions. Although no geographical correlation/association/pattern between resistance was found in Tapachula, more specific and finer environmental characteristics must not be discarded in future studies. A previous study used mitochondrial ND4 haplotypes to determine gene flow patterns among 38 Ae. aegypti coastal collections in Mexico [39]. Three genetic clusters were identified, the Northeast, Pacific, and Yucatan peninsula. For all sites, genetic distances remained small below geographic distances of 90 km and became large at distances >150 km. The Pacific cluster had the highest gene flow and diversity. A second study in the Yucatan Peninsula showed high gene flow occurring across 27 Ae. aegypti collection sites located up to 150 km of distance. Single nucleotide polymorphism (SNPs) at eleven loci did not vary across sites, suggesting high levels of gene flow. In contrast, insecticide resistance loci, including kdr alleles (I1016 and C1534) were highly variable across sites, indicating that insecticide resistance offsets the homogenizing effects of gene flow [40]. In this study, we assume complete gene flow among collection sites because: 1) Tapachula belongs to the Pacific cluster, 2) Ae. aegypti is well established throughout the year and, 3) collection sites are within 10 km of distance. However, this remains to be tested.

Conclusion
Despite more than 5 years having passed since the removal of pyrethroids from vector control programs in Tapachula, high levels of pyrethroid resistance and kdr-associated alleles persist in Ae. aegypti populations. Future resistance surveys will reveal if pyrethroid resistance is maintained in mosquito populations. We observed that, after 5 years of chlorpyrifos use in vector control programs, more than 50% of the sites have moderate to high chlorpyrifos resistance but complete susceptibility to malathion. Since malathion was introduced later in 2017, future studies to evaluate the selection of malathion resistance in the field are needed. Two different analyses were conducted 1) the spatial analysis included all 26 sites and, 2) the quadrant analysis to identify operational sources of heterogeneity. The quadrant analysis doesn't include a geographical component and has limitations. Insecticide resistance varied spatially, most likely as a consequence of the pattern of insecticide use combined with environmental factors. Based on the results of our study, we suggest that both of the studied organophosphates and the carbamate remain viable options for use in the control strategy for this vector. The return to a pyrethroid (at least permethrin and deltamethrin) for outdoor spraying is recommended when the levels of resistance have decreased to RR less than 10-fold and once mechanisms other than kdr have been elucidated for pyrethroid resistance.
Supporting information S1 Table. Lethal concentration to kill 50% (LC 50 ) for five insecticides at 26 Aedes aegypti sites in Tapachula, Chiapas, Mexico. The LC 50 is in micrograms (μg) of active ingredient per bottle. The 95% confidence intervals around the LC 50 are enclosed in parentheses. p values higher than 0.05 indicate the observed data fit the binary logistic regression model.