Landscape genetic structure and evolutionary genetics of insecticide resistance gene mutations in Anopheles sinensis
© Chang et al. 2016
Received: 4 November 2015
Accepted: 14 April 2016
Published: 23 April 2016
Anopheles sinensis is one of the most abundant vectors of malaria and other diseases in Asia. Vector control through the use of insecticides is the front line control method of vector-borne diseases. Pyrethroids are the most commonly used insecticides due to their low toxicity to vertebrates and low repellency. However, the extensive use of insecticides has imposed strong selection pressure on mosquito populations for resistance. High levels of resistance to pyrethroid insecticides and various mutations and haplotypes in the para sodium channel gene that confers knockdown resistance (kdr) have been detected in An. sinensis. Despite the importance of kdr mutations in pyrethroid resistance, the evolutionary origin of the kdr mutations is unknown. This study aims to examine the evolutionary genetics of kdr mutations in relation to spatial population genetic structure of An. sinensis.
Adults or larvae of Anopheles sinensis were collected from various geographic locations in China. DNA was extracted from individual mosquitoes. PCR amplification and DNA sequencing of the para-type sodium channel gene were conducted to analyse kdr allele frequency distribution, kdr codon upstream and downstream intron polymorphism, population genetic diversity and kdr codon evolution. The mitochondrial cytochrome c oxidase COI and COII genes were amplified and sequenced to examine population variations, genetic differentiation, spatial population structure, population expansion and gene flow patterns.
Three non-synonymous mutations (L1014F, L1014C, and L1014S) were detected at the kdr codon L1014 of para-type sodium channel gene. A patchy distribution of kdr mutation allele frequencies from southern to central China was found. Near fixation of kdr mutation was detected in populations from central China, but no kdr mutations were found in populations from southwestern China. More than eight independent mutation events were detected in the three kdr alleles, and at least one of them evolved multiple times subsequent to their first divergence. Based on sequence analysis of the mitochondrial COI and COII genes, significant and large genetic differentiation was detected between populations from southwestern China and central China. The patchy distribution of kdr mutation frequencies is likely a consequence of geographic isolation in the mosquito populations and the long-term insecticide selection.
Our results indicate multiple origins of the kdr insecticide-resistant alleles in An. sinensis from southern and central China. Local selection related to intense and prolonged use of insecticide for agricultural purposes, as well as frequent migrations among populations are likely the explanations for the patchy distribution of kdr mutations in China. On the contrary, the lack of kdr mutations in Yunnan and Sichuan is likely a consequence of genetic isolation and absence of strong selection pressure. The present study compares the genetic patterns revealed by a functional gene with a neutral marker and demonstrates the combined impact of demographic and selection factors on population structure.
The development of insecticide resistance is regarded as an empirical model of an adaptive trait  whose origin and spread depend upon the selection pressure applied to the target population, the fitness cost associated with the resistant allele and the gene flow among populations . Geographical landscape is often considered as an important factor limiting gene flow among populations [3–6]. However, the impact of gene flow barriers imposed by landscape differences and local selection on the spatial distribution of resistance alleles is largely unknown.
Resistance to insecticides in mosquitoes involves multiple mechanisms, including the target-site insensitivity caused by mutations in the para sodium channel gene (knockdown resistance kdr) and the detoxification of insecticides in mosquitoes involving detoxification enzymes (i.e. P450 monooxygenases, glutathione-S-transferases and esterases) that metabolise the insecticide before it reaches its target (metabolic resistance) [7–19]. Mutations in the para sodium channel gene have been shown to be associated with the use of pyrethroids, the most commonly used class of insecticide for malaria control worldwide . The spectrum of mutations at the para sodium channel gene is highly conserved across insect species, indicating convergent evolution . There are two hypotheses that relate to the distribution of kdr mutations: either the mutations arose once and then spread, or have evolved independently in different populations across a broad landscape. For instance, in Anopheles gambiae, the main African malaria vector, non-synonymous mutations at position L1014 of the para sodium channel gene occurred at least 4–5 times independently [1, 20]. Some of the mutation events at the para sodium channel gene were derived from single mutational steps from the common ancestor whereas others were the result of two mutational steps [1, 20]. Multiple origins of kdr mutations were also found in Musca domestica and in the aphid Myzus persicae [21, 22]. An alternative situation was seen in the agricultural pest Bemisia tabaci, where the L925I and T929V mutations appeared to have originated once and then spread through migration and global trade [2, 23].
In An. sinensis, the most abundant and important malaria vector species in southeast Asia [24–32], four non-synonymous mutations at codon L1014 of the para sodium channel gene are known, including the L1014F [31, 33–37], L1014S [17, 35], L1014C [31, 33, 34, 36, 37] and L1014W  mutations. There could be various factors influencing the distribution of kdr allele variants in An. sinensis across China and Southeast Asia. For example, the complexity of landscape could impose barriers to gene flow between the mosquito populations. Selection based upon the intensity and duration of pyrethroid use could also determine allelic distribution of kdr. On the other hand, genetic variation based on neutral markers reflects population structure resulting from demographic factors. In this case, the cytochrome c oxidase genes of the mitochondrial genome are used as the neutral reference. Based on the unique characteristics of maternal inheritance, no recombination, high variability relative to nuclear DNA, and low effective population size, the mitochondrial DNA (mtDNA) has been a marker of choice for studies of genetic diversity and population structure.
The central goal of this study was to understand how local selection and migration have shaped the distribution of kdr mutations. The specific objectives were to (i) determine the kdr allele distribution in multiple sites over a large geographic region; (ii) elucidate the evolutionary origin(s) of kdr mutations; and (iii) examine the role of landscape and demographic history of An. sinensis on the evolution of kdr mutations using mitochondrial DNA sequence data.
Sample collection and preparation
Summary of Anopheles sinensis specimen collection sites in China
adult & larva
adult & larva
Amplification and sequencing of the mitochondrial cytochrome c oxidase genes
The COI and COII genes of the mitochondrial genome were amplified, respectively, using degenerate primers 1809 F (5'-CMC TTT CAT CTG GAA TTG CT-3') and 2708R (5'-AAA AAT GTT GAG GGA ARA ATG TTA-3') for COI (~900 bp) and primers COIIF (5'-TCT AAT ATG GCA GAT TAG TGC A-3') and COIIR (5'-ACT TGC TTT CAG TCA TCT AAT G-3') for COII (~800 bp). Reaction mix consisted of 10.5 μl 2X SYBR® Green PCR Master Mix (Thermo Fisher Scientific Inc., Waltham, MA, USA), 1 μl of template DNA (5–20 ng), 0.5 μl of 10 μM each primer, and 8.5 μl water. Conditions for polymerase chain reaction (PCR) were 95 °C for 3 min, 35 cycles of 94 °C for 30 s, 55 °C for 30 s, and 72 °C for 60 s, followed by an additional extension step at 72 °C for 6 min. Products of PCR were visualised on 1 % agarose gel and then purified and sequenced from both ends with the ABI Big Dye Terminator Cycle Sequencing Kit.
Amplification and sequencing of a 1285 bp fragment of the para sodium channel gene
We designed a pair of primers (Forward primer 63 F: 5'-GAC GTT CGT GCT CTG CAT TA-3', Reverse primer 1347R: 5'-GAG CGA TGA TGA TCC GAA AT-3') based on the An. sinensis sodium channel gene sequence (GenBank accession: DQ334052) to amplify a 1285 bp fragment that encompasses the kdr target codon position L1014, an upstream intron (intron-1) and a downstream intron (intron-2). Both primers are located in the coding region. Two internal primers, 566 F (5'-CAC TGC TGT AAA ACC CTG TGT-3') and 855R (5'-CTG TTT GCT AGG CAG TTT GC-3') were used for sequencing. We amplified the 1285 bp fragment from individual specimen and a total of 273 specimens were examined. Reaction mix consisted of 10.5 μl 2X SYBR® Green PCR Master Mix (Thermo Fisher Scientific Inc., Waltham, MA, USA), 1 μl of template DNA (5–20 ng), 0.5 μl of 10 μM each primer, and 8.5 μl water. PCR cycles involved an initial denaturing step at 95 °C for 3 min, 35 cycles of 94 °C for 1 min, 55 °C for 1 min, and 72 °C for 2 min, followed by an additional extension step at 72 °C for 6 min. Amplifications were performed in a Bio-Rad MyCycler Thermal Cycler. Products of PCR were visualised on 1 % agarose gel and then purified and sequenced from both ends with ABI Big Dye Terminator Cycle Sequencing Kit.
This study addressed the central question as to whether the patchy kdr allele distribution in An. sinensis populations from southern China to central China was shaped by geographic isolation and/or by insecticide selection. We asked two questions in particular, i.e., (i) What is the kdr allele frequency distribution pattern and the phylogenetic relationship of kdr haplotypes? and (ii) What is the spatial genetic structure of An. sinensis mosquito populations in China? We used the para sodium channel gene sequence data to address the first question used mtDNA sequence data to address the second.
Phylogenetic relationships of kdr haplotypes
Sequenced 1285 bp fragments of the para sodium channel gene encompassing the kdr codon L1014, intron-1 and intron-2 were aligned in Bioedit version 7.0 . DNA polymorphism including the number of segregating sites, number of haplotypes, haplotype diversity and nucleotide diversity were calculated using DnaSP v5 . The two intron regions were tested for neutrality by Tajima‘s D , Fu and Li’s D* and F*  and F tests . Partition of the genetic variation among and within populations was examined by AMOVA test implemented in Arlequin v3.5 . A Bayesian approach based on a priori predictions from the coalescent theory implemented in Phase 2.1.  was used to reconstruct haplotypes from the population genotypic data in mosquitoes heterozygous at more than one site within the amplified 1285 bp fragment. These haplotypes were confirmed by cloning of the 1285 bp fragment in a subset of 20 individuals. Amplified PCR products from the selected 20 individuals were cloned into pCR2.1TOPO TA vectors (Invitrogen, Carlsbad, CA, USA). Three to five clones were sequenced for each individual to identify haplotypes.
A statistical parsimony network was constructed using the TCS v1.21 software  to infer the genealogical relationships among haplotypes based on an analysis of 1285 bp fragments in the samples. The software calculates the frequencies of the haplotypes in the sample and estimates haplotype outgroup probabilities, which correlate with haplotype age . The pairwise comparisons of haplotypes was conducted using the parsimony algorithm described in Templeton et al.  to estimate the number of mutational steps between haplotypes using a connection limit at two steps and a connection probability threshold of 0.95.
Spatial genetic structure of An. sinensis
COI and COII gene sequences were aligned using Bioedit version 7.0 . The number of variable sites and haplotypes, as well as haplotype and nucleotide diversity were estimated using Arlequin v3.5  and DnaSP v5 . To determine the genetic structure of the mosquito populations, we used analysis of molecular variance (AMOVA) to partition genetic variation among groups (F CT), among populations within groups (F SC) and among populations among groups (F ST). Pairwise F ST and Ф ST among all populations were calculated. Significance values were adjusted for multiple comparisons using the false discovery rate (FDR) approach . Pairwise genetics distances were used to construct unrooted UPGMA tree in Phylip 3.695  to infer population clustering. Isolation-by-distance among populations was explored with Mantel tests  using the Isolation by Distance Web Service v3.23 . Euclidian geographical distances among the 15 populations were calculated based on geographical coordinates using great circle calculations to accurately represent distances on a spherical earth. Pairwise genetic distance between populations was calculated as F ST/(1- F ST) with 10,000 randomizations .
To examine the proportion of total genetic variance between groups of populations and identify possible genetic barriers between them, the spatial genetic structure of haplotypes (combined COI and COII sequences) was analysed using the Spatial Analysis of Molecular Variance (SAMOVA 2.0 software) . This software implements an approach, which combines genetic differentiation and geographical distance to define groups of populations that are maximally differentiated from each other, without the constraint for the geographic composition of the groups. The number of initial conditions was set to 100 simulated annealing processes for each configuration of K groups, with K ranging from 2 to 10, and run for 10,000 iterations. For each K, the configuration with the largest F CT values after the 100 independent simulated annealing processes was retained as the best grouping of populations.
Genetic landscape shape (GLS) interpolation analysis
We performed the genetic landscape shape (GLS) interpolation analysis based on the pairwise genetic and geographical distance matrices using the Alleles-In-Space (AIS) software . AIS is a computer program for the joint analysis of inter-individual spatial and genetic information. Detailed patterns of spatial genetic structure across the complete sampled area were visualised using the ‘genetic landscape shape’ interpolation procedure. This procedure produced a 3-dimensional surface plot where X- and Y-axes correspond to geographical locations and surface heights (Z-axes) represent genetic distances. Genetic structure across the landscape was inferred from measured genetic distances using an inverse distance weighted interpolation across a uniform grid laid over the entire sampling area. GLS was made by using a 80 × 80 grid and a raw genetic distance, with a distance weighting parameter α = 1.
Demographic history of the 15 populations was elucidated based on the neutrality indices Tajima’s D, Fu’s Fs and Ramos-Onsins and Rozas’s R2 statistic in DnaSP. Significantly negative D, Fs or small R2 values indicated recent population expansion, whereas positive or high values indicated population decline. Population expansion parameter τ (τ = 2 t* μ; where t = time in years, μ = mutation rate per locus) was estimated using the moment method, assuming the infinite sites model . Deviations from a model of population expansion were evaluated by computing the statistical significance of sums of squared deviation (SSD) and Harpending's raggedness index (r) over 1000 simulated samples of pairwise nucleotide differences. All P-values were corrected for multiple testing using the false discovery rate (FDR) approach .
Migration rate between populations
Gene flow and migration rates between all pairwise populations were estimated using the Bayesian coalescence-based approach implemented in LAMARC version 2.1.10 . Felsenstein‘84 (F84) substitution model was used with empirical base frequencies [61, 62]. Default prior settings were used for Theta and Bayesian estimation of migration rates. LAMARC analysis consisted of 3 simultaneous searches with automatically adjusted heating temperature using 10 initial chains, with 500 samples and a sampling interval of 20 steps, followed by 2 final chains, with 10,000 samples and sampling interval of 20 steps. One thousand samples were discarded for initial and final burn-in. Migration rate was measured as 4 Nm, a multiplication of LAMARC‘s M and Theta values for the recipient population.
The PCR amplification of para-type sodium channel gene generated a 1285 bp PCR fragment, including a 187 bp exon encompassing the kdr target codon position L1014, its upstream intron 905 bp and exon 74 bp, and its downstream intron 64 bp and exon 55 bp. The PCR amplification of the COI and COII genes of the mitochondrial genome resulted in 900 bp and 774 bp fragments, respectively. In order to elucidate the kdr allele frequency distribution pattern and the phylogenetic relationships of kdr haplotypes, the para sodium channel gene sequence data were analysed for kdr allele frequency distribution, kdr codon upstream and downstream intron polymorphism, population genetic diversity and kdr codon evolution. Because mtDNA is haploid and maternally inherited, isolated populations will drift to different haplotype frequencies faster and achieve approximately twice the level of differentiation comparing to nuclear markers. So, the mtDNA sequence data were analysed for population variations, genetic differentiation, spatial population structure, population expansion and gene flow patterns.
Kdr allele frequency distribution
Three non-synonymous mutations (L1014F, L1014C, and L1014S) were detected at the kdr codon L1014. Overall, the distribution of kdr mutant allele frequencies varied by geographic location. Populations in central China (Anhui, Hubei, Hunan and Jiangsu) reached a near fixation for kdr at codon L1014 mutations, whereas the southern populations (Guangdong, Fujian, Guangxi, Hainan and Guizhou) exhibited a low frequency of kdr mutations (Fig. 1). We did not detect any kdr mutation in the southwestern populations (Yunnan and Sichuan provinces). In central China, L1014F was the predominant mutation and L1014C was also common in three provinces (Henan, Hubei and Anhui) in central China (Fig. 1).
Kdr codon upstream and downstream intron polymorphism
Polymorphism of kdr intron and neutrality test of 15 Anopheles sinensis populations
Fu and Li’s D*
Fu and Li’s F*
The overall haplotype diversity (Hd) and nucleotide diversity (Pi) were 0.82 and 0.07, respectively. Genetic diversity varied significantly among geographical regions (Table 2). In central China, four of the five populations showed low haplotype (0.30 ± 0.05; F (2,11) = 46.3, P < 0.0001) and nucleotide diversity (0.014 ± 0.003, F (2,11) = 47.3, P < 0.0001), whereas in the south and southwest populations, high haplotype (0.92 ± 0.005 and 0.85 ± 0.13) and nucleotide diversity (0.07 ± 0.005 and 0.06 ± 0.03) were detected.
The Tajima’s D and Fu and Li's D* tests indicated no significant departures from neutrality in all the populations. Similarly, only one test of Fu and Li‘s F* indicated significant departure from neutrality in all the populations. However, significant departures from neutrality were detected by Fu’s Fs in eight populations and seven of them were from the south and southwest (Table 2), which indicated a pattern of haplotype diversity that deviated from the expectation under neutrality. Since the Fu‘s Fs statistic is particularly sensitive to demographic effects, it is difficult to conclude whether positive selection or demographic history (e.g., population expansion) account for the observed pattern. Although the AMOVA indicated that most of the variation (~70 %) occurred among populations among regions (F ST = 0.28, P < 0.001), the percentage of variation found among populations within the same region (3.2 %) as well as among the three regions (27.9 %) was also highly significant (F SC = 0.04, P < 0.001; F CT = 0.18, P < 0.001). Both the F ST and Ф ST values (based only on haplotype frequencies) were significant for 61 and 55 out of the 105 pairwise comparisons, respectively, after Bonferroni correction for multiple comparisons (Additional file 1: Table S1). The majority of the differences were detected between the central China populations and the remaining populations.
kdr codon evolution
Mitochondrial DNA sequence variation
A total of 341 sequences for COI and 303 sequences for COII were generated for the 15 populations (Additional file 1: Table S3). COI and COII sequences were aligned giving a total length of 814 and 704 bp, respectively. The COI sequence alignment showed a total of 146 variable sites, of these 78 were parsimony informative. The COII sequence alignment revealed 85 variable sites, including 41 parsimony informative. Nucleotide diversity, number of haplotypes and haplotype diversity were higher in COI (p = 0.88 × 10−2; h = 238; Hd = 0.99; Additional file 1: Table S3) than in COII (p = 0.47 × 10−2; h = 142; Hd = 0.93). Among the 15 populations, An. sinensis from southern China exhibited higher haplotype and nucleotide diversity in COI and COII than those in central China (Additional file 1: Table S3).
Genetic differentiation among populations
AMOVA found a significant population structure in the 15 An. sinensis populations examined, based on the combined COI and COII sequences (F ST = 0.027, P < 0.0001). Population differentiation was significant between the southwest and south/central China regions (F ST = 0.069, P < 0.001). Both the F ST and Ф ST values (based only on haplotype frequencies) were significant for 30 and 26 out of the 105 pairwise comparisons, respectively, after Bonferroni correction for multiple comparisons (Additional file 1: Table S4). The majority of the differences were detected between the two Yunnan populations (YNLH and YNNE) and the remaining populations.
Spatial genetic structure analysis
Demographic history and neutrality test based on mtDNA sequences
Analysis of historical demography. Variation of estimated demographic parameters based on mitochondrial COI and COII sequence data of 15 populations of An. sinensis
Migration and gene flow patterns
Near fixation of kdr mutation was detected in populations from central China, but no kdr mutations were found in populations from Yunnan and Sichuan provinces. Absence of kdr mutations in An. sinensis samples from southwestern China and low frequency of kdr mutations in southern China is likely a consequence of geographical isolation in the mosquito populations combined with the absence of strong local selection.
Geographical isolation of An. sinensis populations from southwestern China
Based on mtDNA sequencing, significant genetic differentiation was detected between Yunnan (southwestern China) and populations from south/central China. The AMOVA analysis indicated that the percentage of genetic variation between these two regions (~28 %) is much higher than that among populations within each of the regions (~3 %). Such regional genetic differentiation was also manifested by the observation of more private alleles in each of the regions than within individual populations of a region. Our data show that An. sinensis in Yunnan has a relatively stable population size over time without drastic demographic changes, unlike populations from central and southern China. The strong genetic differentiation and limited gene flow between the Yunnan and central China populations suggest that these two regions are genetically isolated. SAMOVA analysis further confirmed that there is a genetic barrier restricting gene flow between the Yunnan populations and the other populations examined. These patterns were also consistent with the phylogenetic trees of populations. The two Yunnan populations are located in the Yunnan Plateau, a mountainous area in southwest China with high elevations in the northwest and low elevations in the southeast of Yunnan province. Therefore, both physical distance and heterogeneous landscape could be factors inhibiting gene flow between Yunnan and south/central China. However, based on the analysis of DNA sequences surrounding kdr codon L1014, significant genetic differentiation was detected between the central China populations and the remaining populations, suggesting a strong local selection of the kdr resistance allele in the central China populations.
Effects of insecticide-associated selection
For the two Yunnan populations, neutrality tests (Fu’s Fs statistic) did not provide evidence of selection at the kdr introns. This finding concurs with a relatively low rate of insecticide use in agriculture and indoor spraying in Yunnan as compared to central China. Apart from limited migrations, the absence of strong selective pressure may in part explain the lack of insecticide-resistance alleles in Yunnan. Lower genetic variation was observed in the introns of the kdr locus in populations from central China when compared to the other populations. These populations showed a high frequency of kdr mutation (L1014F), suggesting that low variation at the adjacent introns could be a consequence of a recent selective sweep associated with intense insecticide use following the selection and/or fixation of the resistance genotypes. A selective sweep occurs when an allele rapidly increases its frequency due to positive or directional selection. Through genetic hitchhiking, the frequency of linked alleles in the flanking regions of the locus under selection can also increase, thus reducing genetic variation [65, 66].
It must be noted that among the different neutrality tests, Fu’ F S statistic shows the greatest power in detecting departures from neutrality under a genetic hitchhiking model . The high Fu’s F S estimates for central China could reflect the impact of strong selection by the increased use of insecticide not only on the kdr locus, but also on the other gene regions. In the past half-century, insecticides are known to have been extensively used for agricultural and vector control purposes in central China . Rice is the major agricultural crop in these study sites, with 1–2 harvests per year. Due to severe insect pest damage, insecticide use for pest control has been very intensive, with several rounds of sprays administered during each growing season. Since the mid-1980s, pyrethroids have been the dominant insecticides with pyrethroids-treated areas constituting more than one third of the total insecticide-treated areas in central China [37, 67, 68]. In addition to their agricultural use, pyrethroids have had various public health applications, e.g., as indoor sprays or incense, impregnated in bed nets, or as tools in public sanitation [25, 37, 67, 68]. The large scale agricultural use and indoor residual spraying of pyrethroid insecticides have probably generated selective pressure at the sodium channel gene in local An. sinensis populations. Anopheles sinensis from central China were found with extremely high resistant to pyrethroids, organochlorine, organophosphates and carbamates [31, 33]. Thus, it is not surprising for selection to have operated beyond the upstream and downstream introns of the kdr locus and towards parts of the mosquito’s genome in central China.
Multiple origins of kdr mutations
Information on whether resistance-associated mutations evolved independently vs single emergency and then spatially spread is valuable to elucidate the relative importance of mutation and migration in the spread of insecticide resistance. We detected three non-synonymous mutations at codon L1014 of the sodium channel gene. The majority of samples from central China harbored the kdr mutation L1014F, whereas the L1014S mutation was only found in Guangdong Province, southern China. The L1014C mutation was detected in both southern and central China. This result is consistent with previous study  of kdr genotyping in An. sinensis which found that the L1014F mutation was largely distributed in central China. Analysis of an upstream intron (intron-1) and a downstream intron (intron-2) of the kdr codon suggests at least eight independent origins of kdr alleles in An. sinensis populations from south to central China. Multiple origins of kdr alleles were also reported in the Afrotropical mosquito vector An. gambiae . The kdr haplotypes H01-1014 F and H01-L1014C from the same ancestor were the most widespread in An. sinensis populations from central China and evolved multiple times subsequent to their first divergence. This phenomenon could be explained by the selection pressure of insecticides experienced in these populations, which favors kdr resistant alleles. Therefore, the frequency of haplotype (H01-L1014) carrying the wild-type allele is much lower than that of haplotype (H01-L1014F) harbouring kdr mutations. However, the kdr haplotype L1014S from two distinct ancestors was detected only in the Guangdong population, suggesting a recent origin of this mutation that is unique to the populations in southern China.
Our results indicate multiple origins of the kdr insecticide-resistant alleles in An. sinensis from southern and central China. Local selection related to intense and prolonged use of insecticide for agricultural purposes, as well as frequent migrations among populations are likely the explanations for the patchy distributions of kdr mutations in China. On the contrary, the lack of kdr mutations in Yunnan and Sichuan is likely a consequence of genetic isolation and absence of strong selection pressure. The present study compares the genetic patterns revealed by a functional gene with a neutral marker and demonstrates the combined impact of demographic and selection factors on population structure.
We thank Ying Wang, Jianhua Duan, Zhengbo He, Jiangyan Li, Yiji Li, Lin Zhou, Tielong Xu and Fengyang Fu for assistance with mosquito collection, two anonymous reviewers for critical review and valuable suggestions, Amruta Dixit and Elizabeth Hemming for their editorial assistance. This work is supported by grants from the U.S. National Institutes of Health (U19 AI089672, R03TW008237 and D43TW009527) and from the Key Projects of Science Research in Education Institutes of Anhui Colleges and Universities (KJ2015A146), the Science Research Innovation Team Projects of Anhui Colleges and Universities (2016–40).
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
- Pinto J, Lynd A, Vicente JL, Santolamazza F, Randle NP, Gentile G, Moreno M, Simard F, Charlwood JD, do Rosário VE, et al. Multiple origins of knockdown resistance mutations in the afrotropical mosquito vector Anopheles gambiae. PLoS ONE. 2007;2(11), e1243.View ArticlePubMedPubMed CentralGoogle Scholar
- Gauthier N, Clouet C, Perrakis A, Kapantaidaki D, Peterschmitt M, Tsagkarakou A. Genetic structure of Bemisia tabaci Med populations from home-range countries, inferred by nuclear and cytoplasmic markers: impact on the distribution of the insecticide resistance genes. Pest Manag Sci. 2014;70(10):1477–91.View ArticlePubMedGoogle Scholar
- Coster SS, Babbitt KJ, Cooper A, Kovach AI. Limited influence of local and landscape factors on finescale gene flow in two pond-breeding amphibians. Mol Ecol. 2015;24(4):742–58.View ArticlePubMedGoogle Scholar
- Cushman S, Shirk A, Landguth E. Landscape genetics and limiting factors. Conserv Genet. 2013;14(2):263–74.View ArticleGoogle Scholar
- Igawa T, Oumi S, Katsuren S, Sumida M. Population structure and landscape genetics of two endangered frog species of genus Odorrana: different scenarios on two islands. Heredity. 2013;110(1):46–56.View ArticlePubMedPubMed CentralGoogle Scholar
- Storfer A, Murphy MA, Evans JS, Goldberg CS, Robinson S, Spear SF, Dezzani R, Delmelle E, Vierling L, Waits LP. Putting the ‘landscape’ in landscape genetics. Heredity. 2006;98(3):128–42.View ArticlePubMedGoogle Scholar
- WHO. Global plan for insecticide resistance management in malaria vectors. Geneva: World Health Organization; 2012. Available at http://www.who.int/malaria/publications/atoz/gpirm/en/.Google Scholar
- Pasteur N, Raymond M. Insecticide resistance genes in mosquitoes: their mutations, migration, and selection in field populations. J Hered. 1996;87(6):444–9.View ArticlePubMedGoogle Scholar
- Wondji CS, Priyanka De Silva WA, Hemingway J, Ranson H, Parakrama Karunaratne SH. Characterization of knockdown resistance in DDT- and pyrethroid-resistant Culex quinquefasciatus populations from Sri Lanka. Trop Med Int Health. 2008;13(4):548–55.View ArticlePubMedGoogle Scholar
- Brengues C, Hawkes NJ, Chandre F, McCarroll L, Duchon S, Guillet P, Manguin S, Morgan JC, Hemingway J. Pyrethroid and DDT cross-resistance in Aedes aegypti is correlated with novel mutations in the voltage-gated sodium channel gene. Med Vet Entomol. 2003;17(1):87–94.View ArticlePubMedGoogle Scholar
- Dong K, Du Y, Rinkevich F, Nomura Y, Xu P, Wang L, Silver K, Zhorov BS. Molecular biology of insect sodium channels and pyrethroid resistance. Insect Biochem Mol Biol. 2014;50:1–17.View ArticlePubMedPubMed CentralGoogle Scholar
- Donnelly M, Corbel V, Weetman D, Wilding C, Williamson M, Black W. Does kdr genotype predict insecticide-resistance phenotype in mosquitoes. Trends Parasitol. 2009;25(5):213–9.View ArticlePubMedGoogle Scholar
- Martinez-Torres D, Chandre F, Williamson MS, Darriet F, Berge JB, Devonshire AL, Guillet P, Pasteur N, Pauron D. Molecular characterization of pyrethroid knockdown resistance (kdr) in the major malaria vector Anopheles gambiae s.s. Insect Mol Bio. 1998;7(2):179–84.View ArticleGoogle Scholar
- Nwane P, Etang J, Chouaibou M, Toto J, Koffi A, Mimpfoundi R, Simard Fl. Multiple insecticide resistance mechanisms in Anopheles gambiae s.l. populations from Cameroon, Central Africa. Parasit Vectors. 2013;6(1):41.View ArticlePubMedPubMed CentralGoogle Scholar
- Nwane P, Etang J, Chouaibou M, Toto J, Mimpfoundi R, Simard F. Kdr-based insecticide resistance in Anopheles gambiae s.s populations in Cameroon: spread of the L1014F and L1014S mutations. BMC Res Notes. 2011;4:463.View ArticlePubMedPubMed CentralGoogle Scholar
- Ranson H, Jensen B, Vulule JM, Wang X, Hemingway J, Collins FH. Identification of a point mutation in the voltage-gated sodium channel gene of Kenyan Anopheles gambiae associated with resistance to DDT and pyrethroids. Insect Mol Bio. 2000;9(5):491–7.View ArticleGoogle Scholar
- Verhaeghen K, Van Bortel W, Trung HD, Sochantha T, Keokenchanh K, Coosemans M. Knockdown resistance in Anopheles vagus, An. sinensis, An. paraliae and An. peditaeniatus populations of the Mekong region. Parasit Vectors. 2010;3(1):59.View ArticlePubMedPubMed CentralGoogle Scholar
- Kasai S, Scott J. Over expression of cytochrome P450 CYP6D1 is associated with monooxygenases-mediated pyrethroid resistance in house flies from Georgia. Pestic Biochem Physiol. 2000;68:34–41.View ArticleGoogle Scholar
- Liu N. Insecticide resistance in mosquitoes: impact, mechanisms, and research directions. Annu Rev Entomol. 2015;60(1):537–59.View ArticlePubMedGoogle Scholar
- Etang J, Vicente JL, Nwane P, Chouaibou M, Morlais I, Do Rosario VE, Simard F, Awono-Ambene P, Toto JC, Pinto J. Polymorphism of intron-1 in the voltage-gated sodium channel gene of Anopheles gambiae s.s. populations from Cameroon with emphasis on insecticide knockdown resistance mutations. Mol Ecol. 2009;18(14):3076–86.View ArticlePubMedGoogle Scholar
- Anstead JA, Williamson MS, Denholm I. Evidence for multiple origins of identical insecticide resistance mutations in the aphid Myzus persicae. Insect Biochem Mol Biol. 2005;35(3):249–56.View ArticlePubMedGoogle Scholar
- Rinkevich FD, Hedtke SM, Leichter CA, Harris SA, Su C, Brady SG, Taskin V, Qiu X, Scott JG. Multiple origins of kdr-type resistance in the house fly, Musca domestica. PLoS ONE. 2012;7(12), e52761.View ArticlePubMedPubMed CentralGoogle Scholar
- Alon M, Benting J, Lueke B, Ponge T, Alon F, Morin S. Multiple origins of pyrethroid resistance in sympatric biotypes of Bemisia tabaci (Hemiptera: Aleyrodidae). Insect Biochem Mol Biol. 2006;36(1):71–9.View ArticlePubMedGoogle Scholar
- Chareonviriyaphap T, Bangs MJ, Ratanatham S. Status of malaria in Thailand. Southeast Asian J Trop Med Public Health. 2000;31(2):225–37.PubMedGoogle Scholar
- Cui F, Raymond M, Qiao CL. Insecticide resistance in vector mosquitoes in China. Pest Manag Sci. 2006;62(11):1013–22.View ArticlePubMedGoogle Scholar
- Jung J, Jung Y, Min GS, Kim W. Analysis of the population genetic structure of the malaria vector Anopheles sinensis in South Korea based on mitochondrial sequences. Am J Trop Med Hyg. 2007;77(2):310–5.PubMedGoogle Scholar
- Jung J, Lee E, Kim W. Isolation and characterization of polymorphic microsatellite markers of Anopheles sinensis, a malaria vector mosquito in the East Asia region. Mol Ecol Notes. 2006;6(4):1272–4.View ArticleGoogle Scholar
- Liu XB, Liu QY, Guo YH, Jiang JY, Ren DS, Zhou GC, Zheng CJ, Zhang Y, Liu JL, Li ZF, et al. The abundance and host-seeking behavior of culicine species (Diptera: Culicidae) and Anopheles sinensis in Yongcheng city, People‘s Republic of China. Parasit Vectors. 2011;4(1):221.View ArticlePubMedPubMed CentralGoogle Scholar
- Ma Y, Yang M, Fan Y, Wu J, Ma Y, Xu J. Population structure of the malaria vector Anopheles sinensis (Diptera: Culicidae) in China: two gene pools inferred by microsatellites. PLoS ONE. 2011;6(7), e22219.View ArticlePubMedPubMed CentralGoogle Scholar
- Sinka M, Bangs M, Manguin S, Chareonviriyaphap T, Patil A, Temperley W, Gething P, Elyazar I, Kabaria C, Harbach R, et al. The dominant Anopheles vectors of human malaria in the Asia-Pacific region: occurrence data, distribution maps and bionomic precis. Parasit Vectors. 2011;4(1):89.View ArticlePubMedPubMed CentralGoogle Scholar
- Zhong D, Chang X, Zhou G, He Z, Fu F, Yan Z, Zhu G, Xu T, Bonizzoni M, Wang MH, et al. Relationship between knockdown resistance, metabolic detoxification and organismal resistance to pyrethroids in Anopheles sinensis. PLoS ONE. 2013;8(2), e55475.View ArticlePubMedPubMed CentralGoogle Scholar
- Zhu G, Xia H, Zhou H, Li J, Lu F, Liu Y, Cao J, Gao Q, Sattabongkot J. Susceptibility of Anopheles sinensis to Plasmodium vivax in malarial outbreak areas of central China. Parasit Vectors. 2013;6(1):176.View ArticlePubMedPubMed CentralGoogle Scholar
- Chang X, Zhong D, Fang Q, Hartsel J, Zhou G, Shi L, Fang F, Zhu C, Yan G. Multiple resistances and complex mechanisms of Anopheles sinensis mosquito: a major obstacle to mosquito-borne diseases control and elimination in China. PLoS Negl Trop Dis. 2014;8(5), e2889.View ArticlePubMedPubMed CentralGoogle Scholar
- Kang S, Jung J, Lee S, Hwang H, Kim W. The polymorphism and the geographical distribution of the knockdown resistance (kdr) of Anopheles sinensis in the Republic of Korea. Malaria J. 2012;11(1):151.View ArticleGoogle Scholar
- Tan WL, Li CX, Wang ZM, Liu MD, Dong YD, Feng XY, Wu ZM, Guo XX, Xing D, Zhang YM, et al. First detection of multiple knockdown resistance (kdr)-like mutations in voltage-gated sodium channel using three new genotyping methods in Anopheles sinensis from Guangxi Province, China. J Med Entomol. 2012;49(5):1012–20.View ArticlePubMedGoogle Scholar
- Tan WL, Wang ZM, Li CX, Chu HL, Xu Y, Dong YD, Wang ZC, Chen DY, Liu H, Liu DP, et al. First report on co-occurrence knockdown resistance mutations and susceptibility to beta-cypermethrin in Anopheles sinensis from Jiangsu province, China. PLoS ONE. 2012;7(1), e29242.View ArticlePubMedPubMed CentralGoogle Scholar
- Wang Y, Yu W, Shi H, Yang Z, Xu J, Ma Y. Historical survey of the kdr mutations in the populations of Anopheles sinensis in China in 1996–2014. Malaria J. 2015;14(1):120.View ArticleGoogle Scholar
- Dong X. The mosquito fauna of Yunnan (Volume one), vol. 1. Kunming: Yunnan Publishing Group Corporation, Yunnan Science & Technology Press; 2010.Google Scholar
- Severson DW. RFLP analysis of insect genomes. In: Crampton JM, Beard CB, Louis C, editors. The molecular biology of insect disease vectors. London: Chapman & Hall; 1997. p. 309–20.View ArticleGoogle Scholar
- Gao Q, Beebe NW, Cooper RD. Molecular identification of the malaria vectors Anopheles anthropophagus and Anopheles sinensis (Diptera: Culicidae) in central China using polymerase chain reaction and appraisal of their position within the Hyrcanus group. J Med Entomol. 2004;41(1):5–11.View ArticlePubMedGoogle Scholar
- Joshi D, Park MH, Saeung A, Choochote W, Min GS. Multiplex assay to identify Korean vectors of malaria. Mol Ecol Resour. 2010;10(4):748–50.View ArticlePubMedGoogle Scholar
- Hall TA. A user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT. Nucleic Acids Symp Ser. 1999;41:95–8.Google Scholar
- Rozas J, Sanchez-DelBarrio JC, Messeguer X, Rozas R. DnaSP, DNA polymorphism analyses by the coalescent and other methods. Bioinformatics. 2003;19(18):2496–7.View ArticlePubMedGoogle Scholar
- Tajima F. Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics. 1989;123(3):585–95.PubMedPubMed CentralGoogle Scholar
- Fu YX, Li WH. Statistical tests of neutrality of mutations. Genetics. 1993;133(3):693–709.PubMedPubMed CentralGoogle Scholar
- Fu YX. Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics. 1997;147(2):915–25.PubMedPubMed CentralGoogle Scholar
- Excoffier L, Lischer HE. Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Mol Ecol Resour. 2010;10(3):564–7.View ArticlePubMedGoogle Scholar
- Stephens M, Smith NJ, Donnelly P. A new statistical method for haplotype reconstruction from population data. Am J Hum Genet. 2001;68(4):978–89.View ArticlePubMedPubMed CentralGoogle Scholar
- Clement M, Posada D, Crandall KA. TCS: a computer program to estimate gene genealogies. Mol Ecol. 2000;9(10):1657–9.View ArticlePubMedGoogle Scholar
- Donnelly P, Tavaré S. The ages of alleles and a coalescent. Adv Appl Probab. 1986;18:1–19.View ArticleGoogle Scholar
- Templeton AR, Crandall KA, Sing CF. A cladistic analysis of phenotypic associations with haplotypes inferred from restriction endonuclease mapping and DNA sequence data. III. Cladogram estimation. Genetics. 1992;132(2):619–33.PubMedPubMed CentralGoogle Scholar
- Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Series B Stat Methodol. 1995;57(1):289–300.Google Scholar
- Felsenstein J. PHYLIP (Phylogeny Inference Package) version 3.7a, Distributed by the author. Seattle: Department of Genome Sciences, University of Washington; 2009.Google Scholar
- Mantel N. The detection of disease clustering and a generalized regression approach. Cancer Res. 1967;27(2 Part 1):209–20.PubMedGoogle Scholar
- Jensen JL, Bohonak AJ, Kelley ST. Isolation by distance, web service. BMC Genet. 2005;6:13.View ArticlePubMedPubMed CentralGoogle Scholar
- Rousset F. Genetic differentiation and estimation of gene flow from F-statistics under isolation by distance. Genetics. 1997;145(4):1219–28.PubMedPubMed CentralGoogle Scholar
- Dupanloup I, Schneider S, Excoffier L. A simulated annealing approach to define the genetic structure of populations. Mol Ecol. 2002;11(12):2571–81.View ArticlePubMedGoogle Scholar
- Miller MP. Alleles in space (AIS): computer software for the joint analysis of interindividual spatial and genetic information. J Hered. 2005;96(6):722–4.View ArticlePubMedGoogle Scholar
- Rogers AR. Genetic evidence for a pleistocene population explosion. Evolution. 1995;49(4):608–15.View ArticleGoogle Scholar
- Kuhner MK. LAMARC 2.0: maximum likelihood and Bayesian estimation of population parameters. Bioinformatics. 2006;22(6):768–70.View ArticlePubMedGoogle Scholar
- Kishino H, Hasegawa M. Evaluation of the maximum likelihood estimate of the evolutionary tree topologies from DNA sequence data, and the branching order in hominoidea. J Mol Evol. 1989;29(2):170–9.View ArticlePubMedGoogle Scholar
- Felsenstein J. PHYLIP (Phylogeny Inference Package). Seattle: Department of Genetics, University of Washington; 1993.Google Scholar
- Stephens M, Donnelly P. A comparison of Bayesian methods for haplotype reconstruction from population genotype data. Am J Hum Genet. 2003;73(5):1162–9.View ArticlePubMedPubMed CentralGoogle Scholar
- Librado P, Rozas J. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics. 2009;25(11):1451–2.View ArticlePubMedGoogle Scholar
- Gentile G, Santolamazza F, Fanello C, Petrarca V, Caccone A, della Torre A. Variation in an intron sequence of the voltage-gated sodium channel gene correlates with genetic differentiation between Anopheles gambiae s.s. molecular forms. Insect Mol Biol. 2004;13(4):371–7.View ArticlePubMedGoogle Scholar
- Weill M, Chandre F, Brengues C, Manguin S, Akogbeto M, Pasteur N, Guillet P, Raymond M. The kdr mutation occurs in the Mopti form of Anopheles gambiae s.s. through introgression. Insect Mol Biol. 2000;9(5):451–5.View ArticlePubMedGoogle Scholar
- Huang J, Qiao F, Zhang L, Rozelle S. Farm pesticides, rice production, and human health in China. In: EEPSEA, editor. Research report no2001-RR3. Singapore: Economy and Environment Program for Southeast Asia; 2001. p. 1–58.Google Scholar
- Wang DQ, Xia ZG, Zhou SS, Zhou XN, Wang RB, Zhang QF. A potential threat to malaria elimination: extensive deltamethrin and DDT resistance to Anopheles sinensis from the malaria-endemic areas in China. Malaria J. 2013;12(1):164.Google Scholar