Local selection in the presence of high levels of gene flow: Evidence of heterogeneous insecticide selection pressure across Ugandan Culex quinquefasciatus populations.

BACKGROUND: Culex quinquefasciatus collected in Uganda, where no vector control interventions directly targeting this species have been conducted, was used as a model to determine if it is possible to detect heterogeneities in selection pressure driven by insecticide application targeting other insect species. METHODOLOGY/PRINCIPAL
FINDINGS: Population genetic structure was assessed through microsatellite analysis, and the impact of insecticide pressure by genotyping two target-site mutations, Vgsc-1014F of the voltage-gated sodium channel target of pyrethroid and DDT insecticides, and Ace1-119S of the acetylcholinesterase gene, target of carbamate and organophosphate insecticides. No significant differences in genetic diversity were observed among populations by microsatellite markers with HE ranging from 0.597 to 0.612 and low, but significant, genetic differentiation among populations (FST = 0.019, P = 0.001). By contrast, the insecticide-resistance markers display heterogeneous allelic distributions with significant differences detected between Central Ugandan (urban) populations relative to Eastern and Southwestern (rural) populations. In the central region, a frequency of 62% for Vgsc-1014F, and 32% for the Ace1-119S resistant allele were observed. Conversely, in both Eastern and Southwestern regions the Vgsc-1014F alleles were close to fixation, whilst Ace1-119S allele frequency was 12% (although frequencies may be underestimated due to copy number variation at both loci).
CONCLUSIONS/SIGNIFICANCE: Taken together, the microsatellite and both insecticide resistance target-site markers provide evidence that in the face of intense gene flow among populations, disjunction in resistance frequencies arise due to intense local selection pressures despite an absence of insecticidal control interventions targeting Culex.

Introduction

Contemporary evolution of insecticide resistance in mosquitoes, which threatens the effectiveness of insecticide-based approaches to reduce the burden of vector-borne disease in endemic regions, results from a complex interaction of biological traits (e.g. genetic mechanisms of resistance) and operational implementation of control programs (e.g. choice of insecticide class and selection strength see [1, 2]).

Target-site mutations in either the acetylcholinesterase gene (Ace-1) or the voltage-gated sodium channel (Vgsc) gene, as well as metabolic resistance mediated by upregulation of, or mutations in, detoxification gene families such as the cytochrome P450 (P450s), esterases and glutathione S-transferases (GSTs) are the most well studied genetic mechanisms underpinning resistance [3-5]. Pyrethroids or DDT bind to the VGSC neuronal ion channel resulting in the channel pore being unable to close, so causing repetitive nerve firing, paralysis and death [6]. A number of resistance mutations have been identified in this gene in a range of insect species with an L to F change at codon 1014 frequently identified [7]. By contrast, carbamate and organophosphate insecticides target the enzyme acetylcholinesterase (AChE) at cholinergic synapses, blocking transmission of nerve impulses by inhibition of AChE. A G to S mutation at codon 119 (Torpedo numbering system) of the Ace-1 gene results in resistance and is widespread in organophosphate and carbamate selected populations [8, 9]. In Culex and other vectors and pest insects the intense use of insecticides has been associated with sharp increases of Vgsc-1014F –also called kdr (knock-down resistance) [10-12] and an increased frequency of the Ace1-119S resistant allele.

Despite the recognized role of such target-site mutations in the adaptability of mosquitoes under insecticide selection, many studies have also shown that such resistant-associated mutations can have a range of deleterious effects and fitness disadvantages (see [13]). In Culex and other resistant insect populations where such fitness detriments are seen, the directional selection of adaptive target-site mutations has been linked with an increase in gene copy number. The most well documented example of this is for the duplication of Ace-1, which has been reported in diverse arthropod species such as Aphis gossypii, Tetranychus evansi, Anopheles gambiae as well as in Culex from a range of geographic regions [14-17]. Many studies have proposed that duplication of Ace-1 acts as a compensatory adaptive mechanism to restore the disrupted function of the nervous system driven by the resistant allele so reducing the fitness costs [18-20]. Additionally, duplication of the Vgsc gene was also recently reported in both Cx. quinquefasciatus and Aedes aegypti, indicating a likely adaptive role of gene duplication in the evolution of resistance to pyrethroids [21, 22].

In addition to the underlying genomic basis of resistance, the evolutionary pattern of insecticide resistance can also be influenced by the heterogeneity of environmental selection pressures encountered by mosquito populations. For example, mosquitoes from urban and rural settlements could experience contrasting patterns of selection imposed by direct or indirect vector control interventions, and/or the indirect impact of pesticides and environmental pollution encountered by vectors in rural settings [23]. From a vector control perspective, developing selective control approaches targeting only the primary species of interest is challenging since in most endemic regions, species with roles in vector-borne disease transmission can be co-endemic, heightening the risk of insecticide selective pressure on non-target vector species.

On the African continent, lymphatic filariasis (LF) is predominantly transmitted by Anopheles mosquitoes in rural areas, while in urban and coastal areas of East Africa Culex species are the main vector [24]. In Uganda, a country with high rates of malaria transmission and LF prevalence, Anopheles mosquitoes are incriminated as the primary vector of both diseases [25, 26]. Therefore, vector control operations using insecticide-based approaches, either through insecticide treated nets (ITN) or indoor residual spraying (IRS), chiefly target Anopheles. Whilst there is no state vector control programme targeting Culex populations across Uganda, insecticide resistance may be indirectly selected through the use of insecticides for the control of Anopheles, other vectors or pest species [23, 27].

In this study, we used collections of Cx. quinquefasciatus from Uganda as a model to determine if it is possible to infer heterogeneities in insecticide selection on non-target species through combined study of putatively selected resistant alleles (Vgsc-1014F and Ace1-119S) and neutral markers (microsatellites). To this end, microsatellite markers were developed and utilized to detect likely influence of geographic and colonization events on Cx. quinquefasciatus population diversity and structure [28, 29]. By contrast, changes in the pattern of resistant alleles (Vgsc-1014F and Ace1-119S) was used a proxy of the strength of insecticidal selection across Uganda.

Materials and methods

Ethics statement

Ethical approval was obtained from the Makerere University School of Medicine Research and Ethics Committee. This was an experimental setting for research purposes and was not linked to a malaria control program. Adult household representatives gave written informed consent with the Ethics Statement before entry into their houses for mosquito aspiration. This study used solely mosquito samples, while no human participants were involved.

Area of study and mosquito collection

Cx. quinquefasciatus mosquitoes were collected across a transect from Eastern to Southwestern Uganda, with collection points located in four districts (Fig 1A, S1 and S2 Tables). Collection points were located in Tororo and Kanungu districts (located in the Eastern and Western administrative regions) with populations of 517,000 and 252,000 people respectively and where more than 80% of dwellings are rural, and in Jinja and Kampala (Central region) with populations of 76,000 and 1.5 million people respectively [30, 31]. Samples were from a peri-urban area in Jinja (63% of houses in rural areas) and in the urban centre of Kampala.

Fig 1

Geographic distribution of insecticide target-site mutations Ace1-G119S and Vgsc-L1014F in field-collected Ugandan Cx. quinquefasciatus mosquitoes.

(A) Geographic location of mosquito collections and frequency of target-site mutations. Pie charts depict the relative frequency of Ace1-G119S and Vgsc-L1014F mutations. Image sourced from https://commons.wikimedia.org/wiki/Uganda via the United States Central Intelligence Agency (https://www.cia.gov/library/publications/the-world-factbook/geos/ug.html). The image is originally published under the CC0 public domain license (B) Bar charts show the genotypic frequency of target-site loci.

From each collection site adult male and female mosquitoes were collected by aspiration from inside 10–15 houses (up to 10 mosquitoes per house) between July and August 2012. Then, a sub-set of the samples (S2 Table) were used for microsatellite and target-site mutation genotyping as described below.

Collected mosquitoes are representative of three Ugandan Demographic and Health Survey (UDHS) regions; Central, South-western and Mid-eastern with similar levels of insecticide-treated net (ITNs) coverage i.e. percentage of householders with at least one ITN ranging from 56 to 58.6% reported in 2011 [30] (S1 Table). However, the insecticide usage pattern in Uganda has changed rapidly in recent years, therefore the data provided (S1 Table) reflects very recent use and it does not absolutely reflect the last 5–10 years when resistance may have developed [32, 33]. Therefore, the ITN and IRS figures from 2011 for each UDHS are used herein as a proxy of likely selection strength among the collection sites, although they may not represent precisely the sources and strength of selection in each location. Indeed, recently Abeku et al. [34] have demonstrated that the effect of ITN usage on Kdr allele frequency is impacted by annual rainfall.

Across Uganda, two main pyrethroid based ITNs have been distributed by the government; Permanet 2.0 (polyester coated with deltamethrin) and Olyset Nets (polyethylene with permethrin incorporated). Additionally, a distinct level of ITN usage and Indoor Residual Spraying (IRS) is reported throughout the UDHS regions from where samples were collected (S1 Table). While major IRS programmes have not been reported in Jinja or Tororo, in Kanungu, approximately 45,000 households covering a population of 190,000 were sprayed with λ-cyhalothrin in 2007 as part of the malaria control programme [35]. In Kampala, a programme of environmental management was undertaken in the early 2000s (http://health.go.ug/programs/national-malaria-control-program) but no sustained insecticidal vector control has been used.

DNA isolation and species identification

Genomic DNA from individual mosquitoes was isolated using a DNeasy kit (Qiagen) following the manufacturer’s recommendations. All samples were confirmed as Cx. quinquefasciatus by a diagnostic PCR assay which types diagnostic length polymorphisms in the second intron of the acetylcholinesterase-2 gene to differentiate Cx. quinquefasciatus from Cx. pipiens [36].

Allelic discrimination assays of insecticide target-site mutations

To facilitate design of target-site mutation genotyping assays, partial genomic fragments spanning the location of the G119S mutation in the Ace-1 gene and the L1014F mutation in the Vgsc gene were amplified from genomic DNA isolated from mosquitoes collected from the four populations. A TaqMan assay targeting the Ace-1 gene was designed to genotype the SNP (GGC/AGC) in the first base of the codon at position 119, while a Vgsc-pyrosequencing assay to type the L1014F mutation in exon 20 [6] of the Vgsc gene was developed to detect two synonymous resistant alleles (TTC and TTT) and the wild-type allele (TTA) (S1 Methods). We chose the pyrosequencing method for genotyping Vgsc-L1014F instead of applying the widely used Taqman assay due to the existence of three alleles at this codon, which limits the use of Taqman for genotyping in a single reaction. Further details for PCR amplification of the partial fragment of both genes, as well as the target-site assay primer and probe sequences and reaction conditions are provided in the S1 Methods.

In the studied populations, genotyping of both Ace-1 and the Vgsc indicated the presence of gene duplication. Such gene duplications hinder calculation of many population-genetic statistics such as allele frequencies, linkage disequilibrium (LD) and population differentiation especially when there is polymorphism for copy number. Indeed, the absence of straightforward genotyping and data analysis methods limited our ability to address further the implications of the gene duplication in the evolution of resistance, with similar drawbacks also reported in other studies [15, 37]. Due to the limitations of the genotyping methods (e.g. it is not possible to differentiate the genotypes RRSS from RRRS), our analysis is based solely on the presence/absence of resistant and susceptible alleles as applied to a single copy gene. This assumption could result in an underestimation of insecticide associated-resistant alleles within the studied populations whilst still allowing us to infer a proxy of increased resistance allele frequency in a population under selection pressure.

Detection of Ace-1 gene duplication by haplotype diversity

The absence of Ace-1 homozygous resistant genotypes (AGC/AGC) observed through Taqman (see Results) indicated a likely gene duplication of Ace-1 in all four populations studied. To investigate the possible presence of multiple gene copies, a partial fragment of Ace-1 was PCR amplified from 12 individuals from Kampala. After cloning (S1 Methods), between six and eight colonies from each individual were sequenced in order to detect the presence of >2 alleles, indicative of gene duplication [15].

Sequences from each individual were aligned in CodonCode Aligner software version 4.2.2 with ClustalW [38] and visualized using Jalview [39]. MEGA 5.1 [40] was used to analyse haplotype variability by calculating the number of polymorphic sites and nucleotide diversity (π). Frequencies and relationships between haplotypes were visualised by a minimum spanning network tree generated using the program PopArt available at http://popart.otago.ac.nz.

Microsatellite genotyping

To date few microsatellite markers have been described for Cx quinquefasciatus [41-44], limiting the scope of population genetic studies. Herein, the microsatellite genotyping was conducted using a novel panel of 30 microsatellites loci developed for this study (S2 Methods), which have been combined into five multiplex reactions, thereby enhancing the throughput screening for genetic diversity in Cx. quinquefasciatus.

The newly designed microsatellite markers were isolated by scanning 180 Cx. quinquefasciatus supercontigs downloaded from VectorBase [45] with SciRoKo [46]. Individual mosquitoes were typed using five six-plex PCR reactions followed by genotyping using a Beckman-Coulter CEQ8000 capillary electrophoresis system with a 400 size standard kit. Genotypes were sized using the Beckman-Coulter CEQ 2000 DNA analysis system software and manually verified. Microsatellite genotype data were analysed with the program Micro-Checker [47] to detect possible scoring errors (stutter peaks and allele drop-out) and null alleles.

Microsatellite data analysis

Genetic diversity

Following genotyping of 162 mosquitoes at 30 microsatellite loci, four loci were rejected due to presence of null alleles (see Results). Thus, 26 loci were selected for inferences of genetic diversity and population structuring. Allelic frequencies, observed (HO) and expected (HE) heterozygosities were calculated using GenAIEX 6.5 [48]. Allelic richness, adjusted to the smallest population sample size, was estimated in FSTAT [49], and Polymorphic Information Content (PIC) calculated using Cervus 3.0 [50]. Genotypic frequencies were tested for deviation from Hardy-Weinberg equilibrium (HWE) for each locus by the exact probability test available in GENEPOP 4.3 [51], followed by a sequential Bonferroni correction.

Genetic structure

Genetic differentiation among populations was estimated by pairwise F using Arlequin 3.5.1 [52]. In addition, an AMOVA analysis was also carried out in Arlequin to estimate the level of differentiation among populations from different clusters based on F values. The pattern of migration under an isolation-by-distance model was tested with a Mantel’s test using distance linearized F (F/(1-F)) and geographic distance [53]. Geographic distance corresponding to the great-circle distance, which corresponds to the shortest distance over the earth’s surface, was calculated between pairs of geographic coordinates (using http://www.movable-type.co.uk/scripts/latlong.html).

A Bayesian analysis to infer the population structure without prior information of the geographic distribution of samples was carried out using STRUCTURE 2.3 [54]. To identify the optimal number of clusters (K) in these populations twenty independent runs were conducted for each K value (ranging from K = 1 to K = 7) with 10,000 interactions and 100,000 replications. The most likely K value was calculated for each run by the log likelihood (LnP(D)) method and results compiled using CLUMPP [55]. Population structure was also evaluated by Discriminant Analysis of Principal Components (DAPC) using the adegenet R package [56] [57]. To identify the optimal number of clusters for the DAPC clustering, k-means values were sequentially tested, and then k values compared using Bayesian Information Criterion (BIC), with the lowest value of BIC used as the likely number of clusters

Gene flow

To infer levels of gene flow among populations we utilized a Bayesian approach implemented in BAYESSASS 3.0 [58]. Initially, five independent runs were performed applying the parameters 107 MCMC, 106 burn-in, chain sampling every 2,000 generations, and all mixing parameters set as default (0.10), with the exception of both allele frequencies (DA) and inbreeding coefficient (DF) runs set at 0.30 as suggested by the author. Parallel runs differed only in initial seed value. Mixing and convergence of MCMCs were visually assessed using TRACER 1.6. Among the five independent runs, we chose the one with lowest Bayesian deviation and constant convergence to assess migration patterns.

Detection of loci under selection

The 26 microsatellite loci studied were tested for selection based on the F outlier approach. This method evaluates the relationship between F and H (expected heterozygosity) under the assumption of neutrality. The test was performed under an island model and infinite allele model as implemented in LOSITAN [59]. Analyses were performed applying the “Neutral mean F” to allow the program to run once to determine a first candidate subset of selected loci in order to remove them from the computation of the neutral F, and then the option “Force mean F” was also used to get a simulated mean F close to the observed one found in the real dataset. All other options or settings were left as default and 50,000 simulations were run.

To investigate if differing selection strength is detectable between geographic locations the F outlier method was run, with all populations considered separately and with Central Ugandan populations (largely urban) combined. Thus, four simulations were performed as described above for: (a) all populations, (b) populations defined here as Central Uganda (Jinja and Kampala) vs Tororo, (c) Central Uganda vs Kanungu and (d) Kanungu vs Tororo. Where signatures of genomic selection in the vicinity of likely candidate microsatellite selected markers were detected, we examined the number of genes present in the supercontig encompassing the microsatellite markers and their GO terms using the BioMart tools available at VectorBase (www.vectorbase.org).

Results

Frequency of target-site resistant alleles

Genotyping of two target-site mutations, Ace1-119S and Vgsc-1014F, was conducted utilizing TaqMan allelic discrimination and pyrosequencing assays, respectively (S1 Results). The Ace-1 resistant allele (AGC) was observed in all populations with a frequency ranging from 12–32% (Fig 1A). Despite the relatively high frequency of the 119S allele, we observed a complete absence of homozygous resistant mosquitoes (S1 Fig). The frequency of heterozygotes was almost three times higher in Kampala and Jinja compared to the other two populations (Fig 1B). Among the four populations a moderate He (0.323 ± 0.065) was observed with a marked excess of heterozygotes (FIS = -0.362), whilst within populations, genotype frequencies were not significantly different from Hardy-Weinberg expectation for both the Kanungu and Tororo populations (Table 1). Additionally, no significant difference in resistance allele frequencies was detected between Jinja and Kampala (P = 1.0) and between Kanungu and Tororo (P = 1.0).

Table 1

Ugandan Culex quinquefasciatus populations and genotypes of the target-site mutations Ace1-G119S and Vgsc-L1014F.

	Ace1-G119S								Vgsc-L1014F
			Genotype								Genotype
Population	N	f 119S	GG	GS	SS	He	FIS	HWa	N	f 1014F	LL	FC L	FC FC	FTL	FCFT	FT FT	He	FIS	HWa
Jinja	40	0.32	14	26	0	0.439	-0.481	0.002**	34	0.62	5	13	7	3	5	1	0.611	-0.011	0.945NS
Kampala	46	0.32	17	29	0	0.432	-0.46	0.002**	46	0.62	6	16	7	7	8	2	0.649	-0.05	0.954NS
Kanungu	45	0.12	34	11	0	0.215	-0.139	0.350NS	42	1	0	0	6	0	26	10	0.495	-0.249	0.106NS
Tororo	43	0.12	33	10	0	0.206	-0.132	0.388NS	43	0.97	1	1	19	0	20	2	0.45	-0.085	0.000***

N is the number of mosquitoes analysed; He is the expected heterozygosity and F is the inbreeding coefficient. HWa, P-value of χ tests for Hardy-Weinberg equilibrium.

f 119S and f 1014F, correspond to the allele frequency of the Ace1-119S and Vgsc-1014F resistant alleles.

GG, GS and SS correspond to homozygous and heterozygous genotypes for the Ace1-G119S locus.

FC; codon TTC, FT; codon TTT, L; codon TTA are distinct alleles at the Vgsc-L1014F locus.

NS, not significant;

**P < 0.01;

***P < 0.001.

The Vgsc-1014F locus in Cx. quinquefasciatus has two alternative non-synonymous substitutions (TTT and TTC) at the third position of the codon (L1014F) [60]. For the Vgsc-1014F mutations we detected a high frequency of resistant alleles ranging from 62–100% (Fig 1A). The pyrosequencing genotyping identified 11 individuals harbouring all three kdr alleles instead of the two expected for diploid organisms (S3A Fig). Mosquitoes with tri-allelic genotypes were observed in all populations with the exception of Kanungu, with a frequency of 4.0–10.8% (S3B Fig). Excluding individuals with this tri-allelic pattern to allow the use of standard methods to infer departures from Hardy-Weinberg, our analysis detected significant deviation from HWE only in mosquitoes from Tororo (Table 1).

Significant differences in Vgsc resistance allele frequencies were detected between all populations with the exception of the Jinja and Kampala comparison (P = 0.63). The TTC resistant allele was observed at a high frequency (41–69%) in all populations with the exception of Kanungu, where a higher frequency of the TTT resistant allele (55%) was detected, coupled with an absence of the wild-type allele (Fig 1A). In all populations, we observed a low frequency of susceptible homozygotes (ranging from 0 to 14%), while the majority of resistant alleles were observed in heterozygotes (Fig 1B).

Ace-1 haplotype diversity and identification of allelic copy variation

A 535 bp fragment of the Ace-1 gene was sequenced from 12 individuals (N = 57 sequences with between six and eight colonies sequenced from each individual; GenBank accession numbers: KT591708-KT591765) with nine haplotypes detected displaying a total of 10 polymorphic sites (S3 Table) and haplotype frequencies ranging from 0.017–0.397. The resistant allele was detected in four haplotypes: A, E, H and I accounting for 40% of total haplotypes identified. The minimum spanning network tree (S3C Fig) shows that three resistance bearing haplotypes (A, H, I) differ from each other by a single mutational step with Haplotype A the most common. The remaining 119S haplotype (E) is equally distant from both haplotypes A and B differing by three mutational steps, while it differs by only two mutational steps from haplotype D. For wildtype Ace1-119G, three haplotypes B, C and G were detected with similar frequencies whereas haplotypes D and F were the least frequent. From 12 individuals studied through cloning and sequencing we identified three individuals with more than two distinct Ace-1 haplotypes, indicative of duplication. These three mosquitoes displayed from three to five haplotypes in the neighbor-joining dendrogram (S3D Fig).

A total of 186 alleles were detected across the 30 microsatellite loci (S4 Table, S2 Results) with 19 private alleles detected in 17 of the loci. The number of unique alleles in the sample from Kampala was twice that detected in any other population (S2 Fig); however, the presence of these eight private alleles are unlikely to significantly influence the genetic structure analysis since only two alleles have a frequency higher than 1%. The allelic richness and Polymorphic Information Content (PIC) across markers varied from 2.88 (MCQ2) to 11.14 (MCQ21) and from 0.332 (MCQ2) to 0.825 (MCQ26), respectively (S4 Table).

Hardy-Weinberg exact tests indicated that 26 out of 120 instances did not conform to HW equilibrium after multiple-testing correction. However, only two loci (MCQ28 and MCQ33) showed departure from HW equilibrium across all populations (S5 Table). The number of loci that were not in HW equilibrium in each population ranged from 5 to 8 (S2 Table). For all loci with departure from HW equilibrium, significant deviation was associated with a positive inbreeding coefficient (F), revealing deviation in the direction of heterozygote deficiency. Analyses performed with Micro-Checker indicated a high frequency of null alleles at loci MCQ28 (0.12) and MCQ33 (0.22), which could likely be responsible for the heterozygote deficiency at these loci. Although loci MCQ1 and MCQ34 were within HW expectation in at least one population, these loci failed to amplify in more than 10% of the samples also suggesting the presence of null alleles. We therefore excluded four markers (MCQ1, MCQ28, MCQ33 and MCQ34) of the 30 tested. The microsatellite exclusion level in this study (four of 30 microsatellites or 13%) is consistent with previous studies which have examined genotypes in family lines as an additional verification of utility (e.g. [43, 61]).

After removal of these four markers from the analysis due to the likely presence of null alleles, levels of genetic diversity, allelic diversity, expected heterozygosity (HE) and allelic richness within populations were similar across populations (S2 Fig, S5 Table). The mean number of alleles across all loci varied from 4.96 in samples from Kanungu to 5.42 in samples from Kampala, with the highest allelic richness (Rs = 5.17) detected in Kampala. Nevertheless, no significant difference in the index of genetic diversity was observed among the populations (S2 Fig). Average locus H across all samples ranged from 0.597 in Kampala to 0.612 in Kanungu with a median of 0.604 (S2 Table).

Population structure

Pairwise F estimates were used to measure the genetic differentiation among the four Ugandan populations. The overall F value was relatively low but statistically significant (F = 0.019, P = 0.001), while the pairwise F values between populations ranged from 0.0071 to 0.0293 (S6 Table), with all figures significant at P < 0.05. The DAPC and Structure analysis (S4 Fig) also indicated population structure, which grouped the four populations within three clusters K = 3, with individuals from Jinja and Kampala belonging to the same cluster and this being distinct from the other two samples (Fig 2A and 2B). AMOVA analysis carried out using the four populations showed that the majority of variation exists within populations (83%), while the variation among individuals and among populations are 15 and 2% respectively (S7 Table).

Fig 2

Genetic differentiation estimates among the four Cx. quinquefasciatus populations based on allele frequencies at 26 microsatellite markers.

(A) First and second PCs of the Discriminant Analysis of Principal Components (DAPC) to infer population structure. The graph represents the individuals as dots and the groups as coloured ellipses, which model 95% of the corresponding variability plotted. Populations (ellipses) are plotted within the orthogonal space defined by the first two PCA eigenvalues (inserts). (B) Bayesian analysis of population structure performed using STRUCTURE. Diagrammatic representation of population clusters for the most likely K (K = 3) with each coloured segment corresponding to the proportion of individuals assigned to a hypothetical population or cluster. Each vertical bar represents an individual with the height of the column segments showing the probability of assignment of belonging to one of the three clusters. Dark lines correspond to population assignment based on geographic collection-site.

Levels of migration between the studied populations were measured under Bayesian inference and this substantiated the existence of contemporary migration in the studied populations. The proportion of individuals originating within each population varied from 67.5% to 93.5%, with the highest admixture in Kampala and Kanungu which have the lowest proportions of individuals originating within the respective population (Table 2). Non-symmetrical rates of gene flow between populations, ranging from 0.71% - 14.5% of individuals originating from a different population, are also shown in Table 2. Posterior probability of migration among populations was inferred using the run with the lowest estimate of Bayesian deviance (among five independent runs, S5A–S5E Fig), with mixing parameter adjusted as migration rates; 0.1, allele frequency; 0.3 and inbreeding coefficients; 0.30.

Table 2

Matrix of inferred gene flow between Ugandan Cx. quinquefasciatus populations.

Gene Flow	Jinja	Kampala	Kanungu	Tororo
Jinja	0.9349(0.0334)	0.0080(0.0079)	0.0085(0.0083)	0.0486(0.0322)
Kampala	0.2911(0.0201)	0.6745(0.0076)	0.0079(0.0079)	0.0265(0.0175)
Kanungu	0.1725(0.0342)	0.0079(0.0076)	0.6750(0.0082)	0.1446(0.0339)
Tororo	0.1150(0.0508)	0.0079(0.0077)	0.0092(0.0090)	0.8678(0.0510)

Values in the form mij representing the proportion of individuals in the ith population that originated from the jth population per generation. Values in parentheses are standard deviations of the posterior probability distributions.

Analysis of isolation-by-distance also supports a pattern of migration, with migration between populations positively correlated with the geographic distance (R2 = 0.963, P = 0.047) based on microsatellite markers (Fig 3). In contrast to microsatellite markers, for both insecticide resistance loci Vgsc-1014F and Ace1-119S, no significant association was detected between the frequency of resistant alleles and a sample’s geographical location (Fig 3).

Fig 3

Relationship between geographic distance versus genetic distance (F/1- F), analysed in parallel for the microsatellite, Vgsc-L1014F and Ace1-G119S loci.

Significance was tested using a Mantel test. Grey shading corresponds to the 95% confidence interval. Geographic distance corresponding to the great-circle distance, which corresponds to the shortest distance over the earth’s surface between two points.

Loci under selection

The 26 loci were also tested for selection based on the F outlier approach, from which only marker MCQ21 showed strong signs of divergent selection when comparing all populations, or in comparisons of the two distinct sub-sets of populations based on population genetic structuring—Central Uganda (Jinja and Kampala) compared to Tororo and Kanungu. No signal of selection was observed for MCQ21 when comparing only Tororo and Kanungu (Fig 4). Marker MCQ11 showed a signal of divergent selection only when comparing Central Uganda and Tororo (Fig 4). Supercontigs containing these markers (3.186 for MCQ21 and 3.526 for MCQ11) contain 47 and 22 genes respectively in the vicinity of these markers with a genomic window resolution of 500 kb and 347 kb, respectively (S3 and S4 Results). No genes with recognised links to insecticide resistance were observed (e.g. detoxification gene families such as the cytochrome P450 (P450s), esterases and glutathione S-transferases (GSTs))

Fig 4

F outlier approach for the 26 studied microsatellite loci versus a model of neutral evolution in Ugandan Cx. quinquefasciatus.

Central populations correspond to a cluster of Jinja and Kampala (based on genetic structure analysis, see Fig 3). Each dot represents a microsatellite marker. The markers located in the grey middle area are assumed to be neutral. The markers located below the neutral area are candidates for being under balancing selection and the markers located above the neutral area are candidates for being under directional selection. X-axis: estimated heterozygosity values. Y-axis: F-values. Analyses were performed assuming an infinite allele model.

Discussion

Resistance alleles at the target sites of insecticides have known deleterious effects [13, 62] and are therefore not expected to persist in the absence of selection. Our data which indicate copy number variation at both the Vgsc-L1014F and Ace1-G119S loci and allele frequency differences between populations for both loci suggest differential insecticide selection pressure on Cx. quinquefasciatus populations across Uganda. The significantly higher frequency of Vgsc-1014F in comparison to Ace1-119S for all populations has been reported previously in Culex and other vector mosquitoes on the African continent [63-65], and might result from past and ongoing reliance on pyrethroid-based vector control approaches to mitigate the burden of vector-borne diseases. However, it is important to bear in mind that our allele frequency estimates are likely under-reporting true frequencies, particularly for the Vgsc, as the duplication observed in Ugandan populations imposes limitations on available methods to infer allele frequency as reported elsewhere [37].

Whilst recognising these limitations, our data show a contrasting pattern of Ace1-119S and Vgsc-1014F allele frequencies across populations, likely indicating heterogeneous selection pressure. For both markers, we were not able to identify a clear spatial geographic spread, but we did detect heterogeneity in allele frequencies of these selected markers when comparing Uganda Central to Eastern and Southwest population. Indeed, the highest frequency of Vgsc-1014F was detected outside the Central region, while the highest frequency of Ace1-119S was observed with a pattern of geographic distribution opposite to that of Vgsc-1014F. Therefore, these figures highlight the challenges for implementation of efficient control intervention in highly fragmented habitats including overlapping of urban and rural patches or due to heterogeneous urbanization structure, which could result in a mosaic of landscapes with sharp shifts of insecticide selection pressure and frequency of insecticide-resistance associated alleles [66-68].

Markedly, although our study demonstrates high frequencies of both target-site mutations (Vgsc-1014F and Ace1-119S) in Ugandan Cx. quinquefasciatus populations no vector control programme targeting this species has been conducted. Therefore, these results support a hypothesis that the evolution of insecticide resistance in these populations is driven by control interventions targeting other vector species or pest insects. Indeed, previous studies have suggested that the evolution of insecticide resistance in Cx. quinquefasciatus on the African continent might result from intervention approaches targeting sympatric Anopheles species [27, 69]. For instance, a study of Cx. quinquefasciatus from Zambia showed a significant increase of the Vgsc-1014F mutation and increase of detoxification enzyme activity assays after introduction of ITNs for malaria control [27].

In Uganda in particular, we suggest that variation of Vgsc-1014F frequency between districts might reflect the combination of approaches applied to reduce Anopheles populations, including national distribution of ITN and IRS only in highly endemic and epidemic-prone regions. For example, across Uganda irregular spraying with insecticides such as lambda-cyhalothrin, DDT, alpha cypermethrin and latterly bendiocarb [70] has been reported since 2006 when IRS was re-introduced for vector control [71]. Regional differences in ITN coverage may also contribute to the Vgsc-1014 allele distribution e.g. the Uganda Bureau of Statistics [72] report that householders that own at least one ITN differ by 44% between the East Central and West Nile regions in 2011.

Markedly, we have also noticed a distinct pattern of ITN usage among UDHS, which differ by 14% between Central-South western and 8.7% Central-Mid-eastern, despite a similar ITN coverage of around 58% for all three regions (S1 Table). Additionally, the contrasting pattern of Vgsc-1014 markers among the studied districts might also reflect the drastic shift of ITN distribution, which increased by 26% (Central region), 41.6% (Mid-eastern) and 32.8% (South-western) between 2006 and 2009 [32, 33].

Nevertheless, with our study design it is not possible make a direct association of the UNMCP (Uganda National Malaria Control Program) interventions with the evolution of resistance in Culex populations. Indeed, we recognise that pinpointing the particular source underlying the evolution of resistance remains a difficult task. For example, increased resistance in Cx. quinquefasciatus populations from Ghana and Benin, with high frequencies of Ace1-119S was attributed to the use of organophosphates and/or carbamate insecticides by farmers for pest control or pollutants in mosquito breeding sites [66, 67]. Heavy use of domestic insecticide was also associated with increased insecticide resistance in Cx. quinquefasciatus populations from Ghana [73].

Regardless of the strength of selection imposed by diverse insecticide-based approaches such as ITNs, IRS, pesticides and environment pollutants [23], the evolutionary pattern of increased resistance in Cx. quinquefasciatus from Uganda could also be modulated by variation in the repertoire of adaptive genetic mechanisms, as in the studied population we identified gene duplication for both target-site loci and allelic variation of the Vgsc resistant allele. Interestingly, we identify two alternative resistant alleles (TTT and TTC) in the Vsgc-1014 codon observed at high frequency in all populations, contrasting with previous reports on the geographic distribution of Vgsc-1014F in Culex mosquitoes worldwide [11], which exclusively detect the TTT allele [74-76]. To date, the TTC variant has been detected at high frequency only in Cx. quinquefasciatus from Sri Lanka, which was previously argued as the most likely geographic origin [60], and was only recently detected at very low frequency in African Cx. quinquefasciatus from Zanzibar [69]. Although we identified co-occurrence of both resistant variants of Vgsc it was not possible to determine if heterozygote TTT/TTC increases a mosquito’s level of resistance, or offsets possible deleterious effects.

We also identify a significant difference in frequency of the Vgsc-1014F resistant alleles among populations, with the TTC variant being predominant in all populations with the exception of Kanungu. We do not have a conclusive explanation for the difference in frequency of TTT and TTC among Ugandan populations and between Uganda and other African populations. Nevertheless, a few potential mechanisms could explain this. One possible explanation involves a colonization process with the TTC allele introduced to Uganda from another geographical region with a large founder effect followed by genetic drift that shifted the TTT and TTC allele frequency, although a de novo origin and dispersion cannot be discounted. Alternatively, it is also possible that the TTC allele arose on the African continent but has not been detected in other populations due to the application of specific Culex species Vgsc-1014 genotyping assays, which were mostly designed for detecting the TTT and TTA alleles [77].

Our target-site genotyping of the Ace1-119 and Vgsc-1014 variants also indicated the presence of gene duplication for both genes, suggesting that CNV could be among the mechanisms driving the evolution of insecticide in these populations. Indeed, gene duplication of insecticide-associated genes such as Ace-1, esterase estα21 and estβ21, Rdl (resistance to dieldrin) and detoxification genes (e.g. GSTs, P450s) have been reported as adaptive resistance mechanisms in diverse arthropod species including Tetranychus urticae, Aphis gossypii and A. gambiae [3, 14, 16, 37, 78].

The likely duplication of the Vgsc was indicated by the presence of 6.25% of mosquitoes genotyped by pyrosequencing displaying three alleles simultaneously for Vgsc-1014. For Ace-1, duplication was evidenced by a haplotype analysis, which identified three out of 12 individuals exhibiting either 3 or 5 distinct haplotypes simultaneously. Whilst contamination or DNA from spermatheca could be mistaken for duplications we have shown that the duplication of Vgsc is present in males [79] and we see no triallelic genotypes in the microsatellite data. This duplication could explain the large variation observed within resistant and susceptible alleles as demonstrated by the minimum spanning Network (S3C Fig), which contrasted with an expected low haplotype diversity due to a selective sweep effect on a locus under selection.

The absence of Ace-119S homozygotes in all four populations, was not surprising as the same pattern has been reported previously in both Culex and Anopheles populations [16, 20, 80]. In both species, this pattern has been linked to the deleterious effects of the homozygous resistant genotype due to changes in acetylcholinesterase kinetic properties and mutation fitness cost in an insecticide-free or reduced insecticide exposure environment. Alternatively, the absence of resistant homozygotes might also arise from duplication of the Ace-1 gene involving resistant and susceptible alleles, creating a permanent heterozygosis to partially normalise AChE enzyme activity levels [37, 81]. In contrast to the proposed compensatory adaptive mechanism of the Ace-1 duplication to overcome deleterious effects and fitness disadvantage of the resistant Ace1-119S [19, 20], the likely adaptive role played by the recently reported Vgsc duplication in Cx. quinquefasciatus from the USA (Xu et al. 2011) and in A. aegypti from Brazil [22] is currently unknown.

Our study also detected a contrasting pattern between microsatellite and insecticide selected markers. For both markers linked to resistance, it was not possible to identify a clear spatial geographic pattern, but we did detect heterogeneity in allelic frequencies when comparing Ugandan Central (largely urban) to Eastern and Southwest populations (largely rural). By contrast, no significant difference in genetic diversity indices was observed at neutral markers among the populations, although these data did indicate a structured pattern. Although the low levels of genetic variation of microsatellites across Ugandan populations from distinct geographic locations have also been reported in Cx. pipiens from California, USA [82], heterogeneity in diversity levels across Uganda was expected since vector control interventions have been applied annually, and might be expected to result in population bottlenecks. Such reductions of genetic diversity between populations under different strengths of selection pressure have been reported in Cx. quinquefasciatus from Recife, Brazil after vector control intervention with Bacillus sphaericus, a bio-larvicide [83]. As indicated by these authors, discrepancies between neutral and selected markers might be linked to differences in the marker’s genetic properties that modulate changes in allele frequencies, such as drastic bottleneck effects and differences in intensity of local selection with frequent gene flow [84, 85].

The microsatellite data also indicated significant genetic structure among populations with the three clusters identified by the STRUCTURE and DAPC analyses corresponding to Eastern, Central and Southwest Uganda regions. The level of population structure among our study sites was consistent with findings from similar studies in Culex mosquitoes using molecular markers as diverse as ISSR (Inter-Simple Sequence Repeat) and microsatellites, which reported that population structure reflects geoclimatic, geographic distance and environmental variables [86-88]. Additionally, based on microsatellite markers an intense migration among Ugandan populations with significant patterns of isolation by distance (IBD) was detected, which could prevent loss of genetic diversity and consequently could explain the homogenous microsatellite genetic diversity observed across the studied populations.

The distinct pattern of IBD observed for microsatellite markers and selected markers, also indicated that factors other than geographic distance, such as environmental heterogeneity as well as serial sequential founder effects, may also influence the pattern of spatial autocorrelation of resistance-associated markers [89]. For example, in Aedes rusticus, comparisons among areas treated and not-treated with Bacillus thuringiensis israelensis (Bti) show more intense IBD among treated than among non-treated sites [90]. These authors argue that the difference in observed IBD instead of spatial distribution could result from differences in population size after insecticide application or lower migration capability of selected mosquitoes imposed by a fitness cost of the resistance mechanism.

Therefore, despite the intense gene flow identified in the studied populations, it is possible that local differences in strength of insecticide selection might have created heterogeneous “insecticide adaptive islands” across the geographic regions. Analysis of microsatellite locus neutrality indicated that all markers indeed adhered to a neutral model with the exception of the MCQ21 across all populations and MCQ11 when comparing only Ugandan Central to Eastern populations. For both markers, further analysis of the genomic location of markers does not indicate likely linkage with known insecticide associated genes such as Ace-1, Vgsc or genes encoding detoxification enzymes (e.g. esterase, GSTs and P450s). However, this analysis was limited by the incomplete Cx. quinquefasciatus genome assembly, which restricted the genomic window capable of being applied in the analysis. Due to the contrasting patterns of neutrality for marker MCQ21, we speculate that this marker could be hitchhiking with an insecticide-associated gene nearby this marker, with a genomic location beyond the supercontig resolution.

In conclusion, our study shows that the contemporary pattern of local selection is driven by the evolution of insecticide resistance in Ugandan Cx. quinquefasciatus occurring despite the absence of vector control intervention targeting this species. This study demonstrates that although it is challenging to pinpoint the direct source of insecticide exposure in a background of selection pressure, studies of the evolution of insecticide in non-target species can provide important information about the nature of the insecticidal challenges to which vectors are exposed.

Allelic discrimination assays of insecticide target-site mutations in Culex quinquefasciatus.

(PDF)

Click here for additional data file.

Development of microsatellite multiplex panels for population genetic analysis in Culex quinquefasciatus mosquitoes.

(PDF)

Click here for additional data file.

Cx. quinquefasciatus Ugandan Ace1-119 and Vgsc-1014 genotyping.

(XLSX)

Click here for additional data file.

Cx. quinquefasciatus Ugandan microsatellite genotyping data.

(XLSX)

Click here for additional data file.

Genomic information for supercontig 3.186.

(XLS)

Click here for additional data file.

Genomic information for supercontig 3.526.

(XLS)

Click here for additional data file.

Ace-1 target-site mutation (G119S) genotyping in Cx. quinquefasciatus.

(PDF)

Click here for additional data file.

Genetic diversity estimates across Ugandan Cx. quinquefasciatus populations based on 26 microsatellite markers.

(PDF)

Click here for additional data file.

Detection of copy number variation in the Vgsc and Ace-1 gene in Cx. quinquefasciatus.

(PDF)

Click here for additional data file.

Clustering analysis in field-collected Ugandan Cx. quinquefasciatus using Bayesian assignment implemented in the STRUCTURE software.

(PDF)

Click here for additional data file.

BAYESASS 3.0 trace-plot and analysis parameters from five parallel runs to detect the lowest Bayesian deviance and examine convergence.

(PDF)

Click here for additional data file.

Location of field-caught mosquitoes, demographic index of collection sites and main features of vector control intervention in Uganda 2011.

(PDF)

Click here for additional data file.

Sample sizes and indices of genetic diversity at 26 microsatellite loci for four Ugandan Cx. quinquefasciatus populations.

(PDF)

Click here for additional data file.

Ace-1 haplotype diversity based on sequences obtained from Cx. quinquefasciatus mosquitoes sampled in Uganda.

Polymorphic sites identified in a 535 bp partial fragment (intron 2 and exon 3) of Ace-1.

(PDF)

Click here for additional data file.

Microsatellite locus characterization across all Cx. quinquefasciatus populations.

(PDF)

Click here for additional data file.

Summary of microsatellite variation in different Ugandan populations of Cx. quinquefasciatus.

(PDF)

Click here for additional data file.

Pairwise F estimates (Weir & Cockerham 1984) across Ugandan Culex quinquefasciatus populations.

(PDF)

Click here for additional data file.

Analysis of molecular variance (AMOVA) for four Culex quinquefasciatus populations from Uganda, using 26 microsatellite markers.

(PDF)

Click here for additional data file.