A differential expression of pyrethroid resistance genes in the malaria vector Anopheles funestus across Uganda is associated with patterns of gene flow

Background Insecticide resistance is challenging the effectiveness of insecticide-based control interventions to reduce malaria burden in Africa. Understanding the molecular basis of insecticides resistance and patterns of gene flow in major malaria vectors such as Anopheles funestus are important steps for designing effective resistance management strategies. Here, we investigated the association between patterns of genetic structure and expression profiles of genes involved in the pyrethroid resistance in An. funestus across Uganda and neighboring Kenya. Methods Blood-fed mosquitoes An. funestus were collected across the four localities in Uganda and neighboring Kenya. A Microarray-based genome-wide transcription analysis was performed to identify the set of genes associated with permethrin resistance. 17 microsatellites markers were genotyped and used to establish patterns of genetic differentiation. Results Microarray-based genome-wide transcription profiling of pyrethroid resistance in four locations across Uganda (Arua, Bulambuli, Lira, and Tororo) and Kenya (Kisumu) revealed that resistance was mainly driven by metabolic resistance. The most commonly up-regulated genes in pyrethroid resistance mosquitoes include cytochrome P450s (CYP9K1, CYP6M7, CYP4H18, CYP4H17, CYP4C36). However, expression levels of key genes vary geographically such as the P450 CYP6M7 [Fold-change (FC) = 115.8 (Arua) vs 24.05 (Tororo) and 16.9 (Kisumu)]. In addition, several genes from other families were also over-expressed including Glutathione S-transferases (GSTs), carboxylesterases, trypsin, glycogenin, and nucleotide binding protein which probably contribute to insecticide resistance across Uganda and Kenya. Genotyping of 17 microsatellite loci in the five locations provided evidence that a geographical shift in the resistance mechanisms could be associated with patterns of population structure throughout East Africa. Genetic and population structure analyses indicated significant genetic differentiation between Arua and other localities (FST>0.03) and revealed a barrier to gene flow between Arua and other areas, possibly associated with Rift Valley. Conclusion The correlation between patterns of genetic structure and variation in gene expression could be used to inform future interventions especially as new insecticides are gradually introduced.

including Glutathione S-transferases (GSTs), carboxylesterases, trypsin, glycogenin, and nucleotide binding protein which probably contribute to insecticide resistance across Uganda and Kenya. Genotyping of 17 microsatellite loci in the five locations provided evidence that a geographical shift in the resistance mechanisms could be associated with patterns of population structure throughout East Africa. Genetic and population structure analyses indicated significant genetic differentiation between Arua and other localities (F ST >0.03) and revealed a barrier to gene flow between Arua and other areas, possibly associated with Rift Valley.

Conclusion
The correlation between patterns of genetic structure and variation in gene expression could be used to inform future interventions especially as new insecticides are gradually introduced.

Background
Malaria remains one of the main causes of morbidity and mortality in Sub-Saharan Africa, predominantly in children under 5 years and pregnant mothers [1]. The transmission of these malaria-causing parasites to human is caused by Anopheles mosquitoes of which four species (An. gambiae, An. coluzzii, An. funestus, An. arabiensis,) have been identified as the major malaria vectors in Africa. In Uganda, the main malaria vectors are Anopheles funestus, An. gambiae sensu strict and An. arabiensis [2].
Malaria control in Uganda relies heavily on vector control using long-lasting insecticide nets (LLINs) and indoor residual spraying (IRS) mostly in regions of seasonal transmission [1]. The success of such interventions requires a good knowledge of vector populations particularly their susceptibility status to the main insecticides used for such control programs. The two major causes of insecticide resistance are alterations in the target sites and an increase in the rate of insecticide metabolism [3][4][5]. In Uganda, the previous investigation of pyrethroid resistance has revealed the absence of knockdown resistance (kdr) target-site mutation in An. funestus [4,6,7]. The underlining molecular basis of resistance to the insecticide for this vector in Uganda has been associated with cytochrome P450 over-expression in the eastern part (Tororo) [3] but it remains to establish if the same mechanisms are present country-wide and in neighboring Kenya [3]. It is important to know whether the resistance front is unique, or gene flow is uniform across the region. Previous efforts to characterize the mechanisms of pyrethroid resistance in An. funestus has revealed that resistance is mainly driven by metabolic resistance [8][9][10]. Cytochrome P450s are known to be a primary enzyme family conferring resistance to pyrethroids in malaria vectors [11]. Molecular studies conducted in Malawi and Mozambique have revealed that the duplicated P450 genes, CYP6P9a, and CYP6P9b are the main genes driving pyrethroid resistance in this region [10,[12][13][14][15]. In addition, the studies performed in Malawi have revealed a similarly minor role of CYP6P9a and CYP6P9b [7]. Recently, molecular studies conducted on An. funestus in southern Africa (Malawi), East (Uganda), and West Africa (Benin) have revealed that the duplicated P450 genes (CYP6P9a and CYP6P9b), which were highly overexpressed in southern Africa, were not the most upregulated in other regions, where other genes, including GSTe2 in West (Benin) and CYP9K1 in East (Uganda) [3], were overexpressed. This variation of expression profiles in Africa suggests that the molecular basis of pyrethroid resistance might vary across Africa and within national populations of An. funestus. However, the molecular basis of pyrethroid metabolic resistance in An. funestus across Uganda remains poorly characterized [4,6].
It also remains unknown whether the same control strategy could efficiently control An. funestus populations throughout Uganda and the neighboring countries. This is further the case in the context of ecological landscape changes such as the Rift Valley that spans East Africa and previously shown to restrict gene flow in An. gambiae [16,17]. Indeed, based on microsatellite markers, the magnitude of genetic differentiation (F ST ) between populations on opposite sides of the continent (~6000 km apart) was~0.03, while the corresponding value between populations on either side of the Rift Valley (~400-500 km) was~0.1 [17][18][19][20]. The genetic structure of An. funestus across Uganda remains poorly characterized in the context of the spread of insecticide resistance and impact of the Rift valley although a recent study highlighted a homogeneity between populations from Central and North Uganda [21] Assessing how mechanisms of pyrethroid resistance vary countrywide and whether it is linked to the genetic structure vector populations is an important step in designing nationwide vector control strategies. Furthermore, screening more populations could detect new genes driving such pyrethroid resistance. Identifying the full set of genes involved in resistance will help decipher the molecular basis of resistance and potentially identify resistance markers which can be used in the design of DNA-based molecular diagnostic tools for quick detection and tracking resistance in the field as recently done for P450-mediated metabolic resistance in southern African populations of An. funestus [12,13].
In this present study, using microarray genome-wide transcription analysis, we characterized the molecular basis of pyrethroid resistance in An. funestus in Uganda and Kenya. We also provide evidence, using 17 microsatellite markers, that the genetic structure of the An. funestus in both countries varied and correlates with changes in gene expression.

Study sites and samples
An. funestus mosquitoes were collected between 06.00 AM and 12.00 PM, from four districts in Uganda; Arua (Ar) in North West (3˚1 0 N, 30˚54 0 E), Bulambuli (Bl) in North-East (1˚10 0 N, 34˚23 0 E), Lira (Lr) in North Central (2˚14 0 N, 32˚54 0 E), and Tororo (Tr) in East (0˚41 0 N, 341 0 0 E). A similar collection was made from Kisumu-Siaya district (0˚05 0 S, 34˚15 0 E) in West Kenya. Mosquito collections were carried out between December 2011 and June 2012: between the end of dry season and beginning of the rainy season, with temperatures and relative humidity ranging from 26˚C to 29˚C and 66% to 77% respectively. Indoor resting females were collected randomly in different locations using electric aspirators as described previously [4]. No specific permissions were required for these locations/activities and these field collections did not involve endangered or protected species. The blood-fed F 0 adults were morphologically identified as belonging to the An. funestus group according to the key of Gilles and Coetzee (1987) [22]. They were left to become gravid and forced to lay eggs using the forcedegg laying method (Morgan et al 2010). A PCR assay was performed using the protocol of Koekemoer [23] to confirm that collected F 0 adults were An. funestus s.s. [4]

Resistance profile of different populations districts
The resistance patterns of the five populations districts in Uganda and Kenya to 0.75% permethrin insecticides was determined as described previously by Mulamba et al., [4] following the World Health Organization (WHO) protocol [24]. The Arua, Tororo, and Kisumu population of An. funestus are highly resistant to permethrin (27% mortality, 33% mortality, and 20% mortality, respectively for Arua, Tororo, and Kisumu after 1 hr exposure) [4]. The Bulambuli and Lira population of An. funestus are also resistant to permethrin (49% mortality and 51% mortality respectively for Bulambuli and Lira, after 1hr exposure) and the final mortality was assessed after 24 h [4].

Detection of pyrethroid resistance genes using microarrays and qRT-PCR
Genes associated with pyrethroid resistance in the five locations were detected using the 8 X 60K Agilent microarray chip customarily designed for An. funestus as previously described [10]. This chip designed through the Array program (Agilent) (A-MEXP-2374) and previously described by Riveron et al [10], is made of sets of 15,527 transcripts generated from de novo transcriptome analysis [25], 8,540 Expressed Sequence Tags (ESTs) [26]. It also includes a set of 2850 An. funestus cDNAs from GenBank, a set of P450 genes from the rp1 and rp2 QTL BAC sequence [14,27]; and 13,000 transcripts of the complete An. gambiae genome. Total RNA was extracted from three pools of 10 female mosquitoes per phenotype notably in Lira for which mosquitoes included: Control (unexposed to insecticide, C), Resistant [alive R) after exposure to 0.75% permethrin], and Susceptible (FANG susceptible colony; S). In other locations, the resistant (R) sample was included. The extraction was performed using the Picopure RNA isolation Kit (Arcturus, Applied Biosystems, Carlsbad, CA, USA). The RNA extraction was performed as previously described by Riveron et al., 2017 [3]. The quantity and quality of extracted RNA were assessed using a NanoDrop ND1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and a Bioanalyzer (Agilent, Santa Clara, CA, USA), respectively. The Agilent Quick Amp Labeling Kit (two-color) was used to amplify the complementary RNA (cRNA) from each sample according to the manufacturer's protocol and as described previously [3]. Resistant cRNAs (R) were labeled with cy5 dye, whereas susceptible cRNAs from the strain FANG (S) were labeled with the cy3 dye. The quantity and quality of all cRNAs were assessed using the NanoDrop and Bioanalyzer before labeling. Labeled cRNAs were hybridized to respective arrays for 17 h at 65˚C following the manufacturer's protocol. For each location and comparison, five hybridizations were performed by swapping the biological replicates.
Agilent GeneSpring GX 13.0 software was used to analyze the microarray data. The differentially expressed genes were identified using a threshold of twofold-change (FC) and a statistical significance of P<0.05 with Benjamin-Hochberg correction for multiple testing. The microarray data from this study is deposted in the Array Express under accession no. E-MTAB-9529. Finally, quantitative reverse transcription PCR (qRT-PCR) assays were performed to validate microarray results for the key candidate genes; 1 μg of RNA from key resistance genes, comparing the permethrin-resistant mosquitoes to FANG susceptible (R-S) mosquitoes, was used as a template for complementary (cDNA) synthesis using the superscript III (Invitrogen, Carlsbad, CA, USA) following the manufacturer's guide. The qRT-PCR was carried out as previously described [10,28] with the relative expression level and fold change of each target gene calculated according to the 2 −ΔΔCT method after normalization with the housekeeping genes, the ribosomal protein S7 (RSP7; AFUN007153) and actin5C (AFUN006819) [29].

Genetic population structure of Anopheles funestus in Uganda and Kenya
Randomly field-collected F 0 females confirmed as An. funestus s.s were genotypes at seventeen microsatellites loci genome-widely distributed [30,31].
The mosquitoes from Uganda and Kenya samples (N = 43 in Arua, N = 48 in Bulambuli, N = 43 in Lira, N = 47 in Tororo, N = 26 in Kisumu) were genotyped as previously described by [7]. Briefly, genomic DNA extracted from F0 mosquitoes were used to amplify the 17 microsatellite loci (both di-and tri-nucleotide repeats) using 1.5 μl of reaction Buffer, 0.2 μl of dNTP mix (25 mmol), 0.325 μl of both the forward (included a 19bp tag) and reverse primers, 0.2 μl of Hot Start Taq (Qiagen Inc.), 1 μl of MgCl 2 and 1μl of genomic DNA (15ng/ul). S1 We applied a Bayesian model-based clustering algorithm to infer population structure and to assign individuals (probabilistically) to clusters without a priori knowledge of population units and limits. This approach, implemented in STRUCTURE v 2. The correlation between genetic and geographical distances was assessed by the regression of F ST / (1-F ST ) on the logarithm (ln) of geographical distance [42], and tested using the Mantel test available in GENEPOP.
Kruskal-Wallis test was used to compare the average number of alleles and the average proportions of heterozygosity between populations using GraphPad Prism 5. The Bonferroni correction procedure [43] was applied to evaluate significance when multiple tests were performed.

Transcription profiling of the pyrethroid resistant population of Lira
The triangular hybridization performed in Lira between mosquitoes resistant to permethrin (R), unexposed to insecticide (C), and the FANG susceptible laboratory strain (S) allows a comparative analysis of transcription profile in this location. A high number of the probes were significantly differentially expressed (p<0.05) for the R-S (9263), R-C (4132 with foldchange of 1.5), and the C-S (4128) (Fig 1). However, a Venn-diagram analysis revealed that only 182 probes were commonly differentially expressed in all three groups (Fig 1).
Genes commonly overexpressed in R-C, C-S, and R-S strains. Among these, the probes with the highest over-expression in the R-S comparison (FC79.8) and other included a hypothetical protein ortholog of the AGAP000603 in An. gambiae with an unknown function. The nucleotide-binding protein 2 (nbp2) (Afun008887) gene was also consistently overexpressed in all three comparisons although with the highest fold change seen in R-S (FC: 17.15) ( Table 1). This gene was also significantly overexpressed in R-C and C-S with a much lower FC value (FC = 1.7; FC = 7.2 respectively for R-C and C-S) ( Table 1). The sg2 salivary protein (EE589890.1) corresponding to AFUN016226 was also commonly over-expressed. The sulfotransferase gene (Afun013871) and the aldehyde oxidase (AGAP006225) were other detoxification genes commonly expressed in R-C, R-S, and C-S with a similar fold change (Table 1).
Genes commonly overexpressed in R-S and C-S strains. Several detoxification genes or resistance-related genes were commonly and significantly overexpressed in R-S and C-S strains. Among the most overexpressed genes commonly observed in R-S and C-S were proteases such as a trypsin-related protease (Afun008293), which was the top upregulated with FC 133.61 in R-S and 86.53 in C-S. Several detoxification genes were commonly upregulated in both strains, with cytochrome P450s being the most; notably, CYP9K1 (three probes) overexpressed with FC 13.68 in R-S and 31.80 in C-S. Two other genes CYP9J11 and CYP9J3 were also up-regulated with FC 5.44 in R-S and 4.38 in C-S for CYP9J11 and FC 3.72 in R-S and 2.72 in C-S for CYP9J3. CYP6M7 was also overexpressed with FC 3.90 in R-S and 2.94 in C-S   Table 1). Glutathione-S-transferases (GSTs) were also significantly overexpressed in pyrethroid-resistant mosquitoes from Lira compared to the susceptible FANG strain, notably GSTe3 (Afun008354) (FC 3.8 and 3.6, in R-S and C-S respectively) ( Table 1).

Comparison of expression profiles between Arua, Bulambuli, Lira, Tororo in Uganda, and Kisumu in Kenya
To detect the set of genes associated with permethrin resistance across Uganda and neighboring Kenya, the same 8X60k microarray chip was used to compare mosquitoes alive after permethrin exposure and compared to the full susceptible laboratory strain (R-S). The number of probes that were differentially expressed (>2-fold change, FC) between R and S mosquitoes for each locality and between them is indicated in Fig 2 (P < 0.01). Overall, 3553 probes were differentially expressed in the Arua population, 5903 in the Bulambuli population, 9263 in the Lira population, 7346 in the Tororo population and 5598 in Kisumu. When comparing the four Uganda populations, a total of 1420 probes were commonly differentially expressed (Fig  2A) whereas 1098 probes were commonly differential expressed when Kisumu was compared to Lira, Tororo and Bulambuli (Fig 2B).
Genes commonly up-regulated in the five locations (Arua, Bulambuli, Lira, Tororo, and Kisumu). Among the detoxification genes, several cytochromes P450s were most commonly

PLOS ONE
overexpressed in all the five locations. Most of these P450 genes showed a similar level of expression in all the five localities and included CYP6Z1 (two probes), CYP9K1, and CYP6AG1. However, the overexpression levels of some candidate genes show significant geographical variation between locations. In Arua, CYP6M7 (FC 115.88) has the highest overexpression levels than other locations. Another P450, CYP4H18, although commonly expressed in all five localities, was significantly present in Arua than other localities, suggesting a bigger role for this gene in Arua and supporting a differential expression in Arua compared to other Uganda locations. This shift in expression pattern was also observed for other genes commonly overexpressed in all five locations. These include Afun000500 (glycogen) and Afun007894 (Trypsin) with the highest overexpression in Arua (FC 39.71 and FC 4.96) ( Table 2).

Genes commonly up-regulated only in all locations (Arua, Bulambuli, Lira, and Kisumu) but not in Tororo.
Only CYP6Y2 (two probes) gene was commonly over-expressed in four localities with a similar expression level ( Table 2).
Genes commonly up-regulated in all locations (Arua, Bulambuli, Tororo, and Kisumu) but not in Lira. CYP4H25 was the most overexpressed gene with the highest FC value in Arua (FC, 9.22) followed by Tororo, Bulambuli, and Kisumu with FC values of 5.34, 3.04 and 2.74, respectively ( Table 2).
Genes commonly up-regulated in all locations (Bulambuli, Lira, Tororo, Kisumu) but not Arua). Among the commonly up-regulated detoxification genes, cytochrome P450s were the most predominant with eight genes over-expressed; while only a single glutathiones-transferase GST gene (GSTD3) was over-expressed. Of the P450s, CYP6M4 was the most overexpressed gene with the highest FC value in Tororo (FC 9.17) followed by Lira, Kisumu, and Bulambuli with FC values of 5.25, 4.99 and 3.73 respectively. Most of these P450 genes

PLOS ONE
showed a similar level of expression in all the four localities and included CYP304B1, CYP6M1, CYP4K2, CYP4H16, CYP9J3, and CYP6M3 ( Table 2). The unique glutathione-stransferase (GSTD3) gene commonly up-regulated in the four locations was overexpressed with the highest FC value in Kisumu (FC 11.06) followed by Bulambuli, Tororo, and Lira with FC values of 5.56, 5.07 and 3.80 respectively ( Table 2). Genes common only in two localities. Analysis of the list of genes commonly overexpressed in only two localities revealed that, for Arua and Bulambuli, only CYP325D1 was observed and showed high expression in Arua. For those overexpressed only in Arua and Lira, the P450 CYP6P1 and CYP6P9b were detected, although with higher expression in Arua (FC, 3.47 and 4.06 respectively) than in Lira (FC, 2.87 and 3.02 respectively). The CYP306A1 was also upregulated and showed a similar level of expression in both localities. The list of genes overexpressed only in Arua and Tororo was dominated by the CYP4C36 and CYP307A1 with higher overexpression in Arua for both genes (e.g., FC 19.17

Analysis of genetic population of Anopheles funestus in different locations
The significant differences in the gene expression profiles observed between populations of An. funestus in Uganda (notably between Arua and others) and Kenya could suggest the presence of barriers to gene flow that are affecting the spread of resistant genes. However, knowledge of the population genetic structure of An. funestus in Uganda and Kenya is limited.
Genetic diversity and Hardy-Weinberg equilibrium. Genotypes of An. funestus females were scored at 17 microsatellites loci. These microsatellite loci were highly polymorphic with a number of distinct alleles per locus ranging from 8 (AFUB12) to 20 (FUNO) for the combined five populations ( Table 3). The average number of alleles per locus ranged from 7.1 (Kisumu) to 8.7 (Balambuli) and was not significantly different among populations (P = 0.92). However, Kisumu showed the lowest number of alleles for many loci with the minimum of 3 alleles observed for AFUB12 in this location. Mean observed heterozygosity across all loci ranged from 0.689 (Tororo) to 0.717 (Arua) and was not significantly different among populations (P = 0.85).
A significant heterozygosity deficit was observed in 29 out of 85 tests across the markers after Bonferroni correction at (P< 0.01) across the markers as shown by F IS estimates (Table 3). Some markers such as AFUB12, FUNQ, and FUNR exhibited such heterozygosity deficit in 4 out of 5 locations suggesting that such deviation could be marker related. Kisumu had the lowest number of deficit (3/17) whereas Bulambuli had the highest (10/17) (Table 3).
When the pooled samples were analyzed as a single population, significant deviation from Hardy-Weinberg equilibrium (F IS = 0.135; P<0.01) was observed within each population studied due to significant heterozygote deficiency (Table 3). No linkage disequilibrium was observed in any pair of loci after correction by the Bonferroni procedure (P>0.05) suggesting genetic independence between loci.
Genetic differentiation. The levels of genetic differentiation between pairs of populations were estimated by F ST values. Table 4 shows F ST estimates for all pairwise populations compared. The values of F ST between pairwise population comparisons for all loci ranged from 0 (Bulambuli-Lira) to 0.037 (Arua-Kisumu). The highest significant F ST estimates were obtained between Arua and other localities (Bulambuli, Lira, Tororo, and Kisumu) suggesting the presence of barriers to gene flow between these Arua and these locations (Bulambuli, Lira, Tororo, and Kisumu) ( Table 4). Analysis of patterns of genetic differentiation at individual loci revealed higher F ST estimates between Arua and other populations at some loci such as FUNO (0.27-0.32) and AFND6 (0.04-0.091). These two loci are located on the 2R chromosome at the vicinity of a major QTL (rp1; resistance to pyrethroid 1) previously associated with resistance to pyrethroids and close to a cluster of cytochromes P450s. This high F ST estimates could suggest a difference in selection pressure between Arua and other locations in this genomic region in association with a difference in resistance mechanism supported by the variation in gene expression patterns from the microarray study. In contrast, the four locations of Bulambuli, Lira, Tororo, and Kisumu showed a very low and nonsignificant F ST pairwise estimate, suggesting a high level of gene flow between these populations ( Table 4). The segregation in the Arua population is similar to changes in gene expression profiles of resistance genes in this location. Furthermore, to understand the possible role of geographical distance in generating the genetic distance between Arua and other populations, the Mantel test was performed. The test revealed no significant correlation between the pairwise F ST / (1-F ST ) against the natural logarithm of pairwise geographical distance (R 2 = 0.19, P = 0.806), suggesting that the population genetic structure of An. funestus in Uganda and Kenya did not conform to isolation by  distance model (Fig 4) and suggesting that the difference in Arua may be caused by other factors. Finally, to confirm the estimated genetic differences inferred based on F ST , we performed Bayesian predictions of population structure. The population size was K = 3, predicted using the Evanno method (Fig 5). In agreement with the F ST estimates, the Bulambuli, Lira, Tororo, and Kisumu samples from Uganda and Kenya shared a similar pattern whereas Arua forms a distinct cluster (Fig 5).

Discussion
Insecticide resistance in An. funestus mosquitoes is spreading throughout Africa, threatening the success of malaria control methods. In this study, we characterized the mechanisms driving insecticide resistance in An. funestus population from Uganda and neighboring border town in Kenya. We further investigated the population genetic structure in the malaria vector An. funestus to help predict the pattern of spread of resistance and improve the design of insecticide resistance management strategies.
An. funestus populations in Uganda and neighboring Kenya are resistant to permethrin [4]. Microarray analysis identified several transcript coding for detoxification enzymes (P450s, GSTs, Carboxylesterase, Trypsin, Glycogenin, and nucleotide-binding protein) upregulated in the five localities. Among the commonly up-regulated detoxification genes, cytochrome P450s were the most predominant with several genes over-expressed as previously shown in several  [3,12] suggests that the origin of resistance is not the same across the countries, suggesting that independent selection events of resistance to pyrethroids have occurred in various populations [3].  Arua. Surprisingly, one of the duplicated P450 genes CYP6P9a, which has been shown to play a main role in pyrethroid resistance in southern populations of An. funestus [10,14] was not overexpressed in all the localities examined. The complete absence of overexpression of CYP6P9a indicates that the resistance mechanism in Uganda and neighboring Kenya is different from that observed in Malawi [3].
The reasons for the potential shift in the expression levels of the genes remain unknown but could be due to the nature of the selection that gave rise to the resistance. This shift in gene expression further highlights the genetic plasticity of natural populations of malaria vectors and their ability to adapt to various selection pressures. The finding of geographical differences in the role of key resistance genes in Arua suggests the presence of barriers to gene flow. F ST values have shown a greater genetic differentiation between Arua and other localities (Bulambuli, Lira, Tororo, and Kisumu) suggesting the presence of barriers to gene flow. In addition, An. funestus populations from Uganda and neighboring Kenya are subdivided into two distinct genetic entities: Population around the Rift Valley (Arua) and populations from other side of the valley (Bulambuli, Lira, Tororo, and Kisumu). The existence of these two genetic entities was confirmed by different genetic approaches (i.e. Structure and genetic differentiation). This result confirms the fact that another factor explains the population genetic structure of Anopheles funestus in these localities. This observation suggests the impact of the Rift valley as a barrier to gene flow between populations of An. funestus in this region. Our results corroborate the previous microsatellite study in An. gambiae in Uganda, where Rift valley shows a great barrier to gene flow [16,17]. However, the greater estimates of F ST at two loci located at the vicinity of a known pyrethroid resistance genomic region (rp1) [14], suggests that the genetic differentiation observed between Arua and the other locations in Uganda could also be influenced by a difference in local selection as the FUNO and AFND6 loci have previously be shown to be under selection due to location of the cluster of pyrethroid resistance CYP6 P450 genes in this region [7,47]. It will be useful in the future to investigate the selection sources of pyrethroid resistance in Arua, either from different agricultural practices or local vector control interventions, to establish if these can also account for the difference observed between this population and others.

PLOS ONE
The genetic test of isolation by distance suggested no significant correlation between genetic diversity and geographical distance confirming that this difference could be due to ecological or geographical factors [50]. The cause of this barrier could also be associated with the chromosomal differentiation because the difference of frequency for some inversions such as 3La was reported to impact the genetic structure of An. funestus [51].
When Arua was excluded, our results showed low levels of genetic differentiation between An. funestus populations (Bulambuli, Lira, Tororo, and Kisumu), suggesting that overall there is a high level of gene flow between populations of this species across most of Uganda showing a genes of interest such as insecticide resistance genes or future gene drive constructs, could spread quickly among populations of this vector. A previous study on the genetic structure of

Conclusion
The genetic structure and variation in gene expression could be used to make an informed decision in future interventions especially as new insecticides are needed to control malaria across these countries. In addition, the finding of differences in the molecular basis of resistance of permethrin within a given country means that national resistance management strategies without characterization of underlying resistance mechanisms from localities may be flawed. The similarity of resistance profiles in Bulambuli, Lira, Tororo, and Kisumu, suggests that the same resistance management strategy could be implemented across these localities but might need to be different in Arua. These results have added to our understanding of the dynamics of vector species and will be valuable for planning effective vector control activities based on population genetic structure and resistance genes in Uganda and neighboring Kenya.
Supporting information S1