Signatures of selection and mechanisms of insecticide resistance in Ugandan Anopheles funestus: Insights from embedding translational genomics into the LLINEUP cluster randomised trial

Abstract

In response to the emerging threat of insecticide resistance in malaria vectors, insecticides are being repurposed for vector control or developed de novo. Good stewardship of these finite new resources is essential if disease control programmes are to remain effective. This is dependent on timely data to help guide evidence-based decision-making for National Malaria Control programmes (NMCPs). By embedding genomics into cluster randomized control trials (cRCTs), we can perform surveillance and early detection of insecticide resistance variants to new and repurposed chemicals in natural field conditions, supporting effective stewardship. The LLIN Evaluation Uganda Project (LLINEUP) trial evaluated the efficacy of pyrethroid-piperonyl butoxide (PBO) and pyrethroids-only long-lasting insecticidal nets (LLINs). It was conducted in Uganda between 2017-2020 and was the largest cRCT to date, covering 40% of the country in 104 health sub-districts. We embedded genomic surveillance within LLINEUP to detect and track insecticide resistance variants. At baseline and throughout the trial, we sampled Anopheles mosquitoes with Prokopack aspirators and performed Illumina whole-genome sequencing. We show that An. funestus populations were relatively unaffected by the interventions, compared to An. gambiae s.l., which were markedly reduced six months following LLIN deployment. Standard approaches for describing genetic diversity and population structure e.g. fixation index (F_ST ), Principal Component Analysis (PCA) and Neighbour-Joining (NJ) trees, were consistent with the density observations and suggestive of a single large An. funestus population in Uganda with little genetic differentiation. Genome-wide selection scans revealed strong signals of selection at the Resistance to pyrethroid-1 (RP1) locus and Cyp9k1, both loci previously implicated in pyrethroid resistance. We report two additional loci, eye diacylglycerol kinase (Dgk) (≅ 13.5Mb on the X chromosome) and O-mannosyl-transferase (TMC-like) (≅ 67.9Mb on 3RL) that showed signals of selection. Known DDT and permethrin resistance-associated variants at the Gste2 locus, L119F and L119V, were also identified. Over the trial period, changes in haplotype frequencies were observed in regions under selection, with more pronounced shifts in the PBO arm. Notably, there were significant reductions in the frequencies of swept haplotypes (measured by delta (Δ) H12) in the Dgk and Cyp6p9a regions, while significant increases in haplotype frequency were observed at Gste2 and Cyp9k1 loci. Our findings reveal the differential impact of the trial on An. gambiae s.l. and An. funestus densities and the differing responses of An. funestus populations to pyrethroid and pyrethroid-PBO selection pressure. These insights underscore the potential value of tailored, species- and region-specific vector control strategies, supported by regional genetic surveillance, to better control insecticide resistance evolution and spread. By embedding genomic surveillance in cRCTs we can facilitate the discovery of putative resistance variants and can provide evidence of their impact on vector control tool efficacy; both of crucial importance to evidence-based deployment of vector control tools by NMCPs.

Fig 1:. Map of Uganda showing LLINEUP sample collection sites

A: (pink) and the rest of the districts in the country (grey area). Box-whisker plots of numbers of An. gambiae s.l. (B) and An. funestus (C) per household during the trial. The Y-axis is truncated to 75 to improve visual clarity. We performed a Negative Binomial Model ($glmmTMB(species\sim round+LLINtyp{e}^{+}(1\mid household),family=nbion2)$ to determine the number of vectors per household.

Fig 2:. Population structure and genetic diversity.

A) Neighbour-joining trees (NJT) showing a lack of population structure in An. funestus population from Eastern and Western Uganda (left) and during the different trial rounds (right). NJTs, B) Nucleotide diversity (π) and C) Tajima’s D were calculated using the genomic region 2R:60-80Mbps to avoid polymorphic chromosomal inversions.

Fig 3.. Selection analyses

A) H12 genome-wide selection scans of the 2RL, 3RL, and X chromosomes. Scans were performed using 500 SNP windows; peaks are labeled by loci likely to be driving resistance. Increase haplotype frequency; on the 2RL chromosome, the signal is centered ~8,600,000Mb corresponding to RP1 locus and another signal ~76,400,000 corresponding to the Glutathione S-transferase epsilon 2 gene; on 3RL chromosome, the signal centered ~67,900,000 corresponding to O-Mannosyl-transferase; on the X chromosome, signal at ~8,400,000Mb corresponding to Cyp9k1 gene and another signal ~13,004,000 corresponding to the Diacylglycerol Kinase gene. B) ΔH12 analysis of the 2RL, 3RL, and X chromosomes. Change in haplotype frequency (ΔH12) in the Cyp6p9a, and Cyp9k1 (2 and X chromosomes), a decrease in haplotype frequencies in the Dgk region on the X chromosome. The colored dots show windows that surpassed permutation-based significance threshold (purple=marginal; green=significant).

Fig 4:. Hierarchical clustering

dendrogram of haplotypes over the A) TMC (3RL: 67,925,77167,931,656 and B) Dgk (X: 13,590,000-13,690,044) genes. Dendrogram leaves are labeled by country of mosquito population origin.