Associated patterns of insecticide resistance in field populations of malaria vectors across Africa.

The development of insecticide resistance in African malaria vectors threatens the continued efficacy of important vector control methods that rely on a limited set of insecticides. To understand the operational significance of resistance we require quantitative information about levels of resistance in field populations to the suite of vector control insecticides. Estimation of resistance is complicated by the sparsity of observations in field populations, variation in resistance over time and space at local and regional scales, and cross-resistance between different insecticide types. Using observations of the prevalence of resistance in mosquito species from the Anopheles gambiae complex sampled from 1,183 locations throughout Africa, we applied Bayesian geostatistical models to quantify patterns of covariation in resistance phenotypes across different insecticides. For resistance to the three pyrethroids tested, deltamethrin, permethrin, and λ-cyhalothrin, we found consistent forms of covariation across sub-Saharan Africa and covariation between resistance to these pyrethroids and resistance to DDT. We found no evidence of resistance interactions between carbamate and organophosphate insecticides or between these insecticides and those from other classes. For pyrethroids and DDT we found significant associations between predicted mean resistance and the observed frequency of kdr mutations in the Vgsc gene in field mosquito samples, with DDT showing the strongest association. These results improve our capacity to understand and predict resistance patterns throughout Africa and can guide the development of monitoring strategies.

Malaria prevalence in Africa has declined over the last 20 y and these gains have, in large part, been due to the widespread implementation of insecticide-based vector control measures (1). In the same time period we have also seen increases in insecticide resistance that have the potential to derail or even reverse the progress made in reducing malaria transmission (2, 3). It is therefore vital to understand the impact of insecticide resistance on malaria transmission.

Only 12 available insecticides, all insect neurotoxins, have WHO Pesticide Evaluation Scheme (WHOPES) approval for use in malaria vector control: 6 of these for insecticide-treated nets (ITNs), 3 for long-lasting ITNs (LLINs), and 11 for indoor residual spraying (IRS) (4). To steward this limited resource, WHO encourages countries to develop insecticide resistance management plans based on local data; however, these data are sparse or nonexistent in many places. Insecticide resistance is often considered simplistically as a binary (dead or alive) response to insecticide exposure, but the resistance phenotype can be generated by a plethora of different mechanisms resulting in populations that differ in the spectrum and level of resistance to different insecticides. Moreover, resistance can vary temporally and spatially within and between insecticides and species. To understand the operational significance of resistance, we must consider the patterns of variation in resistance to each insecticide recommended for use. Understanding cross-resistance between insecticides, especially across different classes, is also crucial for developing effective insecticide resistance management strategies.

Studies of the underlying mechanisms of resistance suggest that common patterns of phenotypic resistance to different insecticides may occur in field populations of mosquitoes. However, the quantitative impact of combinations of mechanisms found in the field remains poorly understood. Three categories of mechanism confer resistance to neurotoxic insecticides in malaria vectors: alterations to metabolic genes or pathways, target site mutation, and cuticular thickening. Metabolic resistance resulting primarily from the amplification or up-regulation of detoxification enzymes (especially esterases, P450 monooxygenases, and glutathione S- transferases) occurs commonly and can confer high levels of resistance (5–7). Some enzymes have been linked to resistance to a specific insecticide or class (8, 9) whereas others confer resistance across insecticide classes (5, 7, 10). Changes to insecticide target sites are also prevalent and frequently associated with phenotypic resistance (11–13). Acetylcholinesterase target site alteration caused by a mutation to the Ace-1 gene typically confers cross-resistance to organophosphate and carbamate insecticides (5, 14). Cross-resistance to pyrethroids and the organochlorine DDT arise from knockdown resistance (kdr) mutations in the Vgsc gene which encodes the para voltage-gated sodium channel target site (15). There is currently less evidence for cuticular thickening as a mechanism of resistance (16–18).

Resistance to multiple insecticides within a chemical class is commonly noted but rarely quantified and differences in resistance to different insecticides within a class are also seen (13). Detection of resistance relationships is impeded by the high measurement error of discriminating dose bioassays that estimate the prevalence of resistant phenotypes. Additional complications arise from the co-occurrence, but uncertain impacts, of multiple mechanisms of resistance within mosquito populations and within individual mosquitoes (12, 19–21) and by poor understanding of the influence of environmental factors (22). We do not yet know how these varied resistance mechanisms—with their differing geographical distributions (23–25)—combine to impact variation in resistance within and among field populations or whether the cross-resistance patterns seen in laboratory studies translate into common geospatial patterns of phenotypic resistance to different types of insecticide in the field. If predictable relationships between resistance across different insecticides in field populations can be found, this will improve our ability to assess and manage the impact of resistance on vector control.

To date, the majority of tests of phenotypic insecticide resistance conducted on African malaria vectors have used Anopheles gambiae complex samples. Using data on the prevalence of insecticide resistance in samples collected from 1,183 locations in 38 African countries, we applied a Bayesian geostatistical modeling approach to analyze the relationships in resistance across the insecticides most commonly used in malaria vector control. For insecticides within the pyrethroid and organochlorine classes, we found consistent relationships in the prevalence of resistance across different insecticide types over a large part of sub-Saharan Africa. Making use of relatively abundant data on the frequency of well-characterized mutations in the Vgsc gene in field-collected samples from the A. gambiae complex, we demonstrated clear associations between the average prevalence of the resistance phenotype for DDT and pyrethroid insecticides and the frequency of Vgsc mutations. These results develop our capacity to understand and predict phenotypic resistance patterns throughout Africa for commonly used insecticides.

Results

Large-Scale Resistance Phenotype Relationships Across Insecticides.

We investigated associations between resistance phenotypes across different types of insecticides using the results of WHO susceptibility bioassays conducted on mosquito samples from sub-Saharan Africa over the period 2000–2015 (26), using the proportional sample mortality from the bioassay. The distribution of the sample collection locations is heterogeneous across space and time and varies across the different insecticides tested (Fig. 1 and ). To accommodate this sampling heterogeneity we applied a Bayesian geostatistical modeling approach. Due to the lack of samples from central Africa, we first partitioned the data into two separate spatial regions covering western and eastern parts of the continent (Fig. 1 ). We then assessed resistance associations across different insecticide types by comparing a linear model of coregionalization (LMC) that allows interactions between resistances across each insecticide type with a model where resistances to the different insecticides are assumed to be independent. Two measures were used to compare the models: (i) the Watanabe–Akaike information criteria (WAIC) (27, 28) and (ii) the root-mean-square error (rmse) obtained from K-fold out-of-sample validation, calculated separately for each insecticide type. We performed independent model fits for the data subsets within the West and East spatial regions ().

Fig. 1.

The spatiotemporal distribution of the sample collection locations for insecticide resistance bioassays included in our dataset. Rectangles enclose the West and East regions considered in our analysis. (A) Pyrethroid (Py) bioassays. (B) Organochlorine (Och) bioassays. (C) Carbamate (Ca) bioassays. (D) Organophosphate (Oph) bioassays. (E) The number of bioassay records for each time period. The keys in A–D correspond to the insecticides shown in E. (F) The locations of sample collection used to calculate Vgsc allele frequency data. The mutations L1014F and L1014S present at each location are shown.

Associations within an insecticide class.

For pyrethroids, we jointly model resistance observations for the three most commonly tested pyrethroids: deltamethrin, permethrin, and λ-cyhalothrin (Fig. 1). A LMC that allowed interactions between resistances across all three insecticide types performed better than a model where these resistances did not interact according to both the WAIC (Table 1, row 1 and ) and the out-of-sample rmse (, row 1). This indicates that associations between insecticides explain variation in the resistance observations. The predictions of the mean proportional mortality across each location and time show strong associations between the prevalences of resistance to the three insecticides in populations from both West and East Africa (Fig. 2). The Pearson correlation coefficient r between predicted prevalence of resistance to deltamethrin and that to permethrin has a posterior mode of r = 0.91 [credible interval (CI) = 0.90, 0.93] and r = 0.94 (CI = 0.94, 0.95) for the West and East regions, respectively. For both regions, predicted prevalence of resistance to permethrin is typically higher than that to deltamethrin, but with a greater difference apparent in West Africa (Fig. 2). Among cyano-pyrethroids, higher prevalence of resistance to λ-cyhalothrin than to deltamethrin is seen in both West and East Africa [Fig. 2; r = 0.99 (CI = 0.99, 1.0) for the West region and r = 0.98 (CI = 0.980, 0.984) for the East region]. The predicted resistance relationship between permethrin and λ-cyhalothrin is approximately linear for both regions, although the East region shows stronger relative resistance to λ-cyhalothrin in areas with high resistance to both insecticides [Fig. 2; r = 0.90 (CI = 0.89, 92) for the West region and r 0.94 (CI = 0.94, 0.95) for the East region]. Importantly, despite minor differences, the three predicted resistance relationships are similar between the West and East regions, and there was no credible difference between the two regions in the parameters describing the strength of the interactions across insecticides (). This indicates that phenotypic resistance associations across different pyrethroids are relatively homogenous across sub-Saharan Africa.

Table 1.

The difference in the WAIC between a LMC and a model where resistances do not interact across different insecticides

	ΔWAIC (SE)
Insecticide included in model	West region	East region
D, P, L	−134.2 (22.7)	−154.0 (28.3)
Och, D, P, L	−110.8 (18.3)	−44.3 (16.9)
Ca, D, P, L	5.5 (4.3)	−3.0 (11.9)
Oph, D, P, L	32.5 (16.7)	−4.7 (3.6)
Och, Ca	−4.2 (5.9)	0.7 (0.4)
Och, Oph	7.9 (5.5)	−0.17 (1.1)
Ca, Oph	−9.9 (19.5)	8.1 (5.7)

A negative ΔWAIC value indicates that the LMC performs better and a positive value indicates that a model with no interactions across insecticides performs better. We consider ΔWAIC values that are lower in magnitude than (or approximately equal to) the SE to be inconclusive. Substantial differences are highlighted in boldface type. Independent models are fitted to the West and East regions. The insecticide bioassays included in the models are denoted as follows: bendiocarb (Ca), DDT (Och), deltamethrin (D), λ-cyhalothrin (L), permethrin (P), and organophosphates including fenitrothion for the West region and both fenitrothion and malathion for the East region (Oph).

Fig. 2.

Relationships between predicted mean proportional bioassay mortalities for the three pyrethroid types. (A) Deltamethrin vs. permethrin. (B) Deltamethrin vs. λ-cyhalothrin. (C) Permethrin vs. λ-cyhalothrin. Points show the predicted mean at each location and time for the West region (blue) and the East region (red). Color intensity indicates the width of the posterior CI of the predicted mean.

Associations across insecticide classes.

We investigated associations between resistance to DDT and to pyrethroids by extending the LMC developed for the three pyrethroid types to include interactions between DDT and pyrethroids. This LMC performed better than a model where resistance to DDT did not interact with resistance to pyrethroids according to both the WAIC (Table 1, row 2) and the out-of-sample validation based on withheld DDT resistance observations (, row 2). Predicted prevalence of resistance to DDT is much higher than that for the three pyrethroids (Fig. 3). The form of the relationships is similar between the West and East regions, although for the latter there are more outliers for which resistance to pyrethroids is higher (Fig. 3 and ). Correlations between predicted prevalences of resistance to DDT and to pyrethroids are lower for the East than for the West African region. For the West region the correlations between prevalences of resistance to DDT and permethrin, deltamethrin, and λ-cyhalothrin were r = 0.91 (CI = 0.9, 0.92), r = 0.91 (CI = 0.89, 0.92), and r = 0.92 (CI = 0.91, 0.93), respectively. For the East region these correlations were r = 0.65 (CI = 0.61, 0.7), r = 0.56 (CI = 0.51, 0.6), and r = 0.6 (CI = 0.56, 0.64), respectively. The predicted strength of the interaction between DDT and pyrethroids was also greater for the West compared with the East region ().

Fig. 3.

Relationships between predicted mean proportional bioassay mortalities for DDT and the three pyrethroid types. (A) DDT vs. permethrin; (B) DDT vs. deltamethrin; (C) DDT vs. λ-cyhalothrin. Points show the predicted mean at each location and time for the West region (blue) and the East region (red). Color intensity indicates the width of the posterior CI of the predicted mean.

We used the same approach to investigate associations between resistance to bendiocarb and pyrethroids. There was no evidence that a LMC that included interactions between resistance to the two insecticide classes performed better than a model that excluded these interactions (Table 1, row 3; , row 3; and ). In other words, statistical evidence for broad-scale carbamate–pyrethroid associations was absent for this pair of classes, despite local evidence from Anopheles coluzzii for cross-acting metabolic mechanisms (5). We also assessed associations between resistance to organophosphates and resistance to pyrethroids, using a modified LMC for the East region to account for the two different organophosphate insecticides, fenitrothion and malathion [; we note that due to the small number of available observations, we did not assess relationships in resistance between these two organophosphate insecticides (Fig. 1)]. There was no evidence that the LMCs performed better than models that excluded interactions between resistance to organophosphates and resistance to pyrethroids (Table 1, row 4 and , row 4); i.e., once again broad-scale resistance associations were absent (e.g., ). Finally we investigated relationships between resistance across the remaining pairs of insecticide classes, including organophosphate–bendiocarb, organophosphate–DDT, and bendiocarb–DDT resistance associations, using LMCs with only one interaction term. In all cases, our model comparison did not support choosing the LMC over a model that excluded interactions between resistance across different insecticide classes (Table 1 and , bottom three rows), again indicating a lack of broad-scale resistance interactions across insecticides.

Associations Between the Prevalence of Resistance Phenotypes and Allele Frequencies.

We investigated associations between our predicted values of the prevalence of insecticide resistance and the observed frequency of kdr mutations in the Vgsc gene in field samples of mosquito species from the A. gambiae complex. This is the only specific resistance mechanism for which substantial spatial and temporal genetic data are currently available and encompasses the two well-known Vgsc point mutations, L1014F and L1014S. We obtained 215 observations of Vgsc mutation prevalence in samples from West Africa and 101 observations from East Africa (Fig. 1) spanning 2005–2015 (). As expected, L1014F was the most common form of mutation in the West and L1014S in the East, although both forms were found in both regions (Fig. 1). We obtained estimates of the mean proportional mortality following exposure to DDT and the three pyrethroid types at the locations and times corresponding to our Vgsc mutation frequency observations using our predicted insecticide resistance surfaces (Figs. 2 and 3). We fitted ordinary least-squares (OLS) linear regression models to the predicted mean mortality corresponding to each insecticide, using the frequency of each type of Vgsc mutation as the two covariates (). These models were fitted to the total set of observations combined across the East and West regions (n = 316).

We found a significant association between the frequency of each Vgsc mutation and predicted mean resistance to DDT (Fig. 4; adjusted R2 = 0.61; ). There were also significant, but less strong, associations between the frequency of Vgsc mutations and mean resistances to deltamethrin, permethrin, and λ-cyhalothrin (). However, for these three pyrethroids, the linear regression models did not detect significant differences between the effects of the two mutations (L1014F and L1014S) on predicted resistance.

Fig. 4.

The mortality to DDT predicted by a fitted linear relationship with observed Vgsc allele prevalences (n = 316). Estimates of mean mortality following exposure to DDT obtained from our Bayesian geostatistical models (Fig. 3) were used as the response variable. Color intensity indicates the relative posterior precision of the mean mortality estimates. The diagonal line indicates equivalence between the data and the fitted values.

Discussion

Using data from 1,183 locations across Africa and statistical methods that accommodate spatiotemporal heterogeneity in the data sampling, we detected consistent relationships between phenotypic resistance to different insecticides from the pyrethroid and organochlorine classes, as well as associations with widely screened Vgsc target site mutations, across a broad spatial and temporal scale. Models for each pyrethroid and the organochlorine, DDT, performed better when data for multiple insecticides were included, indicating that associated patterns of resistance across insecticides improve prediction for individual insecticides. Pyrethroid insecticides are of particular importance because permethrin and deltamethrin are two of the three insecticides approved for use in LLINS and are the most commonly used insecticides in ITNs (4).

There were also clear relationships between the modeled prevalence of resistance to these pyrethroid insecticides and to DDT, but with greater variation. Laboratory experiments and local-scale field studies have shown that mutations to the Vgsc target site gene and some enzymes involved in metabolic resistance confer cross-resistance to both groups, which may explain the broad-scale field associations we detect (5, 7, 11, 13). Our demonstration of associations between the prevalence of resistance to DDT and to pyrethroids and the prevalence of Vgsc mutations in field-sampled mosquitoes certainly supports this hypothesis. The predicted associations between the pyrethroid phenotypes and the DDT phenotype were stronger in West Africa than in East Africa. A number of factors may explain differences between these regions. Each of the underlying mechanisms of resistance has a different geographical distribution, including regional differences in frequency of the two Vgsc 1014 alleles, and different combinations of mechanisms are found in different locations (23–25). In addition, there are differences in the geographical distributions of sibling species within the A. gambiae species complex (29) and these species differ in their resistance phenotypes and underlying mechanisms of resistance (11, 13, 21). While we think that these two explanatory factors, and their interaction, probably explain much of the predicted regional variation, there is greater prediction uncertainty for the East region because it covers a larger spatial area and has a sparser spatial coverage of sampling locations. There is also a lower proportion of observations where the prevalence of resistance is high in the East compared with the West region, so predicted resistance interactions are more difficult to quantify.

The same associations among prevalences of phenotypic resistance were not found when organophosphate and carbamate insecticides were included. This was unexpected because mutations in the Ace-1 gene are known to confer resistance to both organophosphates and carbamates where present in West Africa (5, 14). This result may reflect the lower power to detect associations for these classes because data volumes are lower () and perhaps most importantly, variation in mortality was less, especially for the organophosphates to which resistance remains limited. In addition, while organophosphate insecticides were pooled for this analysis, cross-resistance between insecticides within this class has been found to be variable, which could cause variable cross-class interactions in field populations (30–32).

Our results demonstrate broad-scale field population-based associations between Vgsc allele frequency and phenotypic resistance to DDT and pyrethroids. Local studies of the influence of different Vgsc alleles on levels of resistance in field populations have not shown consistent results (33), but our findings provide evidence that Vgsc allele frequencies are a useful partial diagnostic for resistance. Our findings are consistent with genetic studies that have shown that Vgsc mutations have emerged and spread as insecticide resistance has emerged and spread (10, 23–25, 34). The association found here is strongest with DDT resistance, which is consistent with field studies that have found stronger links between Vgsc alleles and DDT resistance compared with pyrethroid resistance (13, 35) and stronger relationships between metabolic resistance enzymes and pyrethroid resistance compared with DDT resistance (36, 37). Importantly, our results demonstrate the potential utility of genetic markers. Given an appropriate sampling procedure, use of DNA diagnostics to predict resistance phenotypes could provide a means of increasing resistance data resources to support monitoring programs.

An important limitation of our results is that they consider resistance in the A. gambiae species complex as a whole. The available bioassay data are predominately linked to the species complex because live mosquitoes are required and sibling species can be distinguished only by molecular analysis. As noted above, individual species are known to differ in the underlying mechanisms of resistance, phenotypic levels of resistance, and their geographical distributions. These factors can be partitioned only when sufficient single-species data are available for analysis, and we advocate, wherever possible, that field studies should type the mosquitoes tested.

In this study, we have considered resistance to seven insecticides of importance to current malaria vector control. In addition, pirimiphos-methyl, α-cypermethrin, bifenthrin, cyfluthrin, and etofenprox are approved and available for use (4). Despite the call by WHO in 2012 for countries to undertake routine resistance monitoring to support operational activities (38), the combined pan-African dataset for 2000–2015 is still too sparse to assess resistance associations comprehensively across the full array of insecticides available for malaria control. Due to the high variability of resistance prevalence across the A. gambiae complex and the measurement error associated with susceptibility bioassays, a comprehensive spatial coverage of samples, as well as replicate sampling, is required to robustly assess variation in resistance phenotypes and any associations with genetic variation.

In conclusion, despite regional variation in species and specific resistance mechanisms, we have demonstrated clear broad-scale associations in the average prevalence of resistance among the pyrethroids and between the pyrethoids and DDT and provided a verification of Vgsc mutations as useful predictors of DDT and pyrethroid resistance. Associations across pyrethroids have the potential to improve our predictions of the prevalence of resistance to insecticides within this class. In resource-limited environments, testing a single pyrethroid insecticide can provide control programs with useful estimates of resistance to other pyrethroids. However, caution is needed because there is uncertainty in these predictions and where resources are sufficient the best approach is still to conduct resistance monitoring for all insecticides being used for malaria control, or being considered for future control, at multiple locations (39).

Materials and Methods

Insecticide Resistance Bioassay Data.

Insecticide resistance bioassay data were obtained from the collated data for bioassays conducted on mosquito samples collected globally described by Coleman et al. (26). The data record the number of mosquitoes in the sample and the proportional sample mortality resulting from the bioassay, as well as information about the mosquitoes tested, the sample collection site, and the bioassay conditions and protocol (26). This information was used to select a subset of records for inclusion in our study (). In summary, we consider only bioassays conducted using standard WHO susceptibility tests on mosquito samples belonging to the A. gambiae complex, conducted over the period 2000–2015 on samples collected within two separate spatial regions selected because they had sufficient data (Fig. 1). The final dataset contained a total of 5,595 data points () but does not provide values for all of the insecticides of interest at each of the 1,183 locations so associations among these insecticides could not be analyzed directly. Further, the spatiotemporal distribution of data for each insecticide is not uniform or random and each distribution incorporates sampling biases (Fig. 1 and ). For these reasons we evaluated the relationships among the insecticides by incorporating these datasets from field populations into a geostatistical model. We found that performing independent model fits to each separate spatial region improved prediction accuracy (see and the next section for further explanation of how this was assessed).

Spatiotemporal Bayesian Statistical Models.

We assess resistance associations across different insecticide types using a Bayesian geostatistical model that extrapolates predictions across space and time by estimating the spatiotemporal correlation structure in the observations. Let denote the proportional mortality record for a bioassay using insecticide type A conducted on a sample collected at geographic coordinates s and sampling time t. To improve posterior validation (assessed using the methods described below) we applied two transformations to these observations, the empirical logit transformation followed by the inverse hyperbolic sine transformation (40). We then assumed that these transformed observations, denoted arise from a spatiotemporal process governed by the measurement equation

where is a constant/intercept, is a Gaussian process modeled as a spatiotemporal Gaussian Markov random field (GMRF) (1, 41, 42), and is Gaussian white noise To jointly model bioassay mortality observations for two insecticides A and B where the observation locations and times differ between the two insecticides, we used a LMC (41, 43, 44). The model of the logit-transformed mortality observations for insecticide B is then given by

where is a constant coefficient and and are defined in the same way as above, replacing A with B in the subscript. We define a Bayesian hierarchical formulation for the joint model given by [] and [], using a vector of prior probability distributions for the hyperparameters where and are the parameters of and (). We implement these models using the R-INLA package (www.r-inla.org) to obtain estimates of the posterior distributions of and (45, 46). To quantify the association between resistance across insecticide types A and B we draw 1,000 simulations from the posterior distribution of the mean mortality for each insecticide type, and across a set of locations j that are defined by the mesh nodes of the GMRFs and (43). For each posterior draw we calculate the Pearson correlation coefficient, r, between the predicted mean mortality values for each insecticide type and estimate the posterior distribution of r based on all posterior draws. We use similar LMC methods to jointly model bioassay mortality observations for up to four different insecticide types A, B, C, and D ( and ref. 44).

We compare performance of these LMCs that allow interactions between resistances across each insecticide types to the performance of models that assume independence of bioassay mortality across different insecticide types to assess whether a joint model best predicts the patterns of resistance seen in field populations. Independence across insecticides is modeled by using Eq. to represent proportional mortality for all insecticides, replacing A by B, C, or D in the subscript. To compare the predictive accuracy of the two model types we compare the WAIC, which use the full posterior distribution of the predicted parameters (27, 28). We also calculated the SE of the difference in the pointwise WAIC values of each model to give an indication of the uncertainty associated with the difference in the WAIC (47) (). For the organophosphates class where we have a small number of resistance observations for each insecticide type, we use a modified version of Eq. that assumes the same pattern of spatiotemporal variation in resistance for both insecticides but allows resistance to each insecticide to differ by a fixed effect. We perform posterior validation using posterior predictive checks, probability integral transform (PIT) histograms (on out-of-sample data), and 10-fold out-of-sample cross-validation ().

Vgsc Allele Prevalence Data.

We obtained data on the prevalence of Vgsc mutations in mosquito samples from the dataset collated by Coleman et al. (26), which describes observed frequencies of point mutations in the Vgsc gene. We extracted those records that reported the proportion of different Vgsc alleles in A. gambiae complex samples collected within the West and East regions of Africa (Fig. 1) over the period 2005–2015. Although some of these data report individual genotypes, we use these to derive the Vgsc allele frequency, because many of the studies reported only allele frequencies. We then applied a process of screening the data to exclude potentially biased estimates of Vgsc allele frequency (). This conservative approach resulted in a reduction in data volume from 1,429 records to 316, incorporating 215 and 101 observations for the West and East regions, respectively (Fig. 1).

Modeling Associations Between Genotype and Phenotype Prevalence Across Populations.

We used OLS linear regression to assess associations between the predicted mean mortality following exposure to a given insecticide type and the prevalence the two Vgsc mutations, L1014F and L1014S, across the genotyped mosquito samples described above. Values of the predicted mean mortality were obtained from our fitted geostatistical models (Figs. 2 and 3) and were represented on the scale of the transformations applied to the insecticide resistance bioassay data described above. Before model fitting we applied the empirical logit transformation to the Vgsc mutation prevalence observations. We verified the normality of the residuals, using the Shapiro–Wilk test (P = 0.154). Examination of the residuals revealed contemporaneous spatial correlation, so we performed inference on the regression coefficients using cluster-robust SEs (48, 49), setting the cluster radius to 100 km (, section S5).