Insecticide resistance and behavioural adaptation as a response to long-lasting insecticidal net deployment in malaria vectors in the Cascades region of Burkina Faso

The decline in malaria across Africa has been largely attributed to vector control using long-lasting insecticidal nets (LLINs). However, this intervention has prompted widespread insecticide resistance (IR) and been associated with changes in mosquito behaviour that reduce their contact with LLINs. The relative importance and rate at which IR and behavioural adaptations emerge are poorly understood. We conducted surveillance of mosquito behaviour and IR at 12 sites in Burkina Faso to assess the magnitude and temporal dynamics of insecticide, biting and resting behaviours in vectors in the 2-year period following mass LLIN distribution. Insecticide resistance was present in all vector populations and increased rapidly over the study period. In contrast, no longitudinal shifts in LLIN-avoidance behaviours (earlier or outdoor biting and resting) were detected. There was a moderate but statistically significant shift in vector species composition from Anopheles coluzzii to Anopheles gambiae which coincided with a reduction in the proportion of bites preventable by LLINs; possibly driven by between-species variation in behaviour. These findings indicate that adaptations based on insecticide resistance arise and intensify more rapidly than behavioural shifts within mosquito vectors. However, longitudinal shifts in mosquito vector species composition were evident within 2 years following a mass LLIN distribution. This ecological shift was characterized by a significant increase in the exophagic species (An. gambiae) and coincided with a predicted decline in the degree of protection expected from LLINs. Although human exposure fell through the study period due to reducing vector densities and infection rates, such ecological shifts in vector species along with insecticide resistance were likely to have eroded the efficacy of LLINs. While both adaptations impact malaria control, the rapid increase of the former indicates this strategy develops more quickly in response to selection from LLINS. However, interventions targeting both resistance strategies will be needed.

. The predicted reduction in postexposure mortality over time as observed with the DD was also detected at 15X the DD (df = 1, χ 2 = 11.25, p = 0.001; Supplementary Table 1). Long-term increases in IR were evident after controlling for statistically significant variation between villages (all concentrations) and season (only detected in bioassays using the DD; Supplementary Table 1).
The mean sporozoite rate (SR) in An. gambiae s.l. across the study area was 3.48% [95% CI 1.51-5.26%]. After controlling for seasonal and spatial variations, their was evidence of a decline in mean SR ~ 5% to ~ 2% over the study period (z = − 2.5, p = 0.01, Supplementary Table 4, Fig. 5). There was no significant difference in the SR betweenAn. coluzzii and An. gambiae; or between An. gambiae s.l. captured host-seeking indoor versus outdoor (Supplementary Table 4). Annual Entomological Inoculation Rates (EIR) were calculated for the subset of 6  . This decline was in accordance with the observed longitudinal declines in the human biting rate and SR across the study. Despite this decline, residents were still predicted to be exposed to ~29 infective bites per person per year after adjustment for the proportion of bites (~ 85%) that can be prevented through effective use of LLINs.

Discussion
Here we tested for evidence of systematic changes in insecticide and behavioural resistance traits in malaria vectors over 2.5 years following a mass LLIN distribution in Burkina Faso. Consistent with expectation, there was a high prevalence of insecticide resistance (IR) within these malaria vector populations, which further intensified over the period of the study, with post-exposure mortality falling by ~ 50% (from 38% at the start to 17% at the end). In contrast, although baseline levels of "behavioural avoidance" traits were higher than expected (e.g. > 50% outdoor biting and resting), behavioural phenotypes were relatively stable across the study period. There was however a reduction in An. gambiae s.l. abundance (Human Biting Rate) and a shift in vector species Table 3. Total numbers of blood-fed female An. gambiae s.l. caught using Resting Bucket Traps (RBT) in the 12 villages, from October 2016 to December 2018. Results are displayed by species, trapping location (Indoor or outdoors) and host blood source, pooled over village. 'HBI' indicate the 'Human Blood Index'; estimated as the proportion of blood meals taken from humans out of the blood meals whose source could be identified. HBI indicate the proportion of mosquito that blood-fed on human.  www.nature.com/scientificreports/ composition across the study period; with a moderate reduction in the proportion of An. coluzzii relative to An. gambiae. The potential epidemiological impacts of these ecological changes are complex. On one hand, there was a significant drop in the Entomological Inoculation Rate between years, indicating a fall in human exposure and malaria transmission. On the other, the proportion of exposure to An. gambiae s.l. expected to be preventable by LLIN use also fell over the study period (π i , ~ 7% fall); indicating a moderate erosion of protection from residual transmission through time. We hypothesize that the fall in π i was due to a shift in vector species composition, with the slightly more exophagic An. gambiae increasing relative to An. coluzzii through time. Thus, we hypothesize that while both IR and avoidance behaviours are contributing to the limited impact of LLINs in this setting, behavioural traits appeared to be responding more slowly to selection from LLINs. The rise in IR across this study is in line with similar increases observed over short time periods in other African settings 49,50 . This acceleration in IR appears to be a relatively recent phenomenon. For example, postexposure mortality rates of An. gambiae s.l. to the DD of synthetic pyrethroids when first introduced in the 1970s (~ 95% 51 ) was relatively similar to that in the 2010s 52 . The intensification of IR over the last decade is likely due to the initiation of mass LLIN distribution programmes in 2010. For example, post-exposure mortality of An. gambiae s.l. to the DD in one of our study sites (Tiefora) was ~ 39% 53 in 2014 compared to only ~ 15% reported here (2016-18). Additionally, in Tengrela (another study village), post exposure mortality to DD declined from ~ 92 to ~ 19% between 2011 and 2013 37 . In light of the rapid rise in IR, the government of Burkina Faso switched to new net products that combine pyrethroids with other insecticide classes in the most recent (2019) mass ITN distribution campaign (e.g. INTERCEPTOR G2: alpha-cypermethrin and chlorfenapyr; PBO-LLINs: deltamethrin and the synergist piperonyl butoxide).
While insecticide resistance increased significantly over the study, mosquito vector behaviours were generally static. However, some mosquito behavioural traits observed here differ substantially from historical reports in the study area; raising the possibility of more gradual long-term changes that may not be detectable over a few years. For example, previous studies from the Central, Plateau Central and West regions (2001-2015) of Burkina Faso reported An. gambiae s.l. was more likely to bite indoors than outdoors 43,54 , or had an even split between indoor and outdoor biting 46 . In contrast, here An. gambiae s.l. was more likely to host seek outdoors [~ 54%; 95% CI ~ 51-57%]. Additionally, a recent review of An. gambiae s.l. biting behaviour from a range of African countries concluded that in general > 80% of vector bites occur indoors 28 ; substantially higher than observed here. Anopheles gambiae s.l. were also found in indoor and outdoor resting traps with relatively similar frequency (53% versus 47%); in contrast to earlier studies in south western Burkina Faso where indoor resting was more common (~ 67% 55 ). The Human Blood Index (~ 65%) of these An. gambiae s.l. populations was also lower than previously reported in Burkina Faso (> 77% 56 ) and Benin (> 90% 57 ). Thus, although there was no clear evidence of mosquito behavioural change across the two years of this study, comparison with historical and wider regional data suggest that slower acting behavioural adaptations may be occuring in tandem with more rapid rises in insecticide resistance. www.nature.com/scientificreports/ LLINs have also been hypothesized to generate selection on mosquito biting times; with vectors shifting to bite earlier in the evening before people go to bed 58 to avoid this intervention. Here, An. gambiae s.l. still exhibited the characteristic "late night" pattern of biting, with no evidence of a shift over the study. The biting pattern observed here is consistent with early work on An. gambiae s.l. before mass LLIN use which also found biting rates peaked around 00 h 12,59 . Similarly, a study in western Kenya found An. gambiae s.l. and An. funestus continued biting late in the night even in the presence of LLINs 60 . However, other studies have detected substantial shifts in malaria vector biting times in association with interventions, including a recent report from Senegal of An. funestus biting earlier in the morning and even into daylight hours 61 . Adaptations in biting time may arise more rapidly in other settings due to variation in local ecology and background levels of IR.
Some previous studies that have reported shifts in malaria vector behaviour (e.g. 27,62 ) did not identify An. gambiae s.l. to species level, thus preventing interpretation of whether changes were due to within-species adaptations or ecological shifts in species composition (e.g. reduction in the more endophagic An. gambiae relative to more exophagic An. arabiensis). While behavioural phenotypes remained relatively fixed within vector species in this study, there was evidence of a longitudinal shift in malaria vector species composition characterized by a reduction in An. coluzzii relative to An. gambiae. The proportion of outdoor biting was slightly higher in An. gambiae (55%) than An. coluzzii (51%), thus the relative decline in An. coluzzii is consistent with the prediction that LLINs have a greater impact on endophagic vector species. For example, highly anthropophagic and endophagic vector species like An. gambiae 63 and An. quadriannulatus 64 declined more substantially than the more zoophagic and exophagic An. arabiensis following ITN introduction in East Africa. This adds to the growing body of evidence from across Africa showing that LLINs can provoke shifts in malaria vector species communities. Such ecological shifts may occur more rapidly than longer-term behavioural adaptations within species.
Taken in combination, these results suggest that IR may be the first line of defence to LLINs compared to changes in mosquito behaviour, and/or that behavioural phenotypes have less capacity for adaptation than the traits underpinning IR. Differences in the rate of change in IR and mosquito behavioural traits also likely reflect their genetic basis. Mosquito behaviours are likely complex multigenic traits 65 , governed by many genes (some possibly antagonistic) that may limit capacity to respond quickly to selection. There is some evidence that behaviours like human host choice 65 and outdoor/indoor biting or resting 66 have a genetic basis in African malaria vectors. The genetic basis of IR traits is well established 18,67 , with some IR mechanisms linked to a single point mutation (e.g. Knock-down resistance 16 ; indicating this type of resistance can respond rapidly and efficiently to selection. Vector populations in Burkina Faso have a high frequency of target site mutations (e.g. frequency of the L1014F type > 0.8 68 ) which could have provided the foundation for quicker adaptation to insecticides in contrast to slower development of behavioural resistance.
The potential epidemiological impacts of the vector adaptations observed here are complex. On one hand, results indicate that use of effective LLINs during typical sleeping hours (10 pm-5 am) can still effectively prevent the bulk of human exposure to malaria vectors (~ 85%) in this setting. Additionally, vector abundance, sporozoite infection rates and the associated Entomological Inoculation Rate fell between the two years of this study; suggesting a transmission decline. However, this degree of protection expected from LLIN use (85%) is somewhat lower than estimated in other parts of Burkina Faso (90%, in 2002-2004), and Africa (95-99% 69 ). Furthermore, the longitudinal changes in vector populations described here could further diminish LLIN impact in several ways. First, the intensification of IR is likely to erode the community protection afforded to non-net users. Second the proportion of exposure preventable by using LLINs (π i ) was predicted to decline over the study period. Even if the proportion of human exposure preventable by LLIN use remained at current levels, high levels of residual transmission will persist because of the relatively high baseline abundance and malaria infection rates in these vector populations. People in this area are estimated to be exposed to between ~ 29 infective bites per person per year even if they sleep under an LLIN between 10 pm-5 am. The small shifts in P fƖ and π i observed through time could further increase this exposure. Interventions that tackle both insecticide resistant and outdoor biting mosquitoes will thus be needed to tackle residual transmission.

Conclusions
In a longitudinal study in Burkina Faso, we show that IR increases more rapidly than behavioural adaptations following the mass distribution of insecticidal nets. Although mosquito behaviours stayed relatively stable over the study, baseline traits indicated a higher capacity for 'behavioural avoidance' of nets (through outdoor biting and resting) than historical data. This highlights the possibility that behavioural traits are adapting but at a more gradual rate than physiological resistance. As most human exposure to infected mosquitoes still occurs indoors in this and other African settings, effective indoor interventions should still be prioritized including the use of novel insecticide classes on nets, sprays and non-insecticidal approaches. However, the proportion of exposure occurring outdoors is notable and increasing. Supplementary measures tackling human exposure in and outdoors will thus be required to make further progress and tackle residual transmission in this and other high burden African settings.

Methods
Study site. This study was conducted in 12 villages within the Banfora District, in the Cascades Region, south-western Burkina Faso (Supplementary Table 5, Fig. 6). In 2018, the population size of the region was 822,445 70 . The region has a humid savannah climate characterized by a rainy (May to October) and dry season (November to April). Annual rainfall ranged from 600 to 900 mm; with a mean temperature of ~ Tables 6  and 7). Larvae were collected from aquatic habitats and reared to adulthood under standard conditions 73 . Insecticide resistance assays were performed by exposing groups of 22-27 adult female An. gambiae s.l. to the World Health Organization discriminating dose for deltamethrin (DD = 0.05%, a concentration that can kill 99.9% of a susceptible population 74 ), 5 times (0.25%), 10 times (0.5%) and 15 times (0.75%) the DD respectively (Supplementary Table 6). Mosquitoes were exposed to insecticide-treated papers obtained from the WHO/Vector Control Research Unit, University Sains Malaysia 75 , following WHO guidelines for Tube bioassays 74 .

Assessment of human biting rate and behaviours.
Vector surveillance was conducted to measure human biting rates (HBR) and mosquito resting behaviours bi-monthly at all 12 sites for 16 months (Oct 2016-Feb 2018, Fig. 6). Additionally, surveillance continued for a further 10 months (Feb-Dec 2018) in a subset of 6 villages to generate longer-term data. The longer-term study villages were selected to achieve a relatively broad spatial distribution. Host-seeking mosquitoes were sampled twice a month in each village from two different households per day using Human Landing Catches (HLC). Here, men collected mosquitoes landing on their exposed leg using a mouth aspirator and flash torches 76 . Collection took place both inside houses and outdoors (the peridomestic area) from 7 pm to 6 am each night. During each hour of collection, participants caught mosquitoes for 45 min, followed by a 15-min rest break. Participants involved in mosquito collections rotated between indoor and outdoor trapping stations each hour to avoid confounding location with individual differences in attractiveness to mosquitoes.
In addition, Resting Box Traps (RBTs) were placed in and outside houses 77 to sample resting mosquitoes. These RBTs were made locally using 20 L plastic buckets, with their inner surface covered with moistened black cotton cloth to create a high contrast and humid environment. On each night of collections, two RBTs were placed in the same households where HLCs took place (set within different houses in the compound). Inside houses, RBTs were placed on the floor in a relatively shaded corner of the sitting room. Two RBTs were also set www.nature.com/scientificreports/ outdoors at ~ 8 m from houses to capture outdoor resting mosquitoes. RBTs were set up at approximately 7 pm and emptied the following morning (~ 5 am) using electrical aspirators.

Mosquito processing. Mosquitoes collected in HLCs and RBTs were transferred to an insectary in Banfora
town and sorted to species complex level using morphological keys 78 . Of those identified as belonging to the malaria vector group Anopheles gambiae s.l., a subset of 7852 females (~ 20% of total) were selected to provide a representative sample from each month, village, trapping location (indoor vs outdoor) for identification to species level by PCR 79 . The subsampling strategy is described elsewhere 80 . Furthermore, malaria infection in mosquitoes was assessed by testing for the presence of Plasmodium falciparum circumsporozoite protein (CSP) in their head and thoraxes using a monoclonal sandwich Enzyme Linkage Immuno-Sorbent Assay (ELISA) developed by 81 . Additionally, female An. gambiae s.l. from the RBT collections were visually graded according to their repletion status (abdominal condition) into categories of blood-fed, unfed, gravid, and half gravid 82 . Blood-meal identification was carried out on blood-fed An. gambiae s.l. to determine the origin of their bloodmeal (human, cattle or both 83 ) using a direct Enzyme Linkage Immuno-Sorbent Assay (ELISA 84 ). The degree of "anthrophagy" in vectors was estimated in terms of the "human blood index" (HBI) 85 , defined as the proportion of identified blood meals taken from humans. In ELISA tests (for sporozoite and blood meal source identification), two technical replicates of each sample were run in two different microplates at the same time and retested in cases where the first result was ambiguous. The absorbance of the solutions/reactions at the end of each ELISA was measured using microplate reader (ELX808; BIO-TEK) at 450 nm. To avoid any false positives (due to background noise), a sample was considered positive for an assay when its optical density (OD) was twofold higher than the average of the OD of both negative controls. Positive controls were also used in all ELISAs to ensure the procedure was working.  Table 8) to assess longer-term trends and seasonality occurring over the 26 months of the study. For the longer-term trend, a temporal variable was created by assigning a continuous value to each day of collections from the first (1 October 2016) until the last (4 December 2018), hereafter called "Longer_term" (Supplementary Table 8). This was incorporated into models as a continuous covariate to test for evidence of a consistent temporal rise or decline after accounting for seasonal and spatial variation. Seasonality was incorporated by fitting a smoothing spline where each day of the year was classified on a scale running from 1 (set as January 1st) to 365 (December 31st) hereafter called "Season" (Supplementary Table 8). Spatial variation was also assessed by fitting village as a fixed effect. The mean nightly temperature and humidity at collection households were derived from hourly values (7 pm-6 am) obtained from data loggers and included as additional explanatory variables (Supplementary Table 8).

Data analysis.
However, due to the low number of mosquitoes caught in RBTs, it was not possible to build a model for testing for spatial and long-term variation in the abundance of resting mosquitoes. Mosquito host choice was assessed in terms of the Human Blood Index (HBI) defined as the proportion of blood-fed An. gambiae s.l. that tested positive for human blood out of the total from which blood-meals were identified (N = 94). We also calculated the Entomological Inoculation Rate (EIR); defined as the average number of infective bites a person (ibppn) would expect to receive from An. gambiae s.l. in a given location per year) for the six villages that were monitored over 2 years. Here EIR was calculated as the product of the mean nightly human biting rate (ma) and vector sporozoite infection rates (SR) multiplied by 365 days 41 monitored over two years.
Mosquito mortality rate after exposure to insecticide. Insecticide resistance was measured in terms of mosquito mortality 24 h after exposure to deltamethrin (different doses). According to World Health Organisation (WHO) guidelines 74 , for a test to be validated the mortality should be less than 5% in the control groups. If, mortality in the control group ranges between 5 and 20%, the guidelines recommend using the Abbott formula 74 to correct for anticipated background mortality. If the mortality rate in the control group exceeds 20%, results from the test should be discarded. Here, average mortality rates recorded in the control group on the same day as insecticide exposures were performed were all below 5%, and could thus be retained for use without correction. We constructed separate models of mosquito mortality for each of the insecticide concentrations used (Models 1 to 4, Supplementary Table 8). As the response variable is all these analyses was binary (e.g '0' is alive, '1' is dead); all models used a binomial distribution. Each model tested for seasonal and long-term variation in mosquito mortality following insecticide exposure (Supplementary Table 8), and variation between villages.
Anopheles gambiae s.l. species composition, biting and resting behaviours. Models were constructed to test for longer-term (systematic increase or decrease across the 26 month study period) and seasonal variation in species composition (Model 5), HBR (Model 6) and outdoor biting (Model 7 and 8) and resting (Model 12) as described in Supplementary Table 8. Additionally, long-term and seasonal variation in the mean mosquito biting times (both An gambiae s.l. overall, and by species; Model 10-11). Model 6 (Supplementary Table 8 www.nature.com/scientificreports/ HBR included an interaction term between village and year in addition to other main effects. Model for mosquito median biting time includes an interaction between village and location. Vector species composition was defined as the proportion of An. coluzzii within An. gambiae complex. This was defined as binary proportion as only two major vector species (An. coluzzii and An. gambiae) account for almost all within the An. gambiae complex in this region. The proportion of outdoor biting was defined as the proportion of An. gambiae s.l. collected in outdoor HLCs out of the total in outdoor and indoor HLCs. As described in Supplementary Table 8, different distributions were incorporated into models depending on the nature of the response variable (e.g count data fit to Poisson or Negative Binomial model; proportion data fit to a binomial distribution).
Human exposure to Anopheles gambiae s.l. bites. Data on the time and location of biting were combined for estimating two metrics of human exposure to bites from An. gambiae s.l.: the proportion of An. gambiae s.l. caught during hours when typically, ≥ 50% of people were indoors and likely to be in bed (P fƖ ). Based on visual surveys described in 80 over 672 households, we estimated that > 50% of residents were indoors between 10 pm and 5 am consistently across the study period; although more recent, detailed anthropological investigation indicates time indoors may vary substantially between individuals and seasons and be lower assumed here 48 . The P fl (Model 14 in Supplementary Table 8) was estimated following 42 by dividing the total number of s.l. collected indoors and outdoors when ≥ 50% of people were indoors (10 pm to 5 am) by the total number collected that night. The proportion of human exposure to bites that occur indoors (π i ) was then calculated by dividing the number of An. gambiae s.l. collected indoors between 10 pm and 5 am by itself plus those collected outdoors when > 50% of residents were outdoors and thus unprotected by LLINs (between 7-10 pm and 5-6 am (Model 15, Supplementary Table 8). As π i and P fl are proportions, they were modelled following a binomial distribution. Within the R statistical software 86 GAMMs within the 'mgcv package' 87 augmented with the lme4 package 88 known as GAMM4 were used to test for associations between IR, all vector ecological and behavioural metrics and explanatory variables of village, season, "longer-term" trend and environmental factors. Details of fixed and random effects and distribution are given in Supplementary Table 8. Prior to starting the research, the project aims, and objectives were explained to community leaders in each village. Signed informed consent was also obtained from all household owners where mosquitoes were collected, and volunteers who took part in mosquito collections by HLC. We confirm that all methods were performed in accordance with the relevant guidelines and regulations.