Kdr genotyping in Aedes aegypti from Brazil on a nation-wide scale from 2017 to 2018

Insecticide resistance is currently a threat to the control of Aedes agypti, the main vector of arboviruses in urban centers. Mutations in the voltage gated sodium channel (NaV), known as kdr (knockdown resistance), constitute an important selection mechanism for resistance against pyrethroids. In the present study, we investigated the kdr distribution for the Val1016Ile and Phe1534Cys alterations in Ae. aegypti from 123 Brazilian municipalities, based on SNP genotyping assays in over 5,500 mosquitoes. The alleles NaVS (1016Val+ + 1534Phe+), NaVR1 (1016Val+ + 1534Cyskdr) and NaVR2 (1016Ilekdr + 1534Cyskdr) were consistently observed, whereas kdr alleles have rapidly spread and increased in frequency. NaVS was the less frequent allele, mostly found in Northeastern populations. The highest allelic frequencies were observed for NaVR1, especially in the North, which was fixed in one Amazonian population. The double kdr NaVR2 was more prevalent in the Central-west and South-eastern populations. We introduce the ‘kdr index’, which revealed significant spatial patterns highlighting two to three distinct Brazilian regions. The 410L kdr mutation was additionally evaluated in 25 localities, evidencing that it generally occurs in the NaVR2 allele. This nationwide screening of a genetic mechanism for insecticide resistance is an important indication on how pyrethroid resistance in Ae. aegypti is evolving in Brazil.

and 1534 Na V sites were obtained. All evaluated population displayed at least one kdr allele. Six genotypes were substantially observed, all including the S, R1 and R2 (wild-type, kdr in 1534 only and kdr in 1016 + 1534, respectively) haplotypes. The other three genotypes composed by the kdr R3 haplotype (kdr in 1016 site only) accounted for 0.1% of the total samples ( Fig. 1 and Table 1). Therefore, the R3 haplotype was not considered in further analyses.
The most frequent Ae. aegypti genotypes was R1R2, followed by the double kdr homozygote R2R2, with medians of 29.0% and 24.4%, respectively (Table 1). The genotype frequencies per population in their respective geographic region are displayed in Fig. 1. Concerning allelic frequencies, medians were 13.3%, 37.8% and 45.3%, respectively for Na V S, Na V R1 and Na V R2, evidencing the high prevalence of kdr alleles (83.1%) in Brazil. The wild-type allele Na V S was absent from 17 of the 123 evaluated populations, while kdr Na V R1 was present in all localities, ranging from 7.1% (Corumba/MS) to 100% (São Gabriel da Cachoeira/AM) and the kdr Na V R2 ranged from absent (Tucurui/PA, Altamira/PA and São Gabriel da Cachoeira/AM) to 89.3% (Corumba/MS) (Fig. 2). The kdr allelic frequencies were regionalized in the country, as a reflection of their genotypic distribution, as follows: Na V R2 predominated in the Central-West, Southeast and South regions, while Na V R1 was the most frequent allele in the North (Fig. 3). The detailed genotypic and allelic frequencies for each locality are available in the Supplementary Table S1.
The spatial distribution kdr index value analysis concerning the geographic coordinates of their respective localities revealed statistically significant spatial patterns, i.e. the existence of clearly identifiable Brazilian regions with different predicted Ae. aegypti resistance status. A total of 55 candidate models, related to different (a priori) inter-city "links" structures, were generated (see Supplementary Text S1 for methodological details and intermediary results). Figure 5 displays the most significant model, first evidencing a large scale spatial differentiation in the Northeast-Southwest direction, separating Northeast (and to a lesser extent, North) localities from the rest of the country (areas A11 and A12, respectively, represented in Fig. 5a and, second, a three-mode differentiation in the same direction, separating: (i) North and Northeast localities; (ii) a strip of cities oriented  www.nature.com/scientificreports/ and 1534C (Na V R2 allele) and that 1534C also occurs alone (Na V R1 allele). Nevertheless, we cannot exclude the possibility that the underrepresented genotypes could be composed by additional haplotypes circulating under low frequencies. The complete dataset displaying the numbers and frequencies for each SNP in all populations is presented in Supplementary Table S2.  www.nature.com/scientificreports/

Discussion
The molecular nation-wide monitoring of the main pyrethroid resistance marker in 123 Ae. aegypti populations from Brazil, a country comprising 8.5 million km 2 , where the 26 federal states are infested with this vector, is presented herein. To the best of our knowledge, this is the largest kdr survey performed using samples collected in the time span of one year. The regionalized profile of allelic frequency distribution indicates that this sample size is enough to represent the whole country. The frequency of kdr alleles was very high in Brazil (median of 83.1%), representing an increase of 27.8% since a previous survey was performed, evaluating 29 localities where mosquitoes were collected between 2009 and 2012 14 . The kdr Na V R1 allele (1016Val + 1534Cys kdr ) was present in all Brazilian populations evaluated here, and fixed in one of them. This was probably the first kdr allele to spread in Latin America 20 and is present in at least two haplotypes of independent phylogenetic origins 17 . kdr Na V R2 (1016Ile kdr + 1534Cys kdr ) confers a higher level of resistance to pyrethroids compared to Na V R1 17 , which may explain why Na V R2 frequency is rapidly increasing compared to the wild-type Na V S and kdr Na V R1 alleles, as observed in Ae. aegypti populations from São Paulo 9 . It is interesting to note that, out of the 17 populations where the wild-type Na V S was absent, 12 were from the Amazonian region. In many of these localities, chemical control is not only exerted to combat Ae. aegypti, but also targeting other important vectors, such Anopheles mosquitoes, triatomines and phlebotomines, against which pyrethroids are the main choice inside houses and in the peri-domicile area. Therefore, the scenario observed herein may be explained by this constant selection pressure.
We also presented the application of a 'kdr index' which considers the kdr genotype frequencies weighed by their respective knockdown time resistance ratio (KdT RR 95 ). The lowest kdr indexes were predominately found in populations from the Northeast, where the kdr Na V R2 allele was almost absent from the previous survey 24 and is now present in all states in that region, with frequencies ranging from 10.2% (Jardim do Seridó, RN) to 83.8% (Itabaiana, SE). In fact, in a previous study assessing resistance to the pyrethroid deltamethrin amongst the 13 populations with pyrethroid resistance ratios (RR 95 ) below 10, reported that seven were located in the Northeast region. In addition, only two populations in the Northeast were amongst the 24 highly resistant populations (RR 95 > 10) 6 .
The study of the spatial distribution of resistance markers can aid in identifying underlying resistance emergence processes in surrounding areas. In fact, these determinants can exhibit significant "spatial patterns", due to their strong dependency to space, climate, land cover and land use, active and passive mosquito mobility between connected cities and regional strategies for vector control, among others. These patterns may display spatial variations at different scales. As a consequence, a spatial distribution of the resistance level results of the mixture of these different multi-scale spatial patterns is observed. The decomposition of the observed level of resistance into various independent and significant spatial patterns can, consequently, not only exhibit spatial clusters of cities that present significantly high (or low) resistance levels, but also highlight spatial patterns typical of phenomena that may, therefore, be considered as potential determinants. In the present study, we tested a series of analyses considering the spatialized kdr index values, revealing two to three well defined clustered www.nature.com/scientificreports/ populations in the country. This evidenced that the genetic trend of higher or lower pyrethroid resistant (the kdr index) is not random in Brazil. The populations placed in the A12 or A22 areas, in the models with two or three clusters, respectively, are more likely to be resistant to pyrethroids than the rest of the country, considering the target site mechanism. On the other hand, specific conditions may lead to diverse insecticide resistance patterns, including distinct kdr allelic frequencies, as observed in Ae. aegypti populations from five neighborhoods in Merida, Mexico 30 . New studies including spatial patterns details will further knowledge on the dynamics of insecticide resistance distribution and expansion in natural populations. The most frequent genotype in Ae. aegypti from Brazil was R1R2, which was heterozygous for the alleles Na V R1 and Na V R2, in 28.9% of all evaluated samples. The second most frequent, at 24.4%, was the double kdr R2R2, homozygous for Na V R2, followed by R1R1 in 19.9%, homozygous for Na V R1. The sum of all "resistant genotypes" 27 indicated that 73.2% of all Ae. aegypti from Brazil evaluated herein would be resistant to pyrethroid, at least considering the target-site resistance mechanism. It is worth noting, however, that metabolic resistance might also contribute, alone or in synergism with kdr, to confer pyrethroid resistance in the natural populations of this country 6 . Although the organophosphate malathion began replacing pyrethroids in official campaigns oriented by the Brazilian MoH in 2009 31 , the selection pressure of pyrethroids is still heavily present in Ae. aegypti due to household aerosol insecticide products 9,32,33 . Experiments conducted in Merida, Mexico, indicate that 87% of households regularly used pyrethroid-based products against mosquitoes 33 and that Ae. aegypti field populations were resistant to commonly household used products. In addition, an increase in kdr 1016Ile was associated with the employment of pyrethroid surface sprays in houses 34 .
kdr mutations, or at least the kdr Na V R2 allele, have a high fitness cost and tend to decrease in insecticidefree environments 25 . In a very practical example, the first Ae. aegypti lineage infected with Wolbachia released in Rio de Janeiro was actually a result of backcrosses between an original infected Australian lineage and a natural Brazilian population, therefore resulting in a Wolbachia-infected lineage resistant to pyrethroids. However, during several generations in the laboratory awaiting governmental licenses for release, kdr frequencies decreased and the lineage became susceptible 32 . As residents were spraying pyrethroid-based aerosols in the study area, www.nature.com/scientificreports/ Wolbachia dispersion did not occur as expected and the first release failed. Thereafter, a new lineage was again backcrossed with field mosquitoes in order to guarantee similar conditions to the well-adapted field population, which then allowed Wolbachia to expand and stabilize in the study area. The kdr frequency is now monitored in the laboratory lineage and in the target field population in the World Mosquito Program 32 .
The kdr SNP V410L was first detected in Ae. aegypti samples from Northern and Southeastern Brazil 29 , also observed in samples from Mexico since 2002 21 , where it is currently highly disseminated 35 . The V410L SNP was detected in several of our evaluated samples, with strong association of the mutant variation (410Leu) with 1016Ile (and, therefore, with 1534Cys), i.e. the triple kdr allele. In a few samples however, 410L was dissociated from these other kdr mutations, as evidenced by individuals genotyped as LL + VI + FC, VL + VV + CC and VL + II + CC (all in the order 410 + 1016 + 1534), which account for only 1.1% of the total. Similar results were  www.nature.com/scientificreports/ also observed in samples from Mexico, where mutations in 410 and 1016 were always associated, or, otherwise, dissociated in very low frequencies 35 . In order to save consumables in a broad kdr surveillance in Ae. aegypti in Latin American countries, where resources are limited, we recommend that the 1016 site may be genotyped first, and for samples genotyped as 1016 V/V or V/I, a second reaction for the 1534 SNP should be performed. On the other hand, when the samples are detected as 1016 I/I, their genotype are likely to be R2R2, i.e. homozygous for the three kdr mutations. Thus, it would not be necessary to genotype the 410 site at all. For academic purposes however, it would be interesting to track possible fluctuations in the frequency of these rare genotypes over time.
New methods for controlling Ae. aegypti or at least making the species "immune" to arbovirus infections are being tested or have already been implemented as actual control methods 36 . The use of insecticides however may be maintained, due to their potential of rapidly decreasing the density of a targeted population, unless resistance to the applied compound is present. New products with active neonicotinoid-class ingredients have been recently approved by WHO PQT-VC, namely clothianidin, as an indoor residual spray (IRS), and imidacloprid, as an ultra-low volume (ULV) compound 37 . Although these neonicotinoids are slower-acting than pyrethroids and organophosphates, resistant populations to these insecticides seem to not present cross-resistance 38 , indicating them as promising alternative adulticides. In Brazil, the most recent nation-wide Ae. aegypti IR monitoring detected several populations resistant to malathion 39 . The Brazilian MoH is currently planning the implementation of an alternative compound comprising a mixture of pyrethroids and clothianidin to be sprayed inside houses, in addition to another compound comprising a mixture of pyrethroids and imidacloprid as an ULV for outdoor areas. On the other hand, a deeper discussion should consider environmental issues regarding neonicotinoids in non-target organisms 40 . In meanwhile, new arbovirus case records in 2019 in Brazil consisted of over 1.5 million dengue cases, 130 thousand chikungunya cases and 10 thousand Zika cases 1 . In Brazil, as well as in several other countries with similar climate and urban landscape conditions, the implementation of new vector control tools is paramount.

Material and methods
Collections were performed in the context of the national Brazilian insecticide resistance monitoring program for Ae. aegypti between 2017 and 2018. Field agents installed eggtraps randomly distributed according to the total number of houses in each municipality, as follows: 100, 150, 200 or 300 traps for < 50, 50-200, 200-500 and > 500 thousand houses, respectively. The installation and collection methodology of the eggtraps is described elsewhere 39 . A presential meeting was organized with field-personnel representatives of each state and a video with step-by-step procedures was launched (https ://www.youtu be.com/watch ?v=2w89k agSOK M) in order to ensure that collections were made as homogeneously as possible. The palettes were sent to the laboratory, where the eggs were stimulated to hatch and larvae were raised up to adults in order to produce an F1 colony of each locality, as previously described 21 . After 4-5 days in the cages, enough time for copulation, a sampling of this F 0 generation of each population, around 45 insects, preferentially males, were removed and saved for kdr genotyping. These mosquitoes were maintained in 80% ethanol solution or cryopreserved prior to DNA extraction. The DNA was extracted from single mosquitoes using the FastDNA Spin (MP Biomedicals) or Nucleo Spin Tissue (Macherey-Nagel) kits, according to their manufacturers' protocol. The eluted DNA was diluted to 10 ng/µL in extra-pure water and cryopreserved until use.
We performed independent genotyping reactions for each kdr site based on a qPCR approach using the Custom TaqMan SNP Genotyping Assay (ThermoFisher) (see Table 3 for the primer and probe sequences for each assay). Reactions consisted of 1X TaqMan Genotyping Master Mix (ThermoFisher), 1X of the respective Custom TaqMan SNP Genotyping Assay, 20 ng of DNA and ultra-pure water q.s. 10 µL, run in a QuantStudio We ranked the populations in relation to their predicted level of resistance to pyrethroids, based on the kdr index, consisting of the sum of kdr genotype frequencies weighed with their respective resistance ratios, which were previously obtained in a knockdown time test (KdT RR 95 , Brito et al. 27 . The KdT RR 95 values for the genotypes were SS (1), SR1 (2.4), SR2 (1.7), R1R1 (4.6), R1R2 (5.4) and R2R2 (6.7). The formula for the kdr index of a population is: Table 3. Primer and probe sequences for the SNPs V410L, V1016I and F1534C kdr in Aedes aegypti. a Identification of the customized TaqMan SNP Genotyping Assay (ThermoFischer). www.nature.com/scientificreports/ In order to test the existence of significant multi-scale spatial patterns (according to Moran's Index of spatial autocorrelation) regarding the spatial distribution of the kdr frequencies, we considered the aforementioned kdr index and the geographic coordinates of each locality in a series of spatial statistic tests, applying a principal coordinate analysis of neighbor matrices and Moran's eigenvector maps 41 , as detailed in Supplementary Text 1.

Data availability
A table with the genotyping data of all populations can be accessed in the supplementary files. Further data and materials are available from the corresponding author upon reasonable request.