Emergence of knock-down resistance in the Anopheles gambiae complex in the Upper River Region, The Gambia, and its relationship with malaria infection in children

Background Insecticide resistance threatens malaria control in sub-Saharan Africa. Knockdown resistance to pyrethroids and organochlorines in Anopheles gambiae sensu lato (s.l.) is commonly caused by mutations in the gene encoding a voltage-gated sodium channel which is the target site for the insecticide. The study aimed to examine risk factors for knockdown resistance in An. gambiae s.l. and its relationship with malaria infection in children in rural Gambia. Point mutations at the Vgsc-1014 locus, were measured in An. gambiae s.l. during a 2-year trial. Cross-sectional surveys were conducted at the end of the transmission season to measure malaria infection in children aged 6 months–14 years. Results Whilst few Anopheles arabiensis and Anopheles coluzzii had Vgsc-1014 mutations, the proportion of An. gambiae sensu stricto (s.s.) mosquitoes homozygous for the Vgsc-1014F mutation increased from 64.8 to 90.9% during the study. The Vgsc-1014S or 1014F mutation was 80% higher in 2011 compared to 2010, and 27% higher in the villages with indoor residual spraying compared to those without. An increase in the proportion of An. gambiae s.l. mosquitoes with homozygous Vgsc-1014F mutations and an increase in the proportion of An. gambiae s.s. in a cluster were each associated with increased childhood malaria infection. Homozygous Vgsc-1014F mutations were, however, most common in An. gambiae s.s. and almost reached saturation during the study meaning that the two variables were colinear. Conclusions As a result of colinearity between homozygous Vgsc-1014F mutations and An. gambiae s.s., it was not possible to determine whether insecticide resistance or species composition increased the risk of childhood malaria infection. Electronic supplementary material The online version of this article (10.1186/s12936-018-2348-8) contains supplementary material, which is available to authorized users.


Background
Between 2000 and 2015, the prevalence of Plasmodium falciparum infection in sub-Saharan Africa (SSA) has halved due to the mass deployment of long-lasting insecticidal nets (LLINs), and to a lesser extent, indoor residual spraying (IRS) [1]. This has, however, increased selection pressure for insecticide resistance in malaria vectors, particularly resistance to pyrethroids, the only insecticide class used currently for treating bed nets. The strength and distribution of insecticide resistance has increased over time and there is growing concern that this will lead to control failure, which has the potential to reverse many of the gains seen in malaria control [2]. One of several mechanisms through which mosquitoes become resistant to insecticides is through mutations in the insecticide target site. Three different point mutations in the voltage-gated sodium channel gene confer knockdown resistance (kdr) to pyrethroids and organochlorines such as dichlorodiphenyltrichloroethane (DDT) in Anopheles gambiae s.l. [3][4][5].
There has been little longitudinal insecticide resistance monitoring in The Gambia, but the general impression is that levels of insecticide resistance are low, but rising. Shortly after the introduction of permethrin-treated bed nets in The Gambia in the early 1990s, there was little or no resistance to DDT or permethrin [6,7]. Prior to the nationwide DDT IRS campaign in 2009, DDT resistance was found in one site bordering Senegal, but mosquitoes from the site in the Upper River Region (URR) were fully susceptible to permethrin, deltamethrin and DDT [8]. Tests performed during a cluster-randomised controlled trial in the URR which compared the efficacy of LLINs versus LLINs and IRS with DDT against malaria in children found complete susceptibility of An. gambiae s.l. to DDT and permethrin in 2010 and some loss of susceptibility the following year [9]. Another study, conducted in the same year but outside the study area indicates that there may be pockets of high resistance in the URR [10]. Research also suggests that insecticide resistance may be partly responsible for the heterogeneities in malaria transmission across the country [11]. Although malaria has been declining in The Gambia since 2000 [12,13], there is continued moderately high seasonal transmission in the URR despite high vector control coverage [14].
The study aimed to examine risk factors for kdr resistance in An. gambiae s.l. and its relationship with malaria infection in children in The Gambia.

Study site
The study was conducted in the URR (regional capital: Basse Santa Su, 13.3167°N, − 14.2167°W), a rural area of open Sudanian savannah which is divided into north and south banks by the River Gambia. Malaria transmission is highly seasonal being associated with the annual rains which occur from June to October. LLIN use by study children was 55% at baseline and IRS with DDT was implemented in 2009, the year prior to study start.

Data collection
This secondary analysis uses data from a cluster-randomised controlled trial assessing the efficacy of LLINs versus LLINs and IRS with DDT against malaria among children aged 6 months-14 years in the URR of The Gambia [9]. The study design and results are described in full elsewhere [9,15]. In brief, 70 clusters of villages were randomly allocated to receive either LLINs or LLINs plus IRS. Permethrin-treated LLINs (2%; Olyset Nets, Sumitomo Chemicals, Japan) were distributed in both arms at the start of the 2010 transmission season to achieve high coverage of sleeping places. IRS with DDT (2 g/m 2 , DDT 75% wettable powder; Hindustan Insecticides, New Delhi, India) was applied to dwelling rooms at the start of each transmission season in the IRS-LLIN arm. Surveys conducted at the end of the transmission season in 2010 and 2011 measured the prevalence of P. falciparum infection in a cohort of children.
Entomological data were collected in 32 clusters (16 in each arm) in six sentinel rooms per cluster. Clusters were chosen purposively for logistical reasons (Fig. 1). Sampling was performed monthly from June to the end of December in 2010 and 2011, and every 2 months during the intervening dry season. Mosquitoes were collected overnight from sentinel rooms in which an adult slept under an LLIN using a CDC light trap. The epidemiological dataset was restricted to children who resided in the entomological clusters (1543 children in 2010 and 1564 children in 2011).

Mapping and spatial analysis
Digitised maps produced by the Japan International Cooperation Agency under The Japanese Government Technical Cooperation Programme and The Government of the Republic of The Gambia from 2002 were obtained. Global Moran's spatial autocorrelation coefficient, I, was calculated at 1 km intervals between 9 km (the shortest distance at which all sampling point locations had at least one neighbour) and 25 km to examine spatial independence in species distributions. The z-score returned indicated the intensity of clustering. Mapping and spatial analysis was performed using ArcGIS ® software (Release 10.4.1, Environmental Systems Research Institute: Redlands, CA).

Statistical analysis
Transmission seasons were defined as 16 August-31 December 2010 and 15 August 2011-1 January 2012 to avoid the months prior to and during application of IRS and the intervening dry season. Proportions of mosquitoes by species and kdr status over time and by village were calculated. Mixed effect logistic regression models including cluster as a random effect were used to determine the relationship between vector species of individual mosquitoes and Euclidean distance of the cluster from the River Gambia, and secondly, the effect of variables such as year and study arm on kdr status of individual mosquitoes, whilst controlling for species. Variables were tested for departure from linear trend where necessary. Stepwise selection procedures and likelihood ratio tests were used to determine the combination of covariates, which fitted the data best. Mixed effect logistic regression models were also used to look at the effect of i) cluster-level kdr status (prevalence of any Vgsc-1014 mutation i.e. any mutation at the Vgsc-1014 locus, and homozygous Vgsc-1014F mutation by cluster) and ii) prevalence of An. gambiae s.s. by cluster on the odds of P. falciparum infection in individual children at the end of transmission season surveys, adjusting for clustering and confounding variables. These three explanatory variables were fit as linear variables and expressed as the odds ratio for the effect of a 1 and 10% increase in these variables on the prevalence of malaria infection. R 2 and the variance inflation factor (VIF) were calculated to identify colinearity between variables. Goodness of fit of models evaluating the effect of cluster level prevalence of either An. gambiae s.s. or homozygous Vgsc-1014F mutations on malaria infection were compared using the Akaike information criterion (AIC). Models were also run to evaluate the effect of absolute numbers of mosquitoes per cluster with any Vgsc-1014 mutation and the homozygous Vgsc-1014F mutation, and absolute number of An. gambiae s.s. on the odds of P. falciparum infection in individual children. Statistical analyses were performed using Stata 14 (College Station, TX, USA).

Results
A total of 6853 An. gambiae s.l. were caught in the 32 sampling sites over the two transmission seasons. Of these, 6828 (99.6%) were identified to species: 71.3% were An. arabiensis, 15.0% An. gambiae s.s., 12.3% An. coluzzii, and 0.1% hybrid (Fig. 2). Higher numbers were caught during 2010 when there was unusually high rainfall and extensive flooding, compared to 2011 when flooding was limited to areas beside the river. During the 2010 transmission season, 76.1% of An. gambiae s.l. were An. arabiensis, 12.0% An. gambiae s.s., 10.1% An. coluzzii and 0.2% hybrid Twenty-nine mosquitoes were caught during the dry season (of these 25 were An. arabiensis, 2 An. gambiae s.s. and 2 An. coluzzii). Spatial autocorrelation was found in species distributions with peak autocorrelation operating between 9 and 14 km depending on the species and year. Analysis of species distributions over the two transmission seasons, showed that An. gambiae s.s. was more common further away from the river (Odds ratio, OR for every km away from the river = 1.29, 95% CI 1.21-1.38, p < 0.001) (Figs. 3, 4). Conversely, both An. arabiensis and An. coluzzii were more common closer to the river (An. arabiensis OR = 0.88, 95% CI 0.83-0.94, p < 0.001; An. coluzzii OR = 0.91, 95% CI 0.85-0.98, p = 0.01). Similar patterns were found when each year was analysed separately.
In Vgsc-1014 mutations were found in all species sampled but at differing levels. An. arabiensis were predominantly wild-type (73.1% during 2010 and 58.1% during 2011), although the proportion with heterozygous Vgsc-1014S mutations increased from 20.5% in 2010 to 28  which mirrors the increase in the proportion of An. gambiae s.s. in these areas (Fig. 3, 4).
In a multivariable model, species, study arm and year of survey were associated with odds of any Vgsc-1014 mutation (Table 3). Adjusting for year and study arm, An.

Discussion
These findings illustrate the temporal and spatial pattern of the An. gambiae complex and Vgsc-1014 mutations in the URR of The Gambia and their association with malaria infection in children from 2010 to 2011. To the authors knowledge, this is the first study to adopt a landscape approach with intensive entomological sampling to understand factors related to the distribution of kdr in malaria vectors.
As in previous work in the URR [21][22][23][24], An. arabiensis was the most abundant member of the An. gambiae complex and persisted longer into the dry season than the other species. The frequency of hybrids was 0.1-0.2%, a slightly lower proportion than that shown by others in The Gambia [24,25]. Anopheles gambiae s.s. were more common in villages away from the River Gambia which corresponds with previous studies that reported An. gambiae s.s. prefers small rain-dependent larval habitats on free-draining soil covered with open woodland savannah or farmland [26,27]. Anopheles arabiensis was more common near the river suggesting that their aquatic habitats are common in wetlands, such as rain-fed ricefields, adjacent to the river [28]. Indeed, previous studies found An. arabiensis in water bodies along the edge of the alluvial soils, particularly in areas of rice cultivation in the Central River Region [29]. Anopheles coluzzii was also more common closer to the river which supports research showing this species exploits semi-permanent aquatic habitats that are also frequented by An. arabiensis [26,30,31]. Although the literature suggests the distribution of species is likely to be due to differential larval habitats, this result may be because An. gambiae s.s. has less flexible host choice behaviours than An. arabiensis [32][33][34] and so contributes a larger proportion of the mosquito catch further away from the river.
Interestingly, fewer An. arabiensis were caught in villages in the double intervention compared to the single intervention arm, while the opposite pattern was seen with An. gambiae s.s. If the Vgsc-1014F mutation translates into phenotypic resistance, this may have given An. gambiae s.s. a competitive advantage over An. arabiensis An. gambiae s.s. since the latter is more endophagic and endophillic and so is thought to be preferentially killed by LLINs [35][36][37].
Species, year of survey and study arm were associated with odds of any Vgsc-1014 mutation. There was a 1.27 increase in the odds of any Vgsc-1014 mutation in the double intervention compared to the single intervention arm and a 1.80 increase in the odds of any Vgsc-1014 mutation between 2010 and 2011. The odds of any Vgsc-1014 mutation was 18.49 times higher in An. gambiae s.s. compared to An. arabiensis. Taken together this suggests that (i) LLINs and DDT used together provide more selection pressure than LLINs alone and (ii) there was an increase in selection pressure over the two years most likely due to the second IRS round, and (iii) selection pressure favours An. gambiae s.s. since it has a higher frequency of kdr. Several studies have shown an increase in the frequency of kdr mutations following implementation of vector control interventions [38][39][40][41]. High levels of kdr in An. gambiae s.s. compared to An. arabiensis may be explained by a greater propensity for indoor resting and feeding of An. gambiae s.s. and, therefore, potential for increased contact with insecticides on walls or LLINs [17,34,42]. The study used DDT for IRS and pyrethroidtreated LLINs, and as such is it unsurprising that selection pressure for development of Vgsc-1014 mutations was high in the double intervention arm. The IRS in the study was performed by government teams using the insecticide selected by the National Malaria Control Programme but an alternative insecticide class, such as an organophosphate or carbamate would have been a better option to reduce selection pressure [10]. In fact, based on insecticide resistance monitoring, the control programme recently started to implement rotation of IRS insecticides beginning with bendiocarb in 2015 and 2016 and pirimiphos-methyl in 2017.
Older children had a higher odds of P. falciparum infection in both of the end of season surveys. It is unclear why females were at increased risk of infection in the first survey but at lower risk in the second survey and this may be an anomalous result. In 2011, children sleeping under an LLIN had half the odds of being infected compared to children not sleeping under a net, but no significant difference was observed in 2010. This result is probably due to chance because of low numbers of children not using an LLIN in this study. There was no significant association between the cluster level proportion of mosquitoes with any Vgsc-1014F mutations and malaria infection in children. However, univariable analysis did show an association between P. falciparum infection in children and the cluster level proportion of An. gambiae s.s. and homozygous Vgsc-1014F mutations specifically, especially in the second year of the study. However, due to colinearity between An. gambiae s.s. and homozygous Vgsc-1014F mutations these variables could not be combined in the multivariable model. Calculation of the AIC was not able to distinguish between models including these two variables and so it was not possible to say whether high proportions of An. gambiae s.s. or homozygous Vgsc-1014F mutations increased P. falciparum infection in children.
This study has several other limitations. Firstly, the analysis used secondary data which meant that the original study was not primarily designed to measure the spatial epidemiology of kdr resistance. Secondly, it did not verify the phenotype of the mosquitoes in bioassays or investigate other resistance markers, such as those involved in metabolic resistance. Insecticide resistance is typically driven by complex interactions between multiple alleles and this dataset only looks at a few alleles. Thirdly, spatial autocorrelation was present in species distributions and therefore the assumption that clusters were independent was false. This may have inflated the value of test statistics and increased the chance of a type I error.
Increased malaria infection in the study children may be explained by differences in the species distribution in the villages, specifically possible higher efficiency transmission by An. gambiae s.s. Anopheles gambiae s.s. is a more efficient vector than An. arabiensis, although it is less clear whether there is a difference in transmission efficiency between An. gambiae s.s. and An. coluzzii [43][44][45]. Alternatively, heterogeneity in malaria infection could also be due to the impact of kdr. Indeed, previous studies have highlighted a lack of decline in malaria in the URR [14,46] and a study of paired high and low malaria prevalence villages in The Gambia suggested that heterogeneous transmission may be partly due to insecticide resistance [11]. Opondo et al. showed that DDT mortality for An. gambiae s.s. was significantly lower in high prevalence compared to low prevalence villages and that there was a significant association between the Vgsc-1014F mutation in An. gambiae s.s. and resistance to DDT and deltamethrin [11]. This mutation was a strong predictor of insecticide resistance and effectively masked the effect of other mutations in this study such as those associated with metabolic resistance. However, the role of kdr is not clear cut [47,48] and several studies show that pyrethroid LLINs were still able to kill An. gambiae despite high kdr frequencies [49][50][51][52]. Analysis did not show any significant relationship between childhood malaria infection and the absolute number of An. gambiae s.s. mosquitoes, mosquitoes with any Vgsc-1014 mutation or mosquitoes with the homozygous Vgsc-1014F mutation per cluster. This is most likely because of low vector numbers in some of the village clusters.

Conclusions
In conclusion, the homozygous Vgsc-1014F mutation occurred predominantly in An. gambiae s.s. and increased almost to saturation during the course of the study. It also occurred at higher frequencies where IRS was used in addition to LLINs, probably because the kdr mutation confers a selective advantage in the presence of insecticides. There was a 13% increase in the odds of malaria infection in children associated with a 10% increase in the proportion of An. gambiae s.l. carrying the Vgsc-1014F mutation in 2011. Moreover there was a 27% increase in the odds of malaria infection with a 10% increase in the proportion of An. gambiae s.s. mosquitoes in 2010 and a 12% increase in 2011. It was, however, impossible to determine whether resistance or species increased the odds of childhood malaria infection since the homozygous Vgsc-1014F mutation was colinear with An. gambiae s.s.