Diversity and frequency of kdr mutations within Anopheles sinensis populations from Guangxi, China

Background Anopheles sinensis is a major vector of malaria in China and its control is under great threat as the development of insecticide resistance. Voltage-gated sodium channel (VGSC) is the target of several classes of insecticides. Genetic mutations of VGSC have been documented to confer knockdown resistance (kdr) to dichlorodiphenyltrichloroethane (DDT) and pyrethroids in mosquitoes. To control this vector efficiently, it is important to know the resistance-associated genetic mutations, their distribution frequencies and genealogical relations. Methods Three hundreds and thirteen (313) adults of An. sinensis collected from nine locations across Guangxi Zhuang Autonomous Region were used. The partial sequence of the An. sinensis voltage gated sodium channel gene (AS-VGSC) containing codon 1014 was sequenced. PHASE2.1 was used to construct the haplotypes of each individual, and the accuracy of haplotypes was further confirmed by clone sequencing. The genealogical relations of kdr mutations in AS-VGSC was analysed using TCS 2.1 and Network 5.0. Results Sixteen AS-VGSC haplotypes including seven haplotypes carrying non-synonymous mutations at codon 1014, and fifty-five AS-VGSC genotypes were identified from 313 mosquitoes collected from nine geographical locations across Guangxi. The number of haplotypes in each of the nine populations ranged from 5 to 13. The frequency of haplotypes carrying kdr mutations ranged from 2.7 to 80.0 % within the nine populations, of which 1014C was unexpectedly high in the northeast of Guangxi. Genealogical analysis suggested multiple origins of kdr mutations in An. sinensis. Conclusion Diverse haplotypes of AS-VGSC are distributed in Guangxi. The presence of haplotypes carrying mutations at codon 1014 indicates a risk of pyrethroid and DDT resistance. The kdr mutations show differential distribution geographically, with high frequencies occurred in the northeast of Guangxi. Genealogical analysis suggests multiple origins of kdr mutations in An. sinensis populations in Guangxi. These findings have important practical implications for the sustainability of An. sinensis control programmes. Electronic supplementary material The online version of this article (doi:10.1186/s12936-016-1467-3) contains supplementary material, which is available to authorized users.

Guangxi Zhuang Autonomous Region was once a malaria-endemic region. Before 1949, there were more than 5 million malaria patients per year [20]. Thanks to the "National Malaria Control Programme", the "Basically Eliminating Malaria" strategy and the "Action Plan of Malaria Elimination (2010-2020)" implemented since 1955, 2000, and 2010 respectively in China [21,22], the malaria burden has been substantially reduced, with one indigenous and 2068 imported cases being reported from 2010 to 2015 in Guangxi [22]. Although Guangxi has already achieved remarkable accomplishments in eliminating malaria, the risk of malaria re-emergence remains partly due to increasing cross-border population migration and the unique natural environment (e.g. rice fields) suitable for mosquito breeding in Guangxi [22][23][24][25].
As the historical use in vector control and continuing use in agriculture of different insecticides, not only the geographical distribution and density of malaria vectors are likely to change, insecticide resistance is expected to be selected as well. Recent investigations have shown that An. sinensis has replaced Anopheles minimus to be the main malaria vector in Guangxi [24]. However, the status of insecticide resistance in An. sinensis is less understood in this region up to date. As an effort in this direction, the distribution and frequency of VGSC mutations that possibly lead to resistance to DDT and pyrethroids were investigated in nine field populations of An. sinensis collected extensively across Guangxi in this study.

Samples
Anopheles sinensis adults used in the study were caught by light trap (wave length 365 nm) from July to September in 2014, around farmers' houses at different geographical locations across Guangxi. The sampling sites were part of the malaria vector monitoring sites that were set up partly because there were imported malaria cases observed since 2010 in these regions. A brief description of the sampling sites was provided in Additional file 1. The farmers living at the sampling sites usually use mosquito-coil (largely containing pyrethroids) to prevent mosquitoes, and apply insecticides (e.g. pyrethroids) for crop protection.
The collected individual mosquitoes were first morphologically identified using the determination key provided by [26]. Anopheles sinensis adults were put into 100-μl Eppendorf tubes containing 100 % ethanol solution, and kept at 4 °C until use. The accuracy of morphological classification was guaranteed by molecular detection of ten randomly selected adults from each population using the rDNA-ITS2 method [27].

Haplotype identification and genealogical analysis
Genotype data were carefully inspected by using Chromas (Technelysium Pty Ltd, Australia). AS-VGSC haplotypes were constructed by software PHASE2.1 (http:// www.bioinf.man.ac.uk/resources/phase/) from the genotype data. The presence (and accuracy) of haplotypes constructed in heterozygotes was confirmed by sequencing of three to five clones derived from 12 representative heterozygotes. Genealogical relations among haplotypes were analysed using Network 5.0 [30] and TCS2.1 [31].

Identification of polymorphic sites in the AS-VGSC gene
The gene of An. sinensis voltage sensitive sodium channel (KE525266.1)contains 31 exons and 30 introns. A ~325 bp-in-length sequence was sequenced through direct sequencing of each purified PCR product. The partial sequence (267 bp), excluding the primer sequence and covering partial exon 19 (169 bp), intron 19 (64 bp) and partial exon 20 (34 bp), was used for further analysis (Fig. 1). Two synonymous mutations (C/G at PS1, T/C at PS9) were observed. Five polymorphic sites were found in the intron 19 (PS4-8). Non-synonymous mutations were identified in the second (T/G/C) and third nucleotide (G/T) of codon 1014 (PS2 and 3).

Haplotype diversity of AS-VGSC gene
From the sequence data of 313 individuals, sixteen AS-VGSC haplotypes were identified ( Table 1). The number of haplotypes in each of the nine populations ranged from 5 (Guigang) to 13 (Yulin). Three haplotypes (1014L2, 1014L3 and 1014L8) were detected in all the nine populations, while 1014L7 was uniquely distributed in Nanning. 1014L1 was the predominant haplotype for Nanning and Yulin, 1014L2 for Baise and Hechi, and 1014L8 was most abundant in Wuzhou (Table 2).
Point mutations at position 1014 of the domain II of the VGSC protein have been documented to confer kdr [5,6]. In this study, three types of kdr mutations were detected: TTT (1014F), TCG (1014S), and TGT (1014C). Accordingly, seven haplotypes carrying kdr mutations, namely 1014F1, 1014F2, 1014S1, 1014S2, 1014S3, 1014S4 and 1014C1 respectively were identified (Table 1), in the frequencies ranging from 2.7 to 80.0 % in the nine populations (Fig. 2). In the four populations in the northeast (Guigang, Guilin, Hezhou and Liuzhou), the frequencies of kdr mutations were more than 50.0 %, predominated by 1014C1. In contrast, only 1014S in a low frequency (<5 %) was observed in the two western populations (Baise and Nanning). Three types of kdr mutations coexisted in the other three populations (Hechi, Wuzhou, Yilin) in a frequency from 11.6 to 33.3 %.

Genealogical analysis of kdr mutations
To estimate the evolutionary relationship of kdr mutations, sixteen haplotypes identified in this study and ten haplotypes retrieved from GenBank database (Table 1) were used to construct a network using TCS 2.1 and Network 5.0. The two software generated identical results for our dataset, thus only the network produced by Network 5.0 was displayed (Fig. 3).
The network analysis revealed complex reticulate pattern and multiple independent mutation events leading to kdr haplotypes. For example, haplotype 1014F1 was possibly derived from 1014L1, while 1014F2 from 1014L3; 1014S1, S2, S3 and S4 were perhaps the results of an independent mutational step from four different wild haplotypes 1014L5, 1014L2, 1014L3 and 1014L4 respectively.
The genealogical analysis also revealed that a single mutation (T165C) might result in the resistant phenotype from the susceptible one. For example, 1014S2 shared identical intron sequence with 1014L2, L5 and S1, and a single mutation (T165C) might change the susceptible 1014L2 to the resistant 1014S2.

Discussion
Sequence analysis reveals that the sodium channel gene of An. sinensis has diverse genetic mutations. From 313 mosquitoes collected from Guangxi, 16 haplotypes were identified (Tables 1, 2). The number and frequency of AS-VGSC haplotypes are different among the nine populations from Guangxi. Overall, the 1014L1, 1014L2 and 1014C1 rank the top three haplotypes. Some haplotypes, for example 1014L2, 1014L3 and 1014L8, are widespread, and some haplotypes such as 1014L7 are restricted to a particular location (Table 2). Less number of haplotypes detected in Hezhou and Guigang may be related to the relative small sample size and/or strong selective pressure leading to 1014C1 in a high frequency. The high polymorphism of AS-VGSC may be partly explained by their large population size, wide distribution range of An. sinensis [14], diverse natural landscapes and different local insecticide selective pressure.
Three non-synonymous mutations (TTT for 1014F, TGT for 1014C and TCG for 1014S) were identified in this work. However, the TTC-encoding 1014F, previously found in Korea [32] and in Anhui and Sichuan of China [19], was not detected in Guangxi. The presence of L1014W substitution was reported based on the direct sequencing data of only one An.  [17]. This result is questionable and needs reconfirmation, because 1014F/W is not the only possibility encoded by the individual genotype T(T/G)(T/G). In this and another study [33], sequencing of cloned fragment with the T(T/G)(T/G) template resolved into two alleles, TTG and TGT, encoding 1014L/C rather than 1014F/W. Differential kdr distribution patterns were observed in the nine examined An. sinensis populations from Guangxi. Overall, haplotypes carrying 1014C mutation are widespread and present in high frequencies in the northeast, while 1014S is rare and distributed in the west (Table 1; Fig. 2). Notably, four of the nine populations with kdr allele frequencies higher than 50 % were located in the northeast of Guangxi, while only 1014S in a frequency less than 5 % was detected in the two western populations. Geographically, the frequency of kdr mutations decreased towards south and west from northeast. This pattern is in keeping with previous observations that a high kdr mutation frequency was detected in the central China which is located northeast to Guangxi [19,34,40], and kdr alleles were less prevalent in Yunnan and Hainan which are in the west and south to Guangxi respectively [16,18,19,34,40].
Considering the fact that haplotype 1014C1 and 1014F1, detected in samples of this study, are also distributed in central China (Table 1), it is possible that migration of An. sinensis population from central China through Hunan and Guizhou, may contribute to the occurrence of 1014C1 and 1014F1 in the northeast of Guangxi. The obvious distinction in the geographic distribution of each allele between populations in the west and in the northeast is likely a combined consequence of Table 1 Haplotypes identified or/and used in this study a Distribution information is adopted from cited references and this study. The polymorphic sites within intron 19   independent mutational events in different geographic locations and the geographical barriers limiting gene flow imposed by the mountainous landscapes of Guangxi. Theoretically, selective pressure may play a role in shaping the frequency of insecticide resistance-associated mutation in a population. However, it was not possible to assess the contribution of local selection of insecticides to the existing kdr distribution pattern because application history of insecticides was not known for samples used in this study. Previous studies have documented that the predominant kdr allele is L1014 F [4,12,13,16,19,[33][34][35]40], the high frequency of 1014C1 in the five An. sinensis populations in northeast of Guangxi (Fig. 3) is unexpected. 1014C was present in An. sinensis samples collected in Anhui, Guizhou, Henan, Hubei, Hunan, Jiangsu and Shandong of China [16,19,34,40] and in Korea [33], and was also found in VGSC of other mosquito species, such as Culex pipiens pipiens from China [36] and Anopheles albimanus from Mexico and Nicaragu [37]. These observations indicate that 1014C is a conserved mutation in mosquitoes and widely distributed. Why 1014C1 is prevalent in the five populations remains to be investigated. One possibility is that 1014C1 may provide a better protection, or/and a lower fitness cost, than other kdr alleles such as the classical 1014F.
The geographic distribution pattern and the genealogical analysis of kdr haplotypes strongly suggests that kdr mutations are not singly originated. For example, there are at least two independent mutation events giving rise to 1014F or 1014S haplotypes from a wild haplotype through a single mutational step (Table 1; Fig. 3). Multiple origins of resistance alleles via point mutations at the voltage-gated sodium channel gene have also been characterized in other insect species [38,39]. Interestingly, 1014C1 and 1014F1 is co-distributed in Guangxi and share the same intron (Tables 1, 2). The network analysis indicates that only one mutational step is able to change 1014F1 to 1014C1 (Fig. 3). Based on these observations, it is hypothesized that 1014C1 may represent a kdr allele evolved via a sequential progression from 1014F1 and perhaps has gradually replaced 1014F1.

Conclusions
Sixteen AS-VGSC haplotypes were identified in An. sinensis from Guangxi, suggesting diverse genetic mutations in the AS-VGSC gene. Three types of kdr mutations (1014F, 1014S and 1014C) were detected in Guangxi. Genealogical analysis suggests multiple origins of kdr mutations in An. sinensis. Although a direct link between kdr mutations and level of insecticide resistance has not been established due to the unavailability of susceptibility data for the samples used in this study, the high frequency of kdr mutation (>50 %) in Guigang, Guilin, Hezhou and Liuzhou indicates a risk of resistance to pyrethroids and DDT in An. sinensis in these areas. The  heterogeneities in the geographic distribution and multiple origins of kdr alleles highlight the need of a locationcustomized strategy for monitoring and management of insecticide resistance.

Fig. 3
Network showing the genealogy of AS-VGSC haplotypes. Empty circle means the wildtype haplotype (without 1014 mutation), solid circles mean mutant haplotypes with 1014 mutation (black 1014F, grey 1014S, dark green 1014C). The size of each circle is proportional to the corresponding frequency of a certain haplotype identified in Guangxi province. Straight line indicates the possible mutational step. The note above the line referred to the mutation position and base. Haplotypes marked with Asterisk mean that these haplotypes are retrieved from GenBank. H represents the type of intron haplotypes (see Table 1 for details)