Co-occurrence of kdr Mutations V1016I and F1534C and Its Association With Phenotypic Resistance to Pyrethroids in Aedes aegypti (Diptera: Culicidae) Populations From Costa Rica.

Aedes aegypti (Linnaeus, 1762) is considered the most important mosquito vector species for several arboviruses (e.g., dengue, chikungunya, Zika) in Costa Rica. The primary strategy for the control and prevention of Aedes-borne diseases relies on insecticide-based vector control. However, the emergence of insecticide resistance in the mosquito populations presents a significant threat to these prevention actions. The characterization of the mechanisms driving the insecticide resistance in Ae. aegypti is vital for decision making in vector control programs. Therefore, we analyzed the voltage-gated sodium channel (VGSC) gene for the presence of the V1016I and F1534C kdr mutations in Ae. aegypti populations from Puntarenas and Limon provinces, Costa Rica. The CDC bottle bioassays showed that both Costa Rican Ae. aegypti populations were resistant to permethrin and deltamethrin. In the case of kdr genotyping, results revealed the co-occurrence of V1016I and F1534C mutations in permethrin and deltamethrin-resistant populations, as well as the fixation of the 1534C allele. A strong association between these mutations and permethrin and deltamethrin resistance was found in Puntarenas. Limon did not show this association; however, our results indicate that the Limon population analyzed is not under the same selective pressure as Puntarenas for the VGSC gene. Therefore, our findings make an urgent call to expand the knowledge about the insecticide resistance status and mechanisms in the Costa Rican populations of Ae. aegypti, which must be a priority to develop an effective resistance management plan.

environmental management and treatment of breeding sites with larvicides. For the adult populations, several pyrethroids have been used in the past, such as deltamethrin, cyfluthrin, and cypermethrin (Marín Rodriguez et al. 2009, Bisset et al. 2013. Currently, permethrin + esbioallethrin + piperonyl butoxide and lambdacyhalothrin are being used (Ministerio de Salud de Costa Rica 2019). Focal studies have reported a decrease in the susceptibility against pyrethroids in Ae. aegypti from Puntarenas and Limon (Bisset et al. 2013;Troyo 2014, 2016;Calderón-Arguedas et al. 2018). Moreover, biochemical and synergists assays suggest that carboxylesterases and cytochrome P450s could be involved in pyrethroid resistance in these same areas. One of the main mechanisms of pyrethroid resistance in insects is the target-site insensitivity (Moyes et al. 2017), which has been poorly described in Costa Rica (Chaves et al. 2015). Nonsynonymous mutations on the knockdown resistance (kdr) region of the voltage-gated sodium channel (VGSC) gene are responsible for this pyrethroid insensitivity (Du et al. 2016). A study published in 2003 reported seven novel mutations in Ae. aegypti strains from Asia, Africa, the Caribbean, South America, and Oceania (Brengues et al. 2003), but only five have been associated with pyrethroid resistance: V410L, S989P, I1011M, V1016G/I, and F1534C (Du et al. 2013, Hirata et al. 2014, Haddi et al. 2017, Chen et al. 2019.
Molecular screening of the Latin American populations of Ae. aegypti evidenced the emergence of the kdr mutation V1016I in the region (Saavedra- Rodriguez et al. 2007) and its co-occurrence with the I1011M mutation, which was recently found in Panama (Murcia et al. 2019). The first report of the F1534C mutation was published in 2011 in the Ae. aegypti population from Thailand (Yanola et al. 2011). Since then, the F1534C mutation has been documented worldwide (Linss et al. 2014, Kushwah et al. 2015, Hamid et al. 2018, Maestre-Serrano et al. 2019, Yaméogo et al. 2019, except Central America. To the best of our knowledge, kdr mutations have not been studied in pyrethroid-resistant Ae. aegypti from Costa Rica.
To fill the gap about the mechanisms associated with insecticide resistance, we analyzed the VGSC gene for the presence of the V1016I and F1534C kdr mutations in pyrethroid-resistant Ae. aegypti populations and its association with insecticide resistance from two of the most affected provinces in Costa Rica, Limon, and Puntarenas.

Study Area
The two collection sites selected for this study were Barrio El Carmen (9.978, −84.829), located in Canton Puntarenas, a municipality of Puntarenas Province, and Barrio Colina (9.986,041) in Canton Limon, a municipality of Limon Province (Fig. 1). Municipalities of Puntarenas and Limon were selected based on its historical persistence of Aedes infestation rates and high risk for arbovirus infections (Troyo et al. 2008, Calderón-Arguedas et al. 2009, Mena et al. 2011, Rodríguez et al. 2014, Marín and Díaz 2015, Rodríguez and Diaz 2015. Additionally, Puntarenas and Limon are the home of the two most important maritime ports in the Pacific and Atlantic coast, respectively, and tourism is predominant in both municipalities. Vector control is mainly responsive and conducted in epidemic-prone zones, areas with high Aedes infestation rates, and where cases are reported. It includes indoor space spraying with the pyrethroid lamdacyalothrin applied in positive houses for arbovirus infections and a minimum of 100 m around the primary house and complemented occasionally with aerial ultra-low volume spraying with permethrin + esbioallethrin + piperonyl butoxide (Aqua-Reslin), and

Mosquito Sampling and Rearing
Aedes mosquito eggs were collected by the Costa Rican local vector surveillance and control units. In total, 92 ovitraps were placed in Barrio el Carmen (Puntarenas Province) and between 80 and 100 ovitraps in Barrio Colina (Limon Province). In Barrio el Carmen, ovitraps were placed in the intradomiciliary and peridomiciliary of houses and premises (hotels, restaurants, and small businesses), whereas ovitraps were placed inside houses in Barrio Colina. Ovitraps were placed for both sites where previously positive containers for Aedes spp. have been reported by the local vector control program. Ovitraps were visited weekly to collect mosquito eggs during consecutive weeks from April 27 to the end of May 2017. Eggs collected from the ovitraps for each site were pooled, hatched, and the larvae reared into adults under insectary conditions (28 ± 2°C; 70 ± 10% RH; and 12:12 (L:D) h photoperiod). An F 1 generation was obtained to get enough individuals to perform the bioassays. All emerged adults were identified as Ae. aegypti based on morphological characters following a pictorial key (Rueda 2004), and then they were provided access to cotton pads soaked in a 7% sugar solution until bioassays were realized.

Adult Bioassays
As vector control programs in Costa Rica rely mostly on pyrethroids, insecticide susceptibility of Costa Rican Ae. aegypti populations were assessed using the CDC bottle bioassay with a member of each type of pyrethroids, Type I (permethrin) and Type II (deltamethrin; Brogdon and McAllister 1998). For the bioassays, a diagnostic dose of 15 µg per bottle was used for permethrin and 10 µg per bottle for deltamethrin, and both doses were prepared with absolute ethanol and from technical grade insecticides (ChemService Inc., West Chester, PA). From each population, a total of 100 F1 females (2-5 d old, nonblood fed) were exposed to the diagnostic dose of deltamethrin or permethrin for 30 min as diagnostic time in four batches of 25 individuals per bottle. In parallel, 25 F1 females (2-5 d old, nonblood fed) were exposed in one bottle impregnated with ethanol as control per insecticide for 30 min. The New Orleans strain was used as insecticide-susceptible reference strain to validate the appropriate performance of the bioassays, which showed full susceptibility to permethrin (n = 101 females, 2-5 d old, nonblood fed) and deltamethrin (n = 100 females, 2-5 d old, nonblood fed). Based on the incapacitation rate at 30 min, each population was classified according to the standard World Health Organization (WHO) criteria as follows (WHO/GMP 2016): a population with an incapacitation rate > 98% was classified as susceptible (S), an incapacitation rate between 90 and 98% suggested possible resistance requiring further tests, and an incapacitation rate < 90% indicated confirmed resistance (R). Resistant and incapacitated (susceptible) mosquitoes were preserved individually at −20°C for molecular assays.

Kdr Genotyping
To calculate the general allelic frequencies of the F1534C and V1016I mutations of both populations from Costa Rica, 50 unexposed female Ae. aegypti (2 to 5 d old, nonblood fed) per site were genotyped. For this purpose, genomic DNA (gDNA) was extracted from mosquitoes using DNAzol (Thermo Fisher Scientific, Waltham, MA) according to manufacturer's instructions with modifications. Briefly, 100 µl of DNAzol was used to homogenize a whole mosquito, and DNA pellets were resuspended in 50 µl of DEPC water. Subsequently, the genotype for each kdr mutation was determined by a melt curve analysis of a real-time PCR using allelespecific primers for the V1016I mutation (Val1016f, Ile1016f, and Ile1016r) reported by Saavedra- Rodriguez et al. (2007) and for the F1534C mutation (F1534-f, C1534-f, and CP-r) reported by Yanola et al. (2011;Supp  PCR reaction for the V1016I mutation contained 10 µl of 2× PowerUp SYBR Green Master Mix (Thermo Fisher Scientific), 1.0 µl of primers Val1016f (25 µM) and Iso1016f (25 µM) and 0.5 µl of Iso1016r (6.25 µM), 1 µl of gDNA (5-15 ng) in a final volume of 20 µl. The PCR conditions were as followed: 95°C for 2 min, 40 cycles of 95°C for 10 s, 60°C for 10 s, and 72°C for 30 s, with a final stage of 10 s at 95°C in a QuantStudio 3 Real-Time PCR System (Applied Biosystems, Waltham, MA). The melt curve conditions were 65-95°C, with an increase of 0.2°C each 10 s. Based on the melting temperatures, samples were considered as mutant homozygous if they presented a unique peak of 79°C, wildtype homozygous if they psresented a unique peak of 85°C, and heterozygous if they presented both of the previous peaks (Supp Fig. 1 [online only]). In the case of F1534C mutation, each reaction contained 10 µl of 2× PowerUp SYBR Green Master Mix (Thermo Fisher Scientific), 0.066 µl of C1534-f primer (50 pmol/µl), 0.2 µl of primers F1534-f (50 pmol/µl) and CP-r (50 pmol/µl), 1 µl of gDNA (5-15 ng) in a final volume of 20 µl. The PCR conditions were as follows: 95°C for 3 min, 40 cycles of 95°C for 10 s, 57°C for 10 s, and 72°C for 30 s, with a final stage of 10 s at 95°C in a QuantStudio 3 Real-Time PCR System (Applied Biosystems). The melt curve conditions were 65-95°C, with an increase of 0.5°C every 5 s. Based on the melting temperatures, samples were considered as mutant homozygous if they presented a unique peak of 85°C, wildtype homozygous if they presented a unique peak of 80°C, and heterozygous if they presented both of the previous peaks (Supp Fig. 2 [online only]). Previously genotyped Ae. aegypti by sequencing were used as positive controls for the PCR assays.

Association between Phenotypic Resistance to Pyrethroids and kdr Mutations
To observe the association between the presence of V1016I and F1534C mutations to pyrethroid resistance, similar subsample sizes of resistant and incapacitated mosquitoes were genotyped for each location, as has been described earlier. For Limon, 21 permethrinresistant and -incapacitated mosquitoes, and 13 deltamethrinresistant and -incapacitated mosquitoes were genotyped. This was limited by the incapacitation rates of the bioassays and the availability of biological material. For Puntarenas, 32 permethrin-resistant and 38 permethrin-incapacitated mosquitoes, and 37 deltamethrinresistant and 39 deltamethrin-incapacitated mosquitoes were genotyped.

Sequencing of the VGSC Gene
To validate the genotypes obtained from the kdr genotyping, two different PCR assays were used to amplify a partial fragment of domain IIS6 (between 562 and 579 bp, Supp  Table S3 and S4 [online only]), respectively, which represents the different genotypes identified by the real-time PCR. The first assay amplified the domain IIS6 using the primers IIS5-6 F and IIS5-6 R (spanning the codons 989, 1011, and 1016, see Supp Table S2 [online only]) originally reported by Al Nazawi et al. (2017). Each PCR reaction contained 1× GoTaq Hot Start Colorless Master Mix (Promega, Madison, WI), 0.2 µM of each primer (IIS5-6F and IIS5-6R), and 5-10 ng of gDNA in a total volume of 25 µl. The cycling parameters included an initial step of denaturation at 95°C for 2 min, followed by 40 cycles of 94°C for 30 s, 50°C for 30 s, and 72°C for 1 min, with a final extension step of 10 min at 72°C. For the second assay, a region of the domain IIIS6 was amplified using the primers AaNa31F and AaNa31R (spanning the codons 1520 and 1534, see Supp Table S2 [online only]) reported by Harris et al. (2010). Each PCR reaction contained 1× Colorless GoTaq Flexi Buffer (Promega), 1.5 mM MgCl 2 (Promega), 0.2 mM dNTPs (Promega), 0.5 µM of each primer (AaNa31F and AaNa31R), 2.5 U of GoTaq Hot Start Polymerase (Promega), and 5 ng of gDNA in a total volume of 50 µl. The cycling parameters included an initial denaturation at 95°C for 2 min, followed by 35 cycles of 94°C for 30 s, 59°C for 30 s, and 72°C for 1 min, with a final extension step of 10 min at 72°C. Both PCR assays were performed using the SimpliAmp Thermal Cycler (Applied Biosystems).
All PCR products were visualized on a 2% agarose gel stained with ethidium bromide under UV light. Then, 20 µl of PCR product was purified using ExoSAP-IT PCR Product Cleanup Reagent (Thermo Fisher Scientific) according to the manufacturer's instructions. The PCR products purified were sequenced by Psomagen Inc. (Rockville, MD) with the primers used in the PCR assays at a concentration of 5 pmol/µl. DNA sequences were manually edited using FinchTv 1.4.0, and consensus sequences were assembled and aligned using MEGA 7 (Kumar et al. 2016). All codons were analyzed for nonsynonymous substitutions using the reference sequence [AAEL023266] of the Ae. aegypti VGSC gene deposited in VectorBase.

Statistical Analysis
Allelic frequencies and their respective confidence intervals were calculated for all subpopulations (unexposed, resistant, and incapacitated) following Agresti-Coull method (Agresti and Coull 1998).
Chi square for testing Hardy-Weinberg equilibrium (HWE), and the inbreeding coefficient (F is ) were calculated in populations where this was applicable. The association between genotype and phenotype was performed through an odds ratio with 95% confidence analysis using a recessive model, and the Fisher's exact test (P < 0.05 was considered significant).

Results and Discussion
The results of the CDC bottle bioassays presented in Table 1 confirmed the presence of resistance to deltamethrin and permethrin in both Costa Rican Ae. aegypti populations. Limon exhibited higher incapacitation rates for permethrin and deltamethrin, 79 and 87%, respectively, compared with Puntarenas with 35 and 62% of incapacitation, respectively. Incapacitation rates for permethrin were lower than deltamethrin in the two Ae. aegypti populations from Costa Rica. A previous study in 2013 found suspected resistance to deltamethrin (incapacitation rate of 90.5%) in a Puntarenas Ae. aegypti population (Bisset et al. 2013). Therefore, our results suggest that deltamethrin resistance in Puntarenas has increased in the last years. Although, this apparent increase in resistance needs to be analyzed further due to possible variations in test protocols, sample sizes, or collection dates in both studies (Moyes et al. 2017).
To further analyze the resistance found in these populations, kdr genotyping was performed for the 1016 and 1534 loci. The co-occurrence of these mutations has been already reported in different parts of the Americas, and their association with pyrethroid resistance was recently confirmed (Chen et al. 2019). Early detection of insecticide resistance mechanisms is critical for decision making in vector control programs (Corbel et al. 2019). High frequencies of the 1016I and 1534C mutant alleles were found in both populations in randomly sampled individuals (Table 2). A higher frequency for the 1016I allele was found for the Puntarenas population compared with Limon. Genotypes CC/VV, CC/VI, and CC/II were found; however, not one individual analyzed presented the wild-type allele for the 1534 codon, with frequencies of 1.00 for the 1534C allele. These results The kdr mutation V1016I was originally reported as a widespread mutation across many Latin American Ae. aegypti populations, including a population from Guanacaste province, Costa Rica (Saavedra- Rodriguez et al. 2007). Although functional conformation experiments do not confirm that the V1016I mutation  alone can cause pyrethroid resistance, many studies have correlated the presence of this mutation to phenotypic resistance (Deming et al. 2016, Brito et al. 2018, Maestre-Serrano et al. 2019, Ryan et al. 2019. Most recently, the V1016I mutation, in addition to the F1534C, was found to have an additive effect in resistance to both Type I and II pyrethroids (Chen et al. 2019). The F1534C mutation has been strongly associated with DDT and pyrethroid resistance and confirmed by functional experiments in Xenopus oocytes (Du et al. 2013, Hirata et al. 2014. Higher levels of pyrethroid resistance in mosquitoes have been observed in populations with the co-occurrence of the mutations V1016I and F1534C in other countries (Aponte et al. 2013, Kawada et al. 2014, Chapadense et al. 2015, Dusfour et al. 2015, Yaméogo et al. 2019. In Puntarenas, 86% of the individuals analyzed presented this co-occurrence, and 62% of the individuals from Limon (Table 2). Based on recent evidence, fixation of F1534C mutation in the populations analyzed represents a major threat to the insecticide-based vector control in Costa Rica (Vera-Maloof et al. 2015). Because this fixated mutation could favor the emergence of other kdr mutations, like the V1016I or the V410L, multiple kdr mutations could cause an additive effect increasing the levels of insecticide resistance and decreasing the efficiency of vector control programs (Li et al. 2015, Plernsub et al. 2016, Saavedra-Rodriguez et al. 2018. To observe an association between the presence of kdr mutations mentioned above and their phenotypic status, resistant and susceptible mosquitoes were analyzed. These analyses showed contrasting results between both populations, in which Puntarenas presented a clear association with the presence of the V1016I mutation and pyrethroid resistance to permethrin and deltamethrin, with a higher association to deltamethrin compared with permethrin resistance (Table 3). Limon, on the other hand, did not show this association. Cheng et al. (2019) confirmed that the F1534C mutation alone only confers resistance to Type I pyrethroids, such as permethrin, however, double homozygous mutant (CC/II) confers resistance to both, Type I and II pyrethroids (deltamethrin). This could explain the higher association of the V1016I presence with deltamethrin compared with permethrin. Due to the fixation of the 1534C allele, association analysis could not be performed with this locus. However, in agreement with Linss et al. (2014), the V1016I mutation was not present without the F1534C mutation, which also coincides with the sequential evolution of F1534C and V1016 hypothesis proposed by Vera-Maloof et al. (2015) and supported by Cheng et al. (2019).
Notably, our analysis showed that the Puntarenas population was not in HWE for the 1016 locus; in contrast, the Limon population was in HWE. This could indicate an external pressure to the Puntarenas population, which could explain the higher frequency of the mutant allele compared with Limon. Even though the insecticides and vector control strategies used in Costa Rica are the same for all regions, there is a marked difference in the results for the study sites above. This difference could be caused due to the ovitrap locations. In Puntarenas, the traps were located in a highly touristic site populated by several hotels, whereas the traps from Limon were located inside the houses of a residential neighborhood. Because private companies try to decrease mosquito density to accommodate tourist's demands, the Puntarenas population analyzed might be more exposed to insecticide pressure than Limon. Also, the null association between kdr mutations and pyrethroid resistance in the Limon population suggests that other mechanisms, such as metabolic resistance, may be involved as has been observed in other populations (Dusfour et al. 2015, Seixas et al. 2017. Overall, a strong association between the V1016I mutation and pyrethroid (Types I and II) resistance was found in all analyses for Puntarenas. Based on our results, Puntarenas is under stronger pressure, leading to the emergence of higher levels of pyrethroid resistance. In Limon, the lower levels of resistance limited the analysis of genotype/phenotype association, although HWE analyses seem to indicate that there is no existing pressure affecting the V1016I mutation. To accurately assess the insecticide resistance status in Costa Rica, more sites should be sampled and analyzed. The co-occurrence of the F1534C and V1016I mutation could be spread across Central America, making an urgent call to fill the gap of our knowledge about the insecticide resistance status and mechanisms in the Central American populations of Ae. aegypti.

Supplementary Data
Supplementary data are available at Journal of Medical Entomology online.