The positive spatial autocorrelation in counts of the IICC, pyrethroid-resistant genotype at the sampled sites indicates that neighboring sites tend to have similar IR profiles for this species, particularly when sites are approximately 20 km away each other. In a scenario where the status of an Ae. aegypti population is unknown at one site, a vector control district could predict that it will be typically similar to neighboring sites with a known kdr genotype frequency. However, this estimate would merely be a starting point if data on resistance status are not immediately available and would ideally not replace an actual assessment of susceptibility. This is especially true given the limited intensive, local-scale sampling included in this dataset, which was collected with the aim of capturing statewide patterns. For example, while collections in Miami-Dade County were conducted at many sites, most of these collections were pooled to centroids. Across these pools, the average distance between sites ranged from 40 meters to 4.4 kilometers. Within more dispersed clusters, genotype frequencies may have varied across sites, but this variation at the local scale cannot be examined with this dataset because collected mosquitoes were pooled to the cluster centroid. Similarities in IICC counts decrease and approximate a random distribution at the scale of 120 km, indicating that sites from neighboring counties would not be informative in predicting IICC frequencies. For future work in estimating IICC prevalence throughout the state, interpolation methods would likely be informative. However, other studies have found reduced accuracy in multiple interpolation methods when sampling sites are highly clustered, as is the case with this dataset (42). To address this concern and employ interpolation in this scenario, more evenly distributed sampling would need to be conducted, particularly in the central parts of the state that are not currently well represented and generally have a more rural landscape that is distinct from the well-sampled, urban coasts.
We included SaTScan analyses with both circular and elliptic scanning window shape, since the sampled points in this study were predominantly located in coastal areas, which reflects the natural distribution of Ae. aegypti in FL, and an elliptic shape may capture this configuration more adequately. While the results from these two analyses were similar, each approach identified at least two significant clusters that that were unique to that method. This confirms previous findings that suggest considering multiple graphical representations of clusters in SaTScan can yield clusters that would not otherwise be detected (29). Further work in this area would benefit greatly from collecting longitudinal data on kdr genotype frequencies to identify the temporal variation in the identified clusters using the space-time scan statistic available in SaTScan, allowing for assessment of the stability or seasonality of the patterns identified in this study.
The areas identified as significant clusters of the IICC genotype were in coastal cities, including New Port Richey, Naples, and Miami. Each of these clusters comprised four or fewer sites, indicating that these occurrences of high IICC genotype frequency may be the result of neighborhood or household-scale selective pressures. Similar work in Yucatàn State, Mexico found significant differences in IICC frequencies between city blocks (42), and it has been shown that typical household-level insecticide use can result in selection for resistance in Aedes aegypti (44). Clusters of the most susceptible genotype, VVFF, which is homozygous susceptible at both loci, were identified with centers in St. Johns County and Miami-Dade County. In the case of Miami-Dade County, this cluster could represent a localized refuge wherein mosquitoes are somehow protected from the regular applications of pyrethroid treatments, despite heavy use throughout the rest of the county. This persistence of the wild-type genotype, even if only at limited refugia, has important implications for the potential to re-establish susceptibility. Given the fitness cost of pyrethroid resistance to mosquitoes, susceptible mosquito populations would gain an advantage if the selective pressure of pyrethroid applications is removed (45).
The top-ranked beta regression models revealed a complex relationship between vegetation and insecticide resistance in Ae. aegypti in this system. Based on field trials that found lower Ae. aegypti mortality in caged trials in areas with high vegetation density (23), we predicted lower frequencies of the IICC genotype in areas with high vegetation density. There was a negative association between percent non-tree vegetation cover and IICC genotype frequency in three of the original twenty candidate models, all of which were ultimately in the best performing spatial lag models. This indicates that areas with dense ground vegetation may shelter Ae. aegypti from insecticide applications and reduce the selective pressure driving insecticide resistance. However, there were significant positive associations between the variables of January or October EVI, January or July LAI, and percent tree cover, and the outcome of IICC genotype frequency. This means that areas with healthy, dense canopy cover are more likely to have pyrethroid-resistant Ae. aegypti populations. This could be due to overall higher abundance of multiple mosquito species in shaded, sheltered areas (46,47), leading to a response of more liberal applications of insecticides by control districts, pest control companies, and private landowners, and overall intense selection pressure. Understanding the exact nature of the impact of vegetation on insecticide application efficacy and the development of resistance will likely require experimental field trials or fine-scale sampling across several landscape configurations.
The negative association between the intensity of pyrethroid use by control district and the frequency of the IICC genotype was unexpected. This outcome could arise due to limitations in the dataset on insecticide use, which was derived from Chemical Activity Reports submitted by mosquito control districts. In these reports, figures on the total amount of each product used and the number of acres treated were reported, but there was no delineation of the areas within each district that were treated. This means that collections sites used in this study could potentially be outside of the areas that were treated regularly and intensively with pyrethroids. Additionally, districts that have enough financial resources to purchase and apply large amounts of pyrethroids may also be more likely to have the ability to incorporate larvicides, biological control, or public education efforts into their integrated pest management plans, thus mitigating the selective pressure of pyrethroid-use in some areas of the district. Finally, this dataset only included reports from mosquito control districts, while insecticide applications in this and other states are frequently conducted by private pest control companies or landowners (48). Despite these limitations, these data represented the most accurate and current representation of insecticide use in the state, meaning it would be inadvisable to exclude them from consideration for these models. Further modelling of insecticide resistance would benefit from more detailed information on the locations and timing of insecticide treatments.
The results from the original beta regression models identified statistically significant relationships between landscape factors and the outcome of IICC frequency. However, these models only explained a portion (on average, 32%) of the variation in this response variable. The spatial lag models included information on neighboring IICC frequencies, as well as the original landscape factors, to explain variation in the response variable (49). The spatial lag models improved the overall fit and explained a greater amount of the variation present in IICC frequencies than the non-spatial beta regression models, with a maximum R2 value of 0.49. In the spatial lag models, the coefficients of the original variables decreased and the p-values associated with them increased, with some of the associations no longer being statistically significant. This indicates that including IICC frequencies of neighboring sites can dramatically improve our ability to estimate the IICC frequency at a given location and may provide more information than the landscape and insecticide-use variables originally considered. The significant improvement in the models observed with the addition of the spatial lag also indicates that there is likely a diffusion process occurring, meaning immigration between neighboring sites may explain some of the observed spatial dependence (41). While knowledge of this spatial dependence is useful for predicting IICC frequencies if information on nearby sites is available, the original models indicate that information on landscape factors and insecticide use can explain some of the underlying variation present.
Ultimately, the best performing model explained roughly half of the variation in the frequencies of the IICC genotype. Several factors likely contribute to the remainder of that variation. Data on the history insecticide applications and resistance status throughout the county would likely improve these models. For example, if the IICC genotype had reached fixation in a population due to regular applications of pyrethroids in past years, susceptible variants would then only arise due to gene flow, which may not occur if the population is geographically isolated from susceptible strains. Additionally, information on the intensity of insecticide use by individuals and private companies would likely be informative. In a recent study across three regions in North Carolina, researchers found that approximately 31% of survey respondents applied pesticides to control mosquitoes on their properties (50). If Florida residents are applying insecticides at a similar rate, this likely constitutes a selective pressure. To further improve the fit of these models, it would be valuable to restrict sampling to a more limited time period, since variation in IICC frequency has been shown to fluctuate seasonally (51), and include more sampling sites, which would allow the consideration of potential interactions between the explanatory variables.
This study is a demonstration of the potential for spatial analysis methods to explain and estimate frequencies of kdr frequencies in Ae. aegypti populations, with important implications for vector control. While kdr frequencies and the phenotype of pyrethroid resistance have been shown to vary spatially and temporally (51-53), intensive, constant monitoring of resistance is not logistically feasible. This leaves a need for tools to fill in gaps in knowledge both in time and space, although this project only focused on space. Important takeaway points from the perspective of a vector control organization include the following: First, populations within 20 kilometers of each other are likely to have similar IICC frequencies, particularly if these sites have similar landscape features. This means that if resistance is detected at one site, neighboring locations should be assumed to have similar measures, and alternative control strategies should be deployed. A second important point is that some areas in Florida have higher than expected frequencies of the wild-type genotype, meaning these populations are still susceptible to pyrethroids. These control districts may have developed programs that curtail the development of resistance. This highlights the importance of sharing information between control districts on integrated pest management practices that are effective, sustainable, and culturally acceptable within the context of Florida. Finally, we found that pyrethroid use by control districts was negatively associated with the outcome of high IICC frequencies. While there were limitations to this dataset, as discussed above, this indicates that other forms of insecticide treatments, either by landowners or commercial companies, may contribute to selection for kdr alleles, particularly at local scales, as has been shown elsewhere (54). Vector control districts should be aware of the extent to which these largely unregulated applications may be impacting the populations they are aiming to control, and further work should be done to identify and quantify those impacts.
Importantly, all of the above analyses relied on frequencies or counts of the genetic mutations in the voltage-gated sodium channel associated with pyrethroid resistance, rather than the actual phenotype of resistance. While other mechanisms related to the outcome of resistance in this and related species have been identified (16), variation in resistance status in the Ae. aegypti populations studied here is highly correlated with frequencies of these mutations, as was shown in Estep et al. 2018 (18), where the authors found a Pearson correlation coefficient of 0.905 between frequencies of the IICC genotype and the permethrin resistance ratio (3).