An increasing role of pyrethroid-resistant Anopheles funestus in malaria transmission in the Lake Zone, Tanzania

Anopheles funestus is playing an increasing role in malaria transmission in parts of sub-Saharan Africa, where An. gambiae s.s. has been effectively controlled by long-lasting insecticidal nets. We investigated vector population bionomics, insecticide resistance and malaria transmission dynamics in 86 study clusters in North-West Tanzania. An. funestus s.l. represented 94.5% (4740/5016) of all vectors and was responsible for the majority of malaria transmission (96.5%), with a sporozoite rate of 3.4% and average monthly entomological inoculation rate (EIR) of 4.57 per house. Micro-geographical heterogeneity in species composition, abundance and transmission was observed across the study district in relation to key ecological differences between northern and southern clusters, with significantly higher densities, proportions and EIR of An. funestus s.l. collected from the South. An. gambiae s.l. (5.5%) density, principally An. arabiensis (81.1%) and An. gambiae s.s. (18.9%), was much lower and closely correlated with seasonal rainfall. Both An. funestus s.l. and An. gambiae s.l. were similarly resistant to alpha-cypermethrin and permethrin. Overexpression of CYP9K1, CYP6P3, CYP6P4 and CYP6M2 and high L1014S-kdr mutation frequency were detected in An. gambiae s.s. populations. Study findings highlight the urgent need for novel vector control tools to tackle persistent malaria transmission in the Lake Region of Tanzania.

The widespread deployment of primary vector control interventions, principally long-lasting insecticidal nets (LLINs) and indoor residual spraying (IRS), has substantially reduced malaria incidence across sub-Saharan Africa 1,2 . Between 2000 and 2015, 68% of the 1.5 billion malaria cases averted can be attributed to LLINs alone 1 . However, current estimates indicate the rates of decline have begun to stagnate 2 . Tanzania is among the 10 sub-Saharan African countries where malaria burden is concentrated 3 , contributing to 5% of global malaria deaths 2 . Malaria infection varies nationwide with an average prevalence of 7.3% in children under 5 years of age in 2017 4 . Vector control by the National Malaria Control Programme (NMCP) is based on sustaining high LLIN access and use 5 , via universal coverage campaigns supplemented with continuous distribution from school net programmes, antenatal care campaigns and the expanded programme for immunization; and targeted IRS in high transmission areas in the North-West 6 . Effective and sustainable malaria vector control is plagued by a number of challenges, including the evolution of vector behavioural and physiological resistance to current control interventions 7 . In the majority of sentinel districts across Tanzania, Anopheles mosquitoes have demonstrated reduced susceptibly to at least one public health insecticide 8,9 .

Results
Household characteristics. A total of 1593 households were visited during two cross-sectional entomological field surveys, across 86 clusters in Misungwi district, North-West Tanzania on the southern shore of Lake Victoria, between August and December 2018 (Fig. 1A).
The district spans two agro-ecological zones, based on vegetation land cover and rainfall ( Fig. 1), that are divided roughly into northern and southern clusters.
The study had an overall response rate of 86.2% (1372/1592), with consent to participate in the survey given from an adult/head of the household. Ten per cent (164) of dwellings were found vacant, 1.0% (16) were not located, 0.2% (3) not visited due to accessibility and 0.1% (2) were ineligible (no children under 15 years) during the survey period. A small proportion 2.2% (35) refused to participate in the study. The average altitude of study households was 1194.9 m above sea level (Table 1). Similar proportions of houses were classified as improved or unimproved, based on construction with modern or traditional materials, respectively (Table 1; Fig. 3A). Notably, few houses had mosquito proofing materials over the windows (29.9%) and almost no houses had ceilings (97.2%); 40.6% of houses had open eaves (Table 1; Fig. 3A,B). The average household size was 6.6 persons and mean number of room/sleeping place was 2.7 per house. Forty-four per cent of households owned at least one livestock (mostly goats and cattle), which were usually kept outdoors about 20 m away from the house.
LLIN ownership was very high in the study area with the majority of families owning at least one LLIN (94.9%); LLIN access was comparatively lower, however, the majority of households had enough LLINs to cover all of their sleeping places (62.3%). About 54.2% of households were sprayed during the 2015 IRS National Malaria Control Campaign (Table 1). There were no significant differences in household characteristics, including size, altitude, construction materials, however population access to insecticide-treated net (ITN) access was slightly higher in northern than southern clusters (Table 1).
Vector distribution, species composition, relative abundance and seasonality. A total of 23,081 mosquitoes, comprising 23.1% (5329) Anophelines and 76.9% (17,752) Culicines, were collected using Centers for Disease Control and Prevention light traps (CDC-LTs) during two cross-sectional survey rounds between  Table 1. Household characteristics in the study area in Northern and Southern clusters. The district spans two agro-ecological zones, based on vegetation land cover and rainfall ( Fig. 1), that are divided roughly into northern and southern clusters. HH household, ITN insecticide-treated net, N/A not applicable.   Table 2). Amongst sibling species of the An. gambiae complex, there were marked spatial and seasonal fluctuations. Most An. gambiae s.s. (94.3%) were collected from the northern zone and more than 71.5% of An. funestus s.s. from the southern part. An. arabiensis was distributed throughout both study zones, but its abundance was slightly higher in the South (Table 2). Both An. funestus s.s. and An. arabiensis predominated throughout the study period, but An. gambiae s.s. abundance peaked in December in the middle of the rainy season. An analysis of bioclimatic and landcover characteristics across the study area demonstrated several ecological differences between the northern and southern zones, with the former composed mostly of grassland and cropland (91%), with smaller proportions of shrubland and forest (7%) and areas prone to regularly flooding (1%); and the latter with less grassland and cropland (80%) and greater proportions of shrubland and forest (14%) and areas prone to regular flooding (3%) (Supplementary table S1). Furthermore, differences in rainfall were also observed between the northern and southern zones, with villages in the North receiving slightly higher average annual precipitation (959.5 mm) than in the South (911.5 mm).
While overall An. gambiae s.l. density was low, it was closely correlated with seasonal rainfall patterns. Mean An. gambiae s.l. caught per house during the dry season (August and September; average precipitation of 5-7 mm) was 0.14 but rose significantly by two-fold (DR = 1.73 [95% CI 1.08-2.78]; p = 0.02) in the wet season (October, November and December; average precipitation of 147.3-158.5 mm). By comparison, significantly higher An. funestus s.l. densities were observed during the dry months [mean = 4.61; Incidence Rate Ratio (IRR) = 1.58, p = 0.004] ( Table 3).
As summarized in Table 4, the greatest proportions of Anopheles were sampled by indoor CDC-LTs (48.4%) and outdoor tent traps (41.9%). An. arabiensis and An. gambiae s.s. showed similar tendencies of feeding both indoors (54.7% and 45.3% collected in CDC-LTs, respectively) and outdoors ( Phenotypic resistance and underlying molecular and metabolic resistance mechanisms. Wild populations of An. funestus s.l. and An. gambiae s.l. from across the study area were confirmed resistant to the diagnostic concentration of alpha-cypermethrin, with mean 30-min knock-down ranging from 43.7 to 59.4% (Table 5). Similarly, both species were resistant to permethrin, with average 24-h mortality ranging between 38.3 and 56.5% (Table 6). In general, levels of resistance to both pyrethroids were comparable between An. gambiae s.l. and An. funestus s.l., as well as between northern and southern study zones ( Fig. 2C; Tables 5, 6).
Overall, the majority 91.9% (565/615) of An. funestus s.l. mosquitoes tested in bioassays were confirmed by PCR as An. funestus s.s, with a small proportion of An. parensis (7.8%; 48/615). An. gambiae s.l. from bioassays  Mean mortality 24 h after exposure to the standard diagnostic dose of alpha-cypermethrin was predicted for An. gambiae s.l. in 2017 using a geospatial model. The model used data from WHO susceptibility tests conducted from 2005 to 2017 and incorporated associations between resistance and potential explanatory variables such as ITN coverage using three different machine learning approaches. Predicted mean mortality in An. gambiae s.l. for each 5 × 5 km square (Fig. 2C) was high across Tanzania in 2017. Within Misungwi, the lowest mortalities/ highest resistance occurred in the West and North-West.

Discussion
Despite substantial gains achieved in malaria control across Tanzania over the past 20 years, attributable to improved quality and access to diagnostics and treatment and the widespread scale-up of LLINs and targeted IRS, localised transmission persists, especially in the Lake Region. Study findings demonstrated that An. funestus s.l. is becoming a dominant, efficient malaria vector species in Misungwi district, North-West Tanzania in an area with high coverage of standard pyrethroid LLINs and historical IRS activities. A similar phenomenon has  www.nature.com/scientificreports/ recently been reported from south-eastern Tanzania 25 ; however, our study indicates this shift in species composition may not be restricted to the south of the country. Around Lake Victoria, species abundance and transmission intensity vary quite considerably spatially and temporally 28 , with implications for the deployment of effective malaria vector control interventions. These heterogeneities likely reflect differences in climatic conditions such as rainfall and ecological settings, which support the breeding of particular vector species 29 . In our study, overall vector densities were significantly higher in villages located in the southern part of the study district compared to the northern clusters. An. gambiae s.s. and An. arabiensis occurred across Misungwi district, however, An. gambiae s.s. abundance was lower in the South and concentrated mostly in the North. By comparison, An. funestus s.s. was equally distributed throughout the district in sympatry with An. arabiensis and An. gambiae s.s. but found at the highest densities along shorelines and waterways feeding into Lake Victoria. The spatial variation of Anopheles sibling species may be explained by several factors linked to ecological features, including turbidity, water quality, relative humidity, temperature, vegetation type and/or socioeconomic parameters, such as ownership and usage of insecticide-based vector control measures and livestock density, as observed in previous studies conducted on the Kenyan side of Lake Victoria 34 . An. gambiae s.s. are known to breed in rain-dependent temporary habitats 35 , while An. funestus s.s. and An. arabiensis can colonize large permanent aquatic habitats, some with large vegetation, in arid and highland areas 36,37 . Most residents in Misungwi district, especially in the southern clusters, traditionally stored rainwater for domestic purposes and animal husbandry in large, permanent man-made pools, locally called "Rambo", which were filled throughout the year and could serve as potential breeding sites for An. funestus s.s. and An. arabiensis, even in the dry season; the higher density of An. gambiae s.s. during the rainy season is likely due to increased availability of temporary breeding sites 38,39 . In addition, agricultural practices such as irrigated rice paddies create diverse aquatic mosquito breeding habitats that could influence the co-existence and abundance of different vector species in the study area 40,41 . An. funestus s.l. was collected in both seasons but peaked during the dry season, consistent with its ability to develop in habitats that can sustain desiccation 42 . Of concern, both An. gambiae s.l. and An. funestus s.l. malaria vector species were present during different seasons, favoured by distinct climatic and ecological conditions, sustaining malaria transmission throughout the region and the year. Of all mosquitoes sampled in this study, Culex species were the most abundant, contributing to 67.6% of the www.nature.com/scientificreports/ indoor host-seeking population. Previous studies in Tanzania have highlighted that Culex species, commonly referred to as "the house mosquito", predominant in malaria-endemic communities 43 , and when resistant to public health insecticides, can jeopardise community adherence to vector control interventions, due to perceived failure of these strategies 44,45 . This study estimated that each household could receive an average of more than 53 infective bites per year from both major vector species (An. funestus s.s. and An. gambiae s.s.) despite high coverage of LLINs. Comparably high EIRs have also been reported from other rural and peri-urban regions of East Africa, including South-Central Tanzania 46 , coastal Kenya 47 and South-West Ethiopia 48 . Furthermore, the annual EIR was ten times higher in villages located in the southern part of the study district compared to the North. In northern clusters, where An. gambie s.s. and An. funestus s.s. co-existed, even though An. gambiae s.s. was present in very low numbers, these two species generally had equivalent Plasmodium infection rates. Both An. gambiae s.s. and An. funestus s.s. can be highly anthropophilic and endophilic, but the former species may be more aggressive and efficient vector in terms of malaria transmission, possibly due to host competition. In the southern study clusters, malaria transmission was almost exclusively mediated by An. funestus s.s. Only a single P. falciparum-infected An. arabiensis was collected during the study which might be explained by its highly opportunistic behaviour, feeding on both animals and humans; in the absence of the latter host it can display strongly zoophilic feeding preferences for livestock, of which close to half of the households owned 49,50 . This species is also known for its more exophilic tendencies compared to An. gambiae s.s. 50,51 and can easily adapt and feed outdoors in response to insecticidal interventions 52 , especially when human or animal populations are available outside 12 . Our indoor and outdoor collections generally support these behavioural assumptions, with both An. gambiae s.s. and An. funestus s.s. sampled in similar proportions across different traps, with the exception of An. gambiae s.s., which was found at very low densities resting outdoors. The occurrence of highly endophilic and anthropophagic vectors such as An. funestus s.s. host-seeking or resting outdoors could be linked to behavioural divergence among vector populations and/or chromosomal inversion polymorphisms 53,54 , as well as human behavioural changes 55 . However, in our study, more sporozoite-harbouring mosquitoes were collected in houses, suggesting ongoing malaria transmission is still occurring inside, despite high intervention coverage. The strongly exophilic behaviour of An. arabiensis indicated that LLINs and IRS in Misungwi district may have minimal effects against this species; although its relative importance in local malaria transmission appears diminished. Moreover, the presence of infected An. gambiae s.s. and An. funestus s.s. outdoors, coupled with a degree of exophagic behaviour, suggests that additional control tools targeting outdoor vector populations may warrant consideration in the study area 56 .
All three major vector species (An. gambiae s.s., An. arabiensis and An. funestus s.s.), displayed low levels of susceptibility to alpha-cypermethrin and permethrin, primarily due to selection pressure from prolonged use of pyrethroid-based LLINs and likely enhanced by concurrent agricultural pesticide application [57][58][59] . Because vectors for resistance monitoring were collected as adults of unknown age from house walls, and phenotypic resistance declines with age 60,61 , it is probable that our measurements of pyrethroid resistance may be an underestimation. Study findings align with others in the Lake Zone and across Tanzania, demonstrating low mortality among vector populations to the diagnostic doses of pyrethroids 9,58,62 . Data collected across Africa indicated that previously pyrethroid resistance was higher in East African populations of An. funestus s.l. compared to East African populations of An. gambiae s.l. up to 2017, but this difference was not detected in Misungwi in 2018-19 63,64 . It is noteworthy that our bioassays presented higher levels of resistance in An. gambiae s.l. in comparison to those shown in our map produced by geospatial models of phenotypic resistance to alpha-cypermethrin for the year 2017, suggesting pyrethroid resistance may be increasing in An. gambiae s.l. populations in the region. In Misungwi district, population-level frequency of the L1014S-kdr mutation was practically fixed in An. gambiae s.s., while CYP6M2, CYP6P3, CYP6P4 and CYP9K1 were modestly upregulated by comparison to reports from West Africa [65][66][67] . These results indicate that both target site and metabolic mechanisms may be driving pyrethroid resistance in An. gambiae s.s. in this study area. However, further investigation is necessary to identify resistance mechanisms specific to these field populations, including those in An. funestus s.l., which were unavailable to us at the time of this study, for prospective monitoring and to improve our understanding of the specificity of resistance mechanisms to individual interventions and the likelihood of selecting for cross-resistance between active ingredients [65][66][67] . Importantly, intense insecticide resistance may partially explain the persistent malaria transmission in Misungwi district, highlighting the urgent need for novel vector control tools, containing different insecticide classes and combinations. This study was undertaken prior to the phase III evaluation of novel bi-treated LLINs containing a pyrethroid and either a pyrrole (chlorfenapyr), a synergist (piperonyl butoxide; PBO) or a juvenile growth hormone inhibitor (pyriproxyfen; PPF) 31 , which may have the potential to control malaria transmitted by pyrethroid-resistant vector species.
This study was conducted to characterize baseline vector population bionomics and malaria transmission dynamics in Misungwi district, with some limitations. Because mosquito collections spanned five months of the year, encompassing the short rainy season (October-December 2018) and part of the dry season (August and September 2018), vector densities, sibling species composition and sporozoite rates reported in this study may not be representative of the annual variation in this area. Additional studies are ongoing investigating indepth the association between vector spatial distributions and key ecological indices, and to identify insecticide resistance mechanisms in An. funestus s.l., which at the time of study design, was not anticipated to emerge as the major vector species in Misungwi district.

Conclusion
In Misungwi district, North-West Tanzania, An. funestus s.s. is the leading malaria vector species, predominating in southern villages of the study site, across dry and wet seasons. An. gambiae s.s. was present in much lower densities, concentrated mostly in the North during the wet season, potential driving malaria epidemics. Annual www.nature.com/scientificreports/ EIR was high, despite high LLIN usage, but variable within a small geographical area, influenced by vector species diversity and bionomics, with serious epidemiological implications for malaria control. An. gambiae s.s., An. arabiensis and An. funestus s.s. were found similarly resistant to pyrethroids, with high frequencies of target site alleles and overexpression of detoxification genes identified in An. gambiae s.s. Study findings highlight the urgent need for novel vector controls strategies, which incorporate new chemical classes, to control malaria transmitted by pyrethroid-resistant vector populations and sustain gains in malaria control across the Lake Region.

Methods
Study area characteristics. The study was carried in Misungwi district (latitude 2.85000 S, longitude 33.08333 E) in North-West Tanzania on the southern shore of Lake Victoria (Fig. 1A). Misungwi lies at an altitude of 1150 m above sea level, with a population of approximately 351,607 according to the National population and housing census of 2012 68 . The district experiences a dry season typically between June and September and two annual rainy seasons; the long-rainy season between February and May and a short-rainy season between November and December with average annual rainfall ranging between 0.5 and 58.8 mm. The district is geographically divided into two main agro-ecological zones (northern and southern zones), based on the vegetation land cover and amount of rainfall. The local communities practice rice, millet and cotton farming, domestic animal rearing, fishing and have small-scale businesses as sources of income and food. In preparation for the RCT, the study area was sub-divided into 86 clusters, containing 72 villages made up of 453 hamlets from 17 wards (Fig. 2). Detailed criteria and methodology used for cluster formation is described elsewhere 31 . Across Misungwi district, a typical compound was comprised of a house and cattle shed. Houses were generally constructed from both modern and traditional materials and most houses had eave spaces (an opening between the wall and the roof for ventilation) (Fig. 3A,B). The area experiences moderate to high malaria transmission and malaria incidence peaks shortly after the rainy seasons 69 . Recent studies conducted in 2010 and 2017 reported a malaria prevalence of 51.3% across all age groups and 46.3% in school children (7-14 years). LLINs mainly Olyset and PermaNet 2.0 obtained through national bed net distribution campaigns have been the primary malaria control method in the study area 4,70 and IRS was last conducted in this region in 2015. A preliminary survey carried out by our study team in 2018 found An. gambiae s.s. An. arabiensis and An. funestus s.s. as the predominant malaria vector species in the study area.

Environmental characteristics of northern and southern clusters.
To characterize the study area with regards to climatic and environmental conditions, high spatial resolution satellite remote sensing and other geospatial data were downloaded in raster (i.e. gridded) format from publicly available data sources and processed using ArcGIS 10.5.1 (ESRI, Redlands, USA). Data on eight bioclimatic variables at 1 km 2 spatial resolution, representing averages for the years 1970-2000, were obtained from the WorldClim 2 database 71 : annual mean temperature (Bio1), temperature seasonality (Bio4), maximum and minimum temperature of the warmest month (Bio 5 and 6), annual precipitation (Bio12), precipitation of the wettest and driest months (Bio13 and 14), and precipitation seasonality (Bio15). Global elevation data were obtained for the study area from NASA's Shuttle Radar Topography Mission (SRTM) 4.1 at 90-m spatial resolution 72 . Global landcover data were obtained from the European Space Agency GlobCover 2009 project, available at 300-m spatial resolution (© ESA 2010 and UCLouvain; http:// due. esrin. esa. int/ page_ globc over. php). These data identify 22 landcover types, of which 12 were identified in the study area. Zonal mean statistics for the northern and southern clusters were calculated for each bioclimatic variable and elevation using the spatial analyst toolbox in ArcGIS; cluster means were then averaged for each zone (Supplementary Table S1). The proportion of cells within each of the northern and southern zones that were classified as different landcover types were similarly calculated; Supplementary Table S1 shows the results for the five dominant landcover types that were present in the study area [mosaic vegetation (i.e. grassland/shrubland/forest: 50-70%/cropland: 20-50%), herbaceous vegetation (i.e. grassland/savannas), shrubland, broadleaf deciduous forest/woodland, and grassland or woody vegetation on regularly flooded or waterlogged soil], which represent 98.9% and 97.5% of the total area of the northern and southern clusters, respectively. While the available data sources, and hence these estimates, are derived from time periods prior to the study period, we assume that the estimates accurately reflect relative differences in climatic and environmental conditions across the study area.
Indoor and outdoor entomological surveillance. Two cross-sectional entomological field surveys were conducted between August and December 2018 in all 86 clusters, using CDC-LTs (John W Hock Company, USA). Eight households were randomly selected from a census list of households generated during baseline enumeration. CDC-LTs were hung next to the feet of an occupant sleeping under an ITN/untreated net (about 1 m from the ground), between 19:00 and 7:00. A questionnaire was administered to the head of the household to gather information related to the house structure (type of wall and roofing materials, window screening, number of rooms, number of sleeping places, presence of eaves), and coverage and usage of LLINs/untreated nets. Direct observation was also used during data collection to validate participant answers.
Assessment of Anopheles vector feeding and resting behaviours indoors and outdoors was undertaken in 48 clusters between December 2018 and January 2019. Two households per cluster were randomly selected, and each house was installed with a CDC-LT indoors and an occupied Furvela tent trap outdoors 73 . Both Furvela and CDC-LTs were switched on at 19:00 and off at 7:00. Indoor and outdoor resting adult Anopheles were collected from the same houses using a 12 voltage battery-powered Prokopack aspirator 74 , and manual mouth aspirators 75,76 . Systematic sampling of adult resting mosquitoes on the walls, roofs and floors were conducted for up to three minutes depending on the size of the room. Outdoor collections were performed from potential resting sites around the house, such as open resting structures, cow sheds and pit latrines. An. funestus s.l. were assessed in six clusters selected on the basis of high Anopheles population densities. Due to difficulties locating reliable, productive breeding sites for larval sampling, particularly for An. funestus s.l., instead adult female Anopheles resting indoors were collected using both Prokopack and mouth aspirators for resistance testing 74 . Mosquitoes were separated by species complex and supplied with 10% glucose solution for 72 h to allow digestion of blood meal, prior to bioassay testing with permethrin and alpha-cypermethrin. These two pyrethroids were selected on the basis that they are the partner insecticides in next-generation LLINs (Olyset Plus, Royal Guard and Interceptor G2), which were about to be distributed across the study district 77 .
In WHO tube assays, 20-25 gravid female An. gambiae s.l. or An. funestus s.l. of unknown age were exposed to 0.75% permethrin for 60 min 78 . In CDC bottle bioassays, 20-25 gravid female An. gambiae s.l. or An. funestus s.l. of unknown age were exposed to 12.5 μg/ml alpha-cypermethrin for 30 min 79 . For both assays, knock-down was recorded at the end of the diagnostic exposure time (30 or 60 min after exposure, for CDC or WHO bioassays, respectively), and final mortality was scored after 24 h. All mosquitoes tested in bioassays were stored individually for sibling species identification.
Mosquito processing, species identification and sporozoite detection. Adult female mosquitoes collected from the cross-sectional surveys, resistance bioassays and behaviour study were sorted and identified based on their morphology, separating An. gambiae s.l. from An. funestus s.l. and from other genera according to Gillies and Coetzee 80 . At least three female An. gambiae s.l. and three An. funestus s.l. per household/ per collection method was analysed for presence of Plasmodium falciparum CSP using enzyme-linked immunosorbent assay (CSP-ELISA) 81 . All CSP positive samples and a sub-sample of CSP negative An. gambiae s.l. and An. funestus s.l. from the cross-sectional surveys, resistance bioassays and behaviour study, were randomly picked per house and tested for species identification. DNA was extracted from legs/wings and TaqMan PCR assays were performed to distinguish sibling species in An. gambiae 82 or An. funestus complexes 83 .

Identification of insecticide resistance mechanisms.
A subsample of An. gambiae s.s. and An. arabiensis were genotyped for L1014F-kdr and L1014S-kdr mutations associated with pyrethroid and DDT resistance, using TaqMan PCR assays 84 . Blood-fed indoor resting female adult mosquitoes (F0s) were collected using mouth aspirators from three wards, sampled for phenotypic resistance testing. Mosquitoes were held for 3-4 days to allow for blood meal digestion. Individual An. gambiae s.l. were placed into Eppendorf tubes containing moist filter papers and forced to lay eggs, as previously described 23 . The first 3-5 day emerged F1 adults from each parent were stored individually in RNAlater and preserved at − 20 °C for gene expression analysis. Expression profiles for metabolic detoxification genes in a sub-sample of 250-300 F1 wild-caught female An. gambiae s.s. mosquitoes were determined using quantitative reverse transcriptase PCR (qRT-PCR) 85,86 . Individual mosquitoes were first tested for species identification, and only mosquitoes identified as An. gambiae s.s. were analyzed. A minimum of 5 pools of 10 An. gambiae s.s. were analysed for CYP6M2, CYP6P3, CYP6P4 and CYP9K1 gene expression 85 .
Data analysis. Field data were entered into an Open Data Kit (ODK) form. Data analysis was performed using Stata/IC 15.1 (Stata Corp., College Station, USA) 87 . Mean Anopheles caught per night per household, sporozoite rate and their 95% confidence intervals (CIs) were estimated according to study zones (North or South), season (wet or dry) and Anopheles species. Anopheles vector population density and EIR were analysed and compared between study zones and seasons using multilevel negative binomial regression taking into account clustering effect. The EIR was calculated at household level as the average number of CSP-ELISA positive mosquitoes per night and was weighted to account for the proportion of collected Anopheles processed for CSP-ELISA. Logistic regression was used to compare sporozoite rates between the two study zones and seasons. The proportion of households with at least one LLIN was computed by dividing total nets observed and recorded by total households surveyed. Net access was estimated from the proportion of households with enough LLINs over total sleeping places. Unimproved houses were classified as houses with open eaves, unscreened windows, and were constructed from traditional low-quality materials such as a thatched roof, mud and non-plastered walls. Improved houses had closed eaves, with mosquito proof mesh on the windows, and were built with improved modern materials such as an iron sheet as a roof, brick/blocks, with plastered walls. For resistance phenotyping, mean percentage knock-down/mortality post-exposure was calculated and interpreted following the WHO and CDC criteria 78,79 . Susceptibility thresholds were considered at the diagnostic time (24 h and 30 min post-exposure for WHO and CDC bioassays, respectively). Mean mosquito mortality between 98 and 100% indicated full insecticide susceptibility, 90-97% showed suspected resistance that needed confirmation and less than 90% indicated confirmed resistance 78,79 . When control mortality was between 5 and 20%, results were corrected using Abbott's formula 78 . If the control mortality was ≥ 20%, the test was discarded 78 .
Gene expression and fold-change, relative to the susceptible laboratory strain An. gambiae s.s. Kisumu were calculated according to the 2-ΔΔCq method 88 after standardisation with housekeeping genes (elongation factor; EF and 40S ribosomal protein S7; RPS7).
To estimate spatial and temporal trends in pyrethroid resistance using the available field data, which are sparse and have a heterogeneous distribution, a total of 6423 observations of mortality from WHO susceptibility tests from 2005 to 2017 were used to inform a Bayesian geostatistical ensemble model. The model was also informed by a suite of 111 explanatory variables describing potential drivers of selection for resistance such as ITN coverage and produced estimates of mean mortality in An. gambiae s.l. for two regions of Africa 64 . Here we present the results for alpha-cypermethrin resistance in Misungwi in 2017.