Temporally consistent predominance and distribution of secondary malaria vectors in the Anopheles community of the upper Zambezi floodplain

Regional optimisation of malaria vector control approaches requires detailed understanding both of the species composition of Anopheles mosquito communities, and how they vary over spatial and temporal scales. Knowledge of vector community dynamics is particularly important in settings where ecohydrological conditions fluctuate seasonally and inter-annually, such as the Barotse floodplain of the upper Zambezi river. DNA barcoding of anopheline larvae sampled in the 2019 wet season revealed the predominance of secondary vector species, with An. coustani comprising > 80% of sampled larvae and distributed ubiquitously across all ecological zones. Extensive larval sampling, plus a smaller survey of adult mosquitoes, identified geographic clusters of primary vectors, but represented only 2% of anopheline larvae. Comparisons with larval surveys in 2017/2018 and a contemporaneous independent 5-year dataset from adult trapping corroborated this paucity of primary vectors across years, and the consistent numerical dominance of An. coustani and other secondary vectors in both dry and wet seasons, despite substantial inter-annual variation in hydrological conditions. This marked temporal consistency of spatial distribution and anopheline community composition presents an opportunity to target predominant secondary vectors outdoors. Larval source management should be considered, alongside prevalent indoor-based approaches, amongst a diversification of vector control approaches to more effectively combat residual malaria transmission.

Genetic identification: larvae. Over 50% of field-surveyed anopheline larvae were collected for genetic analysis (n = 855/1557), and a viable DNA sequence from COI and/or ITS2 regions was obtained from 658 specimens ( Table 1; the remainder did not amplify using COI or ITS2 primers, potentially due to sub-optimal preservation of DNA in the field). More than two-thirds of these specimens (n = 443) were identified to species based on both COI and ITS2 sequences with ≥ 95% similarity to GenBank accessions. An additional 20% of specimens (n = 129) matched with ≥ 95% identity with GenBank COI sequences, but their corresponding ITS2 sequences returned no match at ≥ 95% similarity with a GenBank sequence. Each remaining specimen (n = 86) was assigned a species identity based on an above-threshold match of either its COI or ITS2 sequence to Gen-Bank accessions, with the exception of 18 specimens (2.7% of the total) whose closest matches were to Anopheles sequences but at similarities below the threshold; these 18 were subsequently designated 'unknown Anopheles species' . Maximum likelihood (ML) tree analysis of the 119 resolved COI haplotypes corroborated the species identities assigned from matches to GenBank sequences by demonstrating statistically well-supported clusters nested with reference specimen sequences of known morphological identity ( Fig. 3; Supplementary Fig. S1). An. gambiae complex species clustered together but were not resolved into species from COI phylogeny.
Within the sample of 658 sequenced larvae, a large majority were identified as An. coustani (80.1%; Table 1; Fig. 4). Two distinct genetic clades could be further separated within An. coustani; the first consisted of specimens www.nature.com/scientificreports/ with both COI and ITS2 sequences matching with ≥ 95% identity to accessions for the species, whilst the second clade (comprising 22% of the total sample) consisted of specimens with above-threshold matches to COI reference sequences but below-threshold matches to any GenBank ITS2 sequence. This second An. coustani clade also forms a well-supported cluster in the ML tree ( Fig. 3; for full ML tree see Supplementary Fig. S1). An. pharoensis was the next most abundant species at 10.9% (Fig. 4), followed by small numbers of An. gambiae sensu lato (2.1%), An. squamosus (2.0%), An. species O/15 23 (1.7%), and 3 specimens (0.5%) designated as An. species UG3, signifying 'unknown group 3' based on similarity to designations made in southern and central Africa by Ciubotariu et al. 44 . An. gambiae sensu lato could be further resolved into An. arabiensis (1.1%) and An. gambiae sensu stricto (1.1%, hereafter An. gambiae) based on ITS2 sequences.
Geographic distribution of larvae. An. coustani was distributed ubiquitously across ecological zones ( Fig. 5c). A significantly lower median per transect point was found in dambos than in the Zambezi and Luena floodplains and the Lui valley, and a significantly higher median in the Lui valley compared to the Zambezi floodplain and its edge ( Fig. 5a; ISM test 26.150, df 4, p < 0.001; SSD adj p < 0.05). Within An. coustani, however, clade 1 exhibited a significantly lower median in the Lui valley (and dambos; median = 0) than all other habitats ( Fig. 6; 0.6-1.592; ISM test 20.036, df 4, p < 0.001; SSD adj p < 0.05). Conversely, the clade 2 median was significantly higher in the Lui valley than all other ecological zones ( Fig. 6; median 1.822 compared to 0-0.369; ISM test 35.518, df 4, p < 0.001; SSD adj p < 0.05). The median proportion of estimated An. coustani populations composed of clade 2 individuals was significantly higher in the Lui valley than all ecological zones except dambos, and the proportion in dambo populations exceeded those of Zambezi floodplain and floodplain edge populations ( Fig. 6; ISM test 34.057, df 4, p < 0.001; SSD adj p < 0.05).
An. pharoensis was also encountered in every ecological zone (Fig. 5d), albeit in fewer transect points (21%) than An. coustani (65%). Its distribution in Zambezi and Luena floodplains was significantly different from dambo and floodplain edge habitats ( Fig. 5b; Kruskal-Wallis H = 21.392, p < 0.001; mean rank former > latter, SSD adj p < 0.05). The primary vector complex An. gambiae s.l. was found in floodplain edge (n = 2 An. gambiae),   Fig. 5e). An. squamosus was found in all ecological zones except for dambos (Fig. 5d), and small numbers of An. species O/15 were present in every ecological zone, occurring in 5% of wet transect points; there was no significant difference in distributions between ecological zones for either of these species. Three specimens of An. species UG3 were found, from samples in floodplain edge and Luena flats ecological zones (Fig. 5e) www.nature.com/scientificreports/ end 10:11 ± 0h40). One hundred and twenty-six trap nights yielded 58 anophelines which were predominantly female (93%), with a considerable majority trapped outside (71%), including all males (n = 3; excluded from subsequent analyses). The distribution of catches amongst the traps was highly clustered; 32 traps contained one or more anophelines, and a single outside trap accounted for 25% of the total catch. Although modest in extent and compromised by the unusually dry wet season, we report and interpret these results as they coincided geographically and temporally with larval surveys. Fifty-four adult female specimens were identified to 6 species from above-threshold matches to COI and/ or ITS2 DNA sequences (Table 2), confirmed by high bootstrap support of clustering within phylogenetic ML tree analysis ( Supplementary Fig. S1). An. species UG3 dominated the adult sample (46%), with the majority www.nature.com/scientificreports/ found outdoors, although a single outdoor trap accounted for half of this outdoor total. Over half of An. species UG3 specimens were identified by matches of both COI and ITS2 sequences to accessions linked to Jones 45 and Ciubotariu et al. 44 , respectively, whilst the remaining specimens were identified from matches to either these ITS2 (24%) or COI (20%) sequences. ITS2 sequences permitted differentiation between species within the An. funestus group that could not be resolved based solely on COI sequences. The typically anthropophilic primary malaria vector species An. funestus (24% of sample) was predominantly trapped indoors (77% of specimens). An. gambiae adults were found exclusively indoors, whilst An. coustani, An. arabiensis and An. rivulorum (An. funestus group member) were found exclusively outdoors. One specimen was assigned 'unknown An. species' . None of the adults were bloodfed or reported positive amplifications from PCR aimed at detecting P. falciparum sporozoites.
Geographic distribution of adult mosquitoes. Adult anopheline catches exhibited a high degree of within-village similarity, whilst varying substantially between villages (Fig. 7). Although An. species UG3 numerically accounted for almost half of the adults sampled, this was a consequence of an unusually high catch in one dambo village; it was found in one-third of villages where Anopheles adults were caught (n = 4 of 12). Samples from two villages in the Luena flats and one in a dambo were exclusively An. funestus, while the species was also found on the floodplain edge and the Lui valley. An. coustani was found in 4 villages, whilst An. gambiae adults were confined to 2 villages in the Lui valley; An. arabiensis and An. rivulorum occurred solely in one floodplain edge and one Zambezi floodplain village, respectively (Fig. 7).   www.nature.com/scientificreports/ parisons between anopheline species distributions in the 2017 dry season (larvae 36 ), after the peak of the 2018 wet season (larvae 36 ), and after the peak of the 2019 wet season (adults and larvae) revealed considerable interannual consistency in the spatial distribution of some anopheline species where sampling effort coincided. For example, An. coustani adults were caught in 6 outdoor traps in 2019 (Fig. 7), and larvae of this species were found in the nearest wet transect point to each of these traps (Fig. 5c)   www.nature.com/scientificreports/ pling periods from 2014 to 2019 were identified as An. coustani, ranging from 49 to 100% of the sample, whilst 13% of the specimens across all periods were An. squamosus (Fig. 8). An. gambiae s.l. mosquitoes were the sole primary vector species detected, occurring only in 2017 when they constituted < 0.5% of the sample. Although anopheline species count data from these studies 37,38 were not published, the data collected during the studies and subsequently (Orba et al., 2021; pers. comm.) are provided in the Supplementary Information to this paper (Supplementary Table 1), along with six COI gene sequences (Supplementary Dataset S1). There was strong support for the morphological identification of adult An. coustani as the DNA sequences obtained from reference samples within the independent dataset (Orba et al., 2021, pers. comm.) clustered with sequences on GenBank and An. coustani sequences from the present study ( Fig. 3; Supplementary Fig. S1). The time series of data indicates the consistent dominance of An. coustani in the region in both dry and wet seasons (Fig. 9).

Discussion
Geographically extensive larval mosquito sampling was undertaken across multiple ecological zones in western Zambia over three years, including one dry and two wet seasons, and data were integrated with adult mosquito catches from this study and an independent series of collections over a longer timeframe 37,38 (Orba et al., 2021, pers. comm.). This combined dataset revealed a marked temporal consistency in the composition of the Anopheles mosquito community across western Zambia, and also consistency in geographic distribution of species across larval habitats in the region. This spatial consistency is evident both at landscape scale, across large ecological zones in a study area over 3500 km 2 , and also at fine scale whereby larvae were frequently encountered within 10-100 s of metres of congruent records from a preceding study in the region 36 . The community consistency was even more notable given the significant variation in ecohydrological conditions between the 2018 wet season and the 2019 season surveyed in the present study. Barotseland experienced an anomalously early and unusually dry wet season in 2019, such that widespread drought in the ensuing dry season was the most intense and severe for at least 20 years 39 , whilst Zambezi River discharge in 2018 was the second highest recorded for almost 30 years 36 . Extensive surveys of Anopheles larvae across diverse habitats in this study, coupled with objective molecular identification to species revealed a paucity of primary vector species, with support for this finding provided by an independent dataset of adult mosquitoes caught in the study area in the same time period 37 S1).
An. coustani was found to be by far the most ubiquitous and abundant species across the region, constituting an estimated 78% of all surveyed anopheline larvae (Fig. 4), compared with 52% of surveyed anophelines in the preceding 2018 wet season and 65% in the dry season of 2017 which comprised An. coustani group larvae 36 . Although An. coustani is often rare in entomological surveys 24 , in recent studies it has been found to dominate local anopheline communities at locations in the western Kenyan highlands 46 and Madagascar 47 , and was the second most prevalent species at locations in Ethiopia 48 and rural southeast Zambia 49 . The species was distributed ubiquitously across our study area, but with higher abundance in the Luena flats corresponding to a significantly Adult An. coustani group mosquitoes maintain a year-round presence in Kenya 53 and Ethiopia, although they reached peak density in the drier months following the rainy season in the latter 48,54 ; An. coustani abundance may have been higher in the current study than in surveys following the 2018 wet season 36 due to the drier conditions experienced in western Zambia in 2019 39 . Although the long-term data from adult trapping in the centre of our study area (Fig. 9) have a limited seasonal resolution, they indicate the persistent presence of the An. coustani; whilst abundance fluctuated between years, An. coustani was consistently predominant, providing strong corroboration for the findings from larval surveillance. The An. coustani group includes multiple species, including An. coustani, An. crypticus, An. tenebrosus and An. ziemanni 51,55,56 ; some are morphologically indistinguishable in some life stages, and the use of molecular methods are advocated particularly to differentiate the cryptic An. crypticus from An. coustani 57 . While COI and ITS2 sequences have been published for the other species, none are available for An. crypticus. In the current study corresponding ITS2 sequences were obtained from 99% of specimens identified as An. coustani based on their COI sequences. This ITS2 dataset revealed two distinct genetic clades within An. coustani in western Zambia: one whose ITS2 sequences matched GenBank An. coustani sequences with > 95% similarity, and another with ITS2 similarity to GenBank accessions of < 95%. There was a high level of spatial overlap between clades, but significant differences in the proportion of the two An. coustani clades between ecological zones raises the possibility that the clades occupy different ecological niches within the study area. Ciubotariu and colleagues 44 also report this phenomenon amongst morphologically-identified adults, reporting a clade whose closest match for known An. coustani ITS2 sequences was ~ 80%. Members of both clades in that study were also reported to overlap spatially in northern Zambia and neighbouring D.R. Congo, and contained specimens positive for human blood meals and P. falciparum confirming their roles as secondary malaria vectors 44 . It is possible that An. coustani clade 2 in this study represents An. crypticus, although we note that our sequences do not match the ITS2 sequences of two 'coustani-like' species reported by Lobo et al. 22 which could also represent the species.
The presence of An. pharoensis in 2019 (reported here) despite its absence in 2018 36 is most likely attributable to a difference in genetic resolution between studies rather than a demographic/biological signal, as species assignment in the previous study was based almost entirely on COI which confers less discriminatory power between An. pharoensis and An. squamosus than ITS2 21,24,58 . The more extensive ITS2 sequencing performed here therefore permitted identification of An. pharoensis, and was strongly supported by phylogeny ( Fig. 3; Supplementary Fig. S1). This artefact does not affect our main findings, as pooling of An. squamosus with An. pharoensis in 2019 did not alter the ecological zone partitioning exhibited by An. pharoensis alone. Significantly www.nature.com/scientificreports/ higher An. pharoensis abundance in Zambezi and Luena floodplain ecological zones in 2019 corresponded with significantly higher prevalence of An. squamosus (potentially including An pharoensis individuals) in these habitats and along the floodplain edge in 2018 36 .
Eight of the nine anopheline species identified in this study have previously been implicated as malaria vectors, whilst the vector status of An. species UG3 and of our An. coustani clade 2 is unknown. In Kenya at least half of 17 species in one study tested positive for P. falciparum or were known vectors 23 , while 12/21 species in a more recent study were potential vectors 24 , and 12/18 in eastern Zambia 22 tested positive for P. falciparum. The number of adult anophelines sampled in the present study was severely limited by the unseasonally dry conditions, and the small sample did not contain bloodfed mosquitoes. Nonetheless, our adult samples were captured in traps baited with CO 2 and lures emulating human odours, so it can be assumed that they were seeking bloodmeals. An. coustani has tested positive for P. falciparum in multiple settings 4,22,46,47,[59][60][61][62] , and marked tendencies for exophagy 46,53,54,59 (as evident in this study) and early biting 48,54 may render it largely unaffected by indoor-focused interventions. Parasite-positive An. pharoensis is implicated in maintaining dry season malaria transmission in Ethiopia 54 and this species also favours outdoor biting 54,59 , as does An. squamosus 63 which was found to have a Plasmodium positivity rate equal to that of primary vectors on Madagascar 47 and an unexpected degree of anthropophily in southern Zambia 64 . Recently An. species O/15 was shown for the first time to harbour malaria parasites following outdoor trapping in the western Kenyan highlands 24 . An. rivulorum is considered to be a less efficient vector than its sibling species An. funestus, but nonetheless carries the malaria parasite 22,65 , whilst An. arabiensis, An. gambiae and An. funestus are classical primary vector species 56 , and the latter two were caught predominantly indoors in the present study.
Comparatively few studies have examined the inter-annual dynamics of anopheline communities across the landscape scales at which larval source management would need to be implemented 66 , given the dispersal capacity of Anopheles vector species 67 . Studies frequently seek to quantify the effect of seasonality on vector communities 48,53,54,68,69 , examine a limited number of locations at high temporal resolution 9,70 or investigate community change in response to an explicit intervention 10,71 . Spatially explicit descriptions of community composition and distribution at relatively large scales of space and time are scarce; previous studies have largely focused on adult mosquitoes and variously reported temporal shifts 9 or consistency 70,72,73 . The current study reports a high degree of consistency in anopheline community composition and distribution over three consecutive years. An. coustani larvae are shown to be widely distributed and to predominate in the community in both dry and wet seasons, while the traditional primary vector species An. arabiensis and An. gambiae comprise a consistently low proportion of the community (< 3% of larvae) and manifest a distinct small scale geographic clustering which is consistent across years. Due to the limited numbers (n = 54), the adult sample is unlikely to be representative of the anopheline community, particularly given the absence of samples from 30% of larval sampling locations and the skewing of the sample by a single anomalous catch constituting 25% of sampled adults. Nonetheless, the adult sample corroborates the An. gambiae complex hotspot in the Lui valley. Despite the high prevalence of An. species UG3 in the adult sample, it represented < 0.5% of the larval sample. The adult sample also contained two species unrepresented in 2019 larval samples (An. funestus and An. rivulorum), perhaps due to the preference of the former species for heavily vegetated water bodies 74 and tolerance of submergence 75 reducing its representation in dipped larval samples 76 . However, an independent time series of extensive catches of adult mosquitoes from the centre of the study area (n > 2400; 37,38 and Orba et al., 2021, pers. comm.) provides substantial corroboration for the anopheline community composition described from the larval surveillance presented here and previously 36 . This dataset confirms the consistent dominance of An. coustani and the notable absence of significant numbers of primary vector species in both dry and wet seasons.
The anopheline community across Barotseland, western Zambia, experiences dramatic seasonal fluctuations in environmental conditions in a highly dynamic ecosystem driven by seasonal flooding and rainfall regimes 40 . Although anophelines were encountered in a lower proportion of water bodies after the 2019 wet season than after the preceding 2018 wet season, the larval community was composed of a higher proportion of late stage larvae, suggesting increased productivity in the drying down phase of the accelerated hydrological year in 2019. This has been reported in other settings, potentially as abundant habitats become smaller and warmer: larval abundance was higher in drying streams than other habitats in Kenya and Tanzania 52 , increased with falling river levels in Sudan 77 , and both larval abundance and adult productivity increased in the early dry season in Kenya 78 . Many of the species encountered in the present study area exhibit preferences for relatively permanent water bodies and community composition was consistent during three survey periods incorporating seasonal and inter-annual disparities in ecohydrological conditions. Nonetheless, nuanced intra-annual variations due to species-specific responses to seasonal changes and the changing importance of different water body types over the hydrological year may occur over finer temporal scales 68 . As climate change is predicted to increase the frequency of meteorological extremes in Zambia 79 , it is important to monitor the response of the Anopheles community to extreme conditions beyond the typical spectrum of seasonal variations. Analysis of historical data has revealed a downward trend in Zambezi river discharge since the 1950s 80 , and Zambia is projected to experience reduced rainfall and increased temperature in the twenty-first century 81 ; a recent drought in southern Zambia reduced the abundance of An. arabiensis by an order of magnitude 69 . A significant and extreme drought in western Zambia followed the larval surveys presented in this study, and expanded entomological surveillance is advocated to characterise the effects of such events on the assemblage of anopheline species in this ecosystem.
In an area of persistent malaria transmission despite long-term indoor vector control efforts, extensive larval sampling of the anopheline community over consecutive years has revealed the numerical and spatial dominance of species widely accepted to be secondary vectors of malaria and a marked paucity of primary vectors. This finding underlines the importance of diversifying vector control approaches to counteract species whose behaviours may permit them to evade widespread use of indoor-centric interventions. The consistent dominance of these secondary vector species was strongly supported by independent surveillance of adult mosquitoes in the area www.nature.com/scientificreports/ over 5 years, and geographic distributions of candidate vector species were robust to seasonal and inter-annual variations in ecohydrological conditions. This temporal consistency of larval community structure in western Zambia suggests that this may be a setting where larval source management strategies could be highly effective, compared to other regions where the larval community fluctuates substantially over time.

Methods
Study area. The study area in Limulunga and Mongu districts in Western Province, Zambia has been characterised hydrologically 40 and five main ecological zones (Fig. 1) described in detail 36 . Briefly, it is dominated by the Barotse or upper Zambezi floodplain (ecological zone 1) which is a heterogeneous grassland-wetland mosaic which receives overbank flow from the Zambezi River and its tributaries 82,83 . The floodplain edge (zone 2) along the base of the eastern escarpment contains persistently wet seepage zones 43  by employing a purposive dipping strategy 27,30,85 to search for larvae within suitable microhabitats. Dip contents were examined in a white plastic tray after a settling period, and counted after morphological differentiation into anophelines and culicines. Up to 12 Anopheles larvae per sampling point were collected and stored individually in 95% ethanol for genetic analyses.
In addition to sampling of ecological zones described above and previously 36 , larval surveys were undertaken along Kambule stream within Mongu town to characterise the anopheline community in this peri-urban setting.
Mosquito sampling: adults. Simultaneous indoor and outdoor sampling of adult mosquitoes was undertaken during the larval sampling period using an adaptation of the BG-Sentinel trap (Biogents, Germany 86 ) optimised for anopheline mosquitoes and known as the BG-Malaria trap 87 . Each trap was suspended from a tripod and inverted with its opening 40 cm above ground level, baited with a cartridge of BG-Lure synthetic attractant (Biogents, Germany) and with CO 2 produced by fermentation of yeast (40 g) and brown sugar (500 g) in 2 l of water and released within the trap. One indoor-outdoor pair of traps was deployed at each sample household to survey both endophagic and exophagic mosquitoes; the indoor trap was located close to the foot of an occupied bed and the outdoor trap positioned in the lee of the house or nearby vegetation, adjacent to outside sitting areas and away from sources of smoke or disturbance. Trap houses were preferentially located on the outer fringes of the village due to the 'edge effect' identified in some studies [88][89][90] , downwind of the village centre as female anophelines are postulated to fly upwind in search of human hosts 91 .
Traps were deployed in late afternoon and ran until the battery and fermentation mixture were disconnected the subsequent morning. Each trap catch bag was closed and labelled, and mosquitoes killed by freezing. Each catch was subsequently examined under a dissecting microscope and screened morphologically by genus, sex and bloodfeeding status. All anophelines were retained and stored individually in 95% ethanol.
Adult mosquito data derived from previous studies. Adult mosquitoes were sampled in several districts across Zambia between 2012 and 2017 by Orba et al. 37 and Wastika et al. 38  DNA barcoding species identification. Genetic analysis of larval and adult samples by mtDNA COI and nuclear ITS2 sequencing followed Cross et al. 36 . The species identity of each specimen was inferred from BLASTn search results yielding COI and/or ITS2 sequence similarity of ≥ 95% to sequences in GenBank (National Center for Biotechnology Information). Adoption of this standard sequence similarity threshold was vital to ensure comparability with previous results 36 and those from other studies 22,23 . Species identities were further assessed by position of clustering within a maximum likelihood inferred phylogeny of the COI dataset and published confirmed species identity reference sequences 22 www.nature.com/scientificreports/ calculated using 100 bootstrap replicates. Additionally, adult female anophelines were screened for the presence of P. falciparum sporozoites by PCR amplification using primers for the protozoan 93 .
Data handling and statistical analyses. Field and molecular data were compiled in a database (MS Access) and linked to geographic locations in QGIS v. 3.18.3-Zürich 94 . Exploratory analyses were undertaken by plotting data points from individual dips together with aggregated summary measures (median values per year/ecological zone) using PlotsOfData 95 . The proportion of the genetically-identified subsample of anopheline larvae composed of each species individually was applied to the total anopheline count at each transect point to obtain the estimated total of each species at each transect point. Statistical analyses were undertaken in SPSS 96 ; comparisons of field statistics (e.g. total anopheline larvae per transect point) and sequence-derived species data (e.g. estimated An. arabiensis per transect point) were made between consecutive wet seasons (2018 36 and 2019) using Odds Ratios 97 and Pearson's Chi-squared test. Sampling in 2017 36 is referenced as the dry season comparator. Within-year comparisons between ecological zones were made using non-parametric statistics: independent samples median (ISM) test and Kruskal-Wallis H, both with stepwise step-down (SSD) post-hoc comparisons and adjusted (adj) p value for multiple comparisons.
Ethical approval. An ethical approval waiver was provided by the University of Zambia's Biomedical Research Ethics Committee (Ref 018-08-17) as the research did not involve human subjects. The Barotseland Royal Establishment granted their approval for entomological surveys to be conducted in and around villages in the study area. District Health Office staff accompanied the field survey team; at the beginning of each day's fieldwork, the survey team checked in with the nearest health facility and sought permission from village chiefs to undertake fieldwork following introductory discussions. After full explanation of the adult trapping procedure in an appropriate language (usually siLozi), written informed consent was obtained from householders who volunteered to participate, and LLINs issued where absent from indoor trapping houses. All methods were performed in accordance with relevant guidelines and regulations.

Data availability
The COI and ITS2 sequence datasets generated by and analysed in the current study are available in the NCBI GenBank nucleotide archive with accession numbers OL619678-OL619792 and OL583776-OL583807 for larval and adult COI sequences, respectively and OL621257-OL621744 and OL621788-OL621838 for larval and adult ITS2 sequences. A full maximum likelihood phylogenetic tree of COI haplotypes from this study and from Orba et al. (2021, pers. comm.) is provided in Supplementary Information (Fig. S1). Anopheles count data and six COI sequences from the latter source are provided in Supplementary Table S1 and Supplementary Dataset S1.