- Open Access
Identification and phylogenetic analysis of voltage-gated sodium channel haplotypes in the malaria vector Anopheles sinensis using a high-throughput amplicon sequencing approach
Parasites & Vectors volume 14, Article number: 499 (2021)
Anopheles sinensis is a dominant vector for malaria transmission in Asian countries. Voltage-gated sodium channel (VGSC) mutation-mediated knock-down resistance (kdr) has developed in many A. sinensis populations because of intensive and long-term use of pyrethroids. Our previous study showed that multiple mutations at position 1014 of the VGSC were heterogeneously distributed in A. sinensis populations across Sichuan, China.
To understand resistance genotypes at the haplotype level and reconstruct the phylogenetic relationship of VGSC haplotypes, a cost-effective next-generation sequencing (NGS)-based amplicon sequencing approach was established to clarify haplotypes containing codon 1014 of the VGSC gene from a total of 446 adults collected in 12 locations of Sichuan, China.
Nineteen (19) haplotypes were identified, including 11 wild 1014L, 6 resistance 1014F, and 2 resistance 1014C haplotypes. We found that resistance haplotypes of A. sinensis VGSC were widely distributed at frequencies ranging from 3.67 to 92.61%. The frequencies of the 1014C haplotype in the southeast of Sichuan (Luzhou, Guangan, and Suining) were relatively higher than those in other sampling locations. Phylogenetic analyses support that kdr-type mutation at position 1014 is not singly originated and resistance 1014C haplotypes evolve from TTT-encoding 1014F.
A cost-effective next-generation sequencing (NGS)-based amplicon sequencing approach has been established in this study. The data revealed the patchy distribution of VGSC resistance haplotypes with overall high frequencies in Sichuan, China. Phylogenetic analyses support multiple origins and sequential evolution (1014L → 1014F → 1014C) for kdr-type mutations in A. sinensis.
Chemical insecticides have been the primary weapon for preventing and controlling vector-borne diseases, including malaria, in the past decades . The World Health Organization (WHO)-recommended pyrethroids are the most commonly used insecticides in China for indoor residual spray (IRS), insecticide-treated bed nets (ITNs), and mosquito coils . Many studies have demonstrated that mutations at position 1014 of the VGSC (numbered after the Musca domestica para-type sodium channel protein, GenBank: CAA65448.1), such as 1014F/H, can reduce the sensitivity of arthropod VGSCs to pyrethroids and cause so-called knock-down resistance (kdr) [3,4,5,6,7,8,9].
Sichuan province is located in the southwest of China (Fig. 1). The natural environment in most parts of Sichuan is suitable for the breeding of malaria vectors. Historically, Sichuan has been a malaria-endemic region in China. Although no indigenous malaria case has been reported since 2011, the increasing imported cases indicate a potential risk for the re-emergence of malaria in Sichuan [10, 11]. The intensive use of insecticides such as deltamethrin and permethrin for mosquito-targeted control and in agriculture has led to the development of insecticide resistance in the primary malaria vector Anopheles sinensis in this region [12, 13].
A recent study detected multiple mutations at position 1014 of the VGSC in Sichuan A. sinensis by Sanger sequencing . However, due to the diverse polymorphisms in the VGSC gene, direct DNA sequencing of PCR product from heterozygous individuals is incapable of haplotype identification. The identification of haplotype is useful for understanding the origin and evolution of kdr-associated mutations and helpful for monitoring distribution and spread of kdr haplotypes. In this context, we adopted the NGS platform to establish an amplicon-based sequencing pipeline to clarify VGSC haplotypes. The specific objectives of this effort were (i) to identify VGSC haplotypes and delineate their geographical distribution in A. sinensis populations across Sichuan and (ii) to elucidate the evolutionary relationship of kdr-type resistance mutations.
Anopheles sinensis collections
A total of 446 adult A. sinensis were caught from 12 sampling locations across Sichuan from August to September in 2018. The sampling locations are shown in Fig. 1.
Anopheles sinensis genomic DNA (gDNA) preparation
The gDNA was isolated from each mosquito using the protocol of Rinkevinch et al. . The gDNA concentration for each mosquito was detected by Nanodrop 2000 (Thermo, USA) and then adjusted to 50 ng/μl with ddH2O. Then, the gDNA of individual mosquitoes from the same sampling location was pooled in equal volume (10 μl). In this manner, 12 pooled gDNA samples representing 12 populations were obtained and stored at – 20 ℃ until subsequent treatments. The procedure is shown schematically in Fig. 2 (I–III).
VGSC amplicon preparation and sequencing
Primers with population-specific tags for amplicon preparation were designed based on the primers used to amplify a 325-bp genomic region that contains codon 1014 and a complete intron of the VGSC (kdr-F: 5′-TGCCACTCCGTGTGTTTAGA; kdr-R: 5′-GAGCGATGATGATCCGAAAT) . The population-specific tagged sequence (4–6 bp) was added at the 5′ end of the primer pairs. Details of the primers with tag sequences are shown in Table 1. The PCR mixture (50 μl) consisted of 1.5 μl gDNA template, 1 μl kdr-F, 1 μl kdr-R, 25 μl high-fidelity DNA polymerase (PrimeSTAR® max, Takara, China), and 21.5 μl ddH2O. The reaction program was 95 ℃ for 3 min, followed by 30 cycles each with 95 ℃ for 30 s, 55 ℃ for 30 s, 72 ℃ for 30 s, and a final extension of 5 min at 72 ℃. PCR products from each pooled gDNA sample were purified by agarose gel DNA purification kit (Accurate Biology, China) and then pooled with equal mass. DNA libraries were prepared using TruSeq DNA PCR-free kits (Illumina, San Diego, CA, USA). DNA sequencing was performed using the Hiseq2500 sequencer (Illumina, San Diego, CA, USA) with PE250 read type by Novogen Co. Ltd (Beijing, China). The procedure is shown schematically in Fig. 2 (IV–VII).
Two FASTQ format files with a total of 1.496 G raw data were obtained after sequencing. Low-quality sequences were filtered using Trimmomatic , and a total of 1.487 G clean data was obtained. Clean data (FASTQ format) were then submitted to pandaseq 2.11  for pair-end assembling. Subsequently, based on population-specific tagged-primer sequences, assembled sequences belonging to the respective populations were separated by the Linux opera system's grep command (Ubuntu 20.04 LTS). Sequence reverse complement was performed by Seqtk (https://github.com/lh3/seqtk). For chimera sequence filtering, sequences with both 5′-end and 3′-end matched population-specific tagged-primer were passed. The deletion of pair-end tag sequences was performed by the Linux opera system's sed command. After that, identical sequences from each population were merged and counted using FASTX-Toolkit (http://hannonlab.cshl.edu/fastx_toolkit/). The haplotype with a count of less than the threshold was removed (threshold = S/2 N; S is the sum of the number of clean sequences in each population, and N is the sample size from each sampling location). The brief pipeline is presented in Fig. 2 (VIII–XII).
Referring to the method of Délye et al. , Pearson’s correlation analysis was used for assessing the agreement in allele frequencies obtained by Sanger sequencing  and the NGS approaches. The Pearson’s correlation coefficient (r) was calculated by function cor () in the R programming language. The haplotypes identified in this study and those retrieved from GenBank [20,21,22,23] were used for constructing a phylogenetic tree by the maximum likelihood (ML) method and Kimura two-parameter model  (1000 bootstrap replicates) in MEGA-X software . Homologous sequences of Anopheles culicifacies were used as outgroup. The genealogical network analysis of haplotypes was performed by Network 5.0 (https://www.fluxus-engineering.com/sharenet.htm) using median joining for network calculation .
Identification of VGSC haplotypes
Fourteen nucleotide polymorphic sites (eight in intron and six in exon) were observed in the amplified fragments (325 bp) of the A. sinensis VGSC gene (Fig. 3). The mutations at nucleotides 203 and 204 led to non-conservative amino-acid residue substitutions (L to F or L to C) at position 1014. In addition, a previously reported mutation at locus 178 (A/T) [20, 22], which causes a conservative substitution (L1006I), was also detected. Based on all the DNA sequences, a total of 19 haplotypes were identified, including 11 wild-type 1014L, 6 kdr-type 1014F, and 2 kdr-type 1014C. Five haplotypes (H8-L, H14-L, H19-L, H11-F, and H18-C) were newly identified in this study (Table 2).
Frequency and distribution of each VGSC haplotype
The number of haplotypes in each sampling location ranged from 6 to 13 (Table 3). H2-L and H3-L were the two common wild VGSC haplotypes distributed in all locations except NJLC, while H12-L and H14-L were rare and only distributed in PZH and GAQF, respectively. The overall frequency of resistance haplotypes was high in Sichuan (Table 3 and Fig. 4). Notably, H1-F was the predominant resistance haplotype across Sichuan except PZH (Fig. 4). H7-C had higher frequencies in LZJY, GAQF, and SNDY than in other sampling locations and was detected in 10 of the 12 locations, while H18-C was only distributed in LZJY at a low frequency (Table 3, Fig. 4).
Phylogenetic analysis of VGSC haplotypes
The maximum likelihood tree (Fig. 5) revealed that VGSC haplotypes could be clustered into two clades and both clades contained resistance haplotypes. Median-joining (MJ) network analysis (Fig. 6) showed that the 19 haplotypes identified in Sichuan could be divided into two groups. Group 1 included 14 haplotypes (H3-L, H4-F, H5-L, H6-L, H8-L, H9-L, H10-L, H11-F, H13-F, H14-L, H16-L, H17-F, H18-C, and H19-L), forming a complex reticulate network. H3-L took a center position in Group 1 and could be evolved to seven haplotypes (H5-L, H6-L, H8-L, H10-L, H16-L, H11-F, and H17-F) by one mutational step. H18-C was located at the far end of Group 1 and could be produced by one mutational step from the resistance haplotype H4-F (Fig. 6). Group 2 included five haplotypes (H2L, H12L, H1F, H15F, and H7C) with a simple evolutionary relationship. The resistance haplotype H1-F could be derived from the wild H2-L and produce two other resistance haplotypes (H15-F or H7-C) by one mutational step (Fig. 6).
In this study, the NGS platform was deployed to understand the type and frequency of kdr mutation at the haplotype level. Correlation analysis showed that the frequency of kdr alleles detected by the NGS-based method herein was in good agreement with that determined by the Sanger sequencing . Pearson’s correlation coefficients, i.e. 0.912 (p-value = 3.54 × 10–5), 0.906 (p-value = 4.84 × 10–5), and 0.954 (p-value = 1.49 × 10–6) for the frequency of 1014L, 1014F and 1014C, respectively (Fig. 7), suggest that the NGS-based method is reliable. By comparison, NGS-based amplicon sequencing is more cost-effective than Sanger sequencing  and can clarify other nucleotide variations in the sequence of interest for haplotype identification.
Using this NGS approach, 19 VGSC haplotypes of A. sinensis were identified in Sichuan. The eight resistance haplotypes display a patchy distribution, and each sampling location has at least one resistance haplotype (Fig. 4). Notably, H1-F is distributed in all sampling locations except PZH and has an extremely high frequency (77.2%) in NJLC, and it also is the predominant resistance haplotype in other locations (Table 3 and Fig. 4). This observation led us to postulate that H1-F may have better protection and/or lower fitness cost than other kdr-type haplotypes in Sichuan. Special attention thus should be paid to this haplotype in future surveillance.
Interestingly, beyond Sichuan, both H1-F and H7-C have been previously documented in Zhejiang, Henan, Anhui, Guangxi, Guizhou, Hubei, and Jiangsu provinces of China [20, 23, 28,29,30] (see Fig. 1 for location). The wide distribution of H1-F and H7-C in locations that are hundreds of miles apart may suggest possible gene flows between these A. sinensis populations, though we could not exclude the possibility of independent origins of the same haplotype in different regions.
The topology of the ML tree (Fig. 5) and MJ network (Fig. 6) strongly suggests that the resistance mutations are not singly originated. The haplotypes carrying a kdr mutation can be convincingly clustered into two major clades or groups. Specifically, the complex reticulate network in Group 1 supports that divergent evolution may have occurred first, followed by diversification of 1014F haplotypes driven by selection pressure (H3-L → H11-F/H17-F; H6-L → H4-F/H13-F), and eventually converged toward H4-F. By comparison, the evolutionary relationships for haplotypes in Group 2 are relatively straightforward and support multiple 1014F haplotype origins (H12-L → KF927156-F; H2-L → H1-F). Multiple origins of kdr mutations have also been documented in other insect species such as the house fly Musca domestica .
Multiple kdr mutations at position 1014 of the VGSC have been reported in several insects . In A. sinensis, three types of kdr mutations (L1014F/C/S) have been documented [14, 23, 30]. The same mutations have been detected in a series of other anophelines [33, 34]. After reviewing the available literature, we found that the classical L1014F appears to be the first selected kdr allele, probably a consequence of DDT application. With global use of pyrethroids, new alternative mutations may emerge . Although the level of pyrethroid resistance conferred by L1014F and L1014H replacement in the house fly is clear , to our knowledge, there are no published data about the effect caused by different mutations at position 1014 on insecticide resistance in A. sinensis. Knowledge about insecticide tolerance and fitness associated with different resistance allele will be of importance to understand the evolution of kdr mutations and helpful for developing pest control strategies.
Our network analysis indicated that only one mutational step is able to change TTT-encoding 1014F to 1014C. This observation confirms a sequential evolution from 1014F to 1014C in A. sinensis as proposed in our previous paper . The codon usage may explain this notion: one mutational step can lead to the L1014F (TTG to TTT) or F1014C (TTT to TGT) substitution, while it requires two base mutations to realize the L1014C (TTG to TGT) replacement. However, this is not the case for the two resistance mutations (L1014F and L1014H) in house flies: the 1014F (TTT) allele could not evolve into the 1014H (CAT) allele (or vice versa) , because it would require two nucleotide substitutions.
This study adopted the NGS platform for cost-effectively clarifying VGSC haplotypes of A. sinensis and provided information about the types and distribution of the kdr-type mutations at the haplotype level in Sichuan, China. Our results support multiple origins of kdr alleles and confirm a scenario of sequential evolution from TTT-encoding 1014F to 1014C in A. sinensis. These findings will advance our understanding about the evolutionary history of VGSC resistance haplotypes in the malaria vector A. sinensis.
Availability of data and materials
All datasets are presented in this article.
Voltage-gated sodium channel
- kdr :
World Health Organization
Indoor residual spray
Insecticide-treated bed nets
Enayati A, Hemingway J. Malaria management: past, present, and future. Annu Rev Entomol. 2010;55:569–91.
Lu B. Integrated mosquito management. 2nd ed. Beijing: Science Press; 1999.
Smith TJ, Lee SH, Ingles PJ, Knipple DC, Soderlund DM. The L1014F point mutation in the house fly Vssc1 sodium channel confers knockdown resistance to pyrethroids. Insect Biochem Mol Biol. 1997;27:807–12.
Tan J, Liu Z, Wang R, Huang ZY, Chen AC, Gurevitz M, et al. Identification of amino acid residues in the insect sodium channel critical for pyrethroid binding. Mol Pharmacol. 2005;67:513–22.
Burton MJ, Mellor IR, Duce IR, Davies TG, Field LM, Williamson MS. Differential resistance of insect sodium channels with kdr mutations to deltamethrin, permethrin and DDT. Insect Biochem Mol Biol. 2011;41:723–32.
Du Y, Nomura Y, Satar G, Hu Z, Nauen R, He SY, et al. Molecular evidence for dual pyrethroid-receptor sites on a mosquito sodium channel. Proc Natl Acad Sci U S A. 2013;110:11785–90.
Usherwood PN, Vais H, Khambay BP, Davies TG, Williamson MS. Sensitivity of the Drosophila para sodium channel to DDT is not lowered by the super-kdr mutation M918T on the IIS4-S5 linker that profoundly reduces sensitivity to permethrin and deltamethrin. FEBS Lett. 2005;579:6317–25.
Vais H, Williamson MS, Goodson SJ, Devonshire AL, Warmke JW, Usherwood PN, et al. Activation of Drosophila sodium channels promotes modification by deltamethrin: reductions in affinity caused by knock-down resistance mutations. J Gen Physiol. 2000;115:305–18.
Scott JG. Life and death at the voltage-sensitive sodium channel: evolution in response to insecticide use. Annu Rev Entomol. 2019;64:243–57.
Xu GJ, Li L, Yu T, Wu XH. Epidemic situation and progress in eliminating malaria in Sichuan. J Prev Med Inf. 2014;30:783–5.
Li L, Liu Y, Xu GJ, Yu T, Zou Y, Wu XH, et al. Epidemiological characteristics of imported malaria in Sichuan Province in 2015. Chin J Parasitol Para Dis. 2017;35:75–85.
Cheng H, Yang W, Kang W, Liu C. Large-scale spraying of bednets to control mosquito vectors and malaria in Sichuan. China Bull World Health Organ. 1995;73:321–8.
Zhang S, Guo S, Feng X, Afelt A, Frutos R, Zhou S, et al. Anopheles vectors in mainland China while approaching malaria elimination. Trends Parasitol. 2017;33:889–900.
Qian W, Liu N, Yang Y, Liu J, He J, Chen Z, et al. A survey of insecticide resistance-conferring mutations in multiple targets in Anopheles sinensis populations across Sichuan. China Parasit Vectors. 2021;14:169.
Rinkevich FD, Zhang L, Hamm RL, Brady SG, Lazzaro BP, Scott JG. Frequencies of the pyrethroid resistance alleles of Vssc1 and CYP6D1 in house flies from the eastern United States. Insect Mol Biol. 2006;15:157–67.
Zhong D, Chang X, Zhou G, He Z, Fu F, Yan Z, et al. Relationship between knockdown resistance, metabolic detoxification and organismal resistance to pyrethroids in Anopheles sinensis. PLoS ONE. 2013;8: e55475.
Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114–20.
Masella AP, Bartram AK, Truszkowski JM, Brown DG, Neufeld JD. PANDAseq: paired-end assembler for illumina sequences. BMC Bioinform. 2012;13:31.
Délye C, Michel S, Pernin F, Gautier V, Gislard M, Poncet C, et al. Harnessing the power of next-generation sequencing technologies to the purpose of high-throughput pesticide resistance diagnosis. Pest Manag Sci. 2020;76:543–52.
Chen S, Qin Q, Zhong D, Fang X, He H, Wang L, et al. Insecticide resistance status and mechanisms of Anopheles sinensis (Diptera: Culicidae) in Wenzhou, an important coastal port city in China. J Med Entomol. 2019;56:803–10.
Qin Q, Li Y, Zhong D, Zhou N, Chang X, Li C, et al. Insecticide resistance of Anopheles sinensis and An. vagus in Hainan Island, a malaria-endemic area of China. Parasit Vectors. 2014;7:92.
Chang X, Zhong D, Fang Q, Hartsel J, Zhou G, Shi L, et al. 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: e2889.
Yang C, Feng X, Huang Z, Li M, Qiu X. Diversity and frequency of kdr mutations within Anopheles sinensis populations from Guangxi. China Malar J. 2016;15:411.
Kimura M. A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences. J Mol Evol. 1980;16:111–20.
Kumar S, Stecher G, Li M, Knyaz C, Tamura K. MEGA X: molecular evolutionary genetics analysis across computing platforms. Mol Biol Evol. 2018;35:1547–9.
Bandelt HJ, Forster P, Röhl A. Median-joining networks for inferring intraspecific phylogenies. Mol Biol Evol. 1999;16:37–48.
Sanger F, Nicklen S, Coulson AR. DNA sequencing with chain-terminating inhibitors. Proc Natl Acad Sci USA. 1977;74:5463–7.
Chang X, Zhong D, Lo E, Fang Q, Bonizzoni M, Wang X, et al. Landscape genetic structure and evolutionary genetics of insecticide resistance gene mutations in Anopheles sinensis. Parasit Vectors. 2016;9:228.
Zhang HW, Liu Y, Hu T, Zhou RM, Chen JS, Qian D, et al. Knockdown resistance of Anopheles sinensis in Henan province. China Malar J. 2015;14:137.
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. Malar J. 2015;14:120.
Rinkevich FD, Hedtke SM, Leichter CA, Harris SA, Su C, Brady SG, et al. Multiple origins of kdr-type resistance in the house fly, Musca domestica. PLoS ONE. 2012;7: e52761.
Dong K, Du Y, Rinkevich F, Nomura Y, Xu P, Wang L, et al. Molecular biology of insect sodium channels and pyrethroid resistance. Insect Biochem Mol Biol. 2014;50:1–17.
Lol JC, Castellanos ME, Liebman KA, Lenhart A, Pennington PM, Padilla NR. Molecular evidence for historical presence of knock-down resistance in Anopheles albimanus, a key malaria vector in Latin America. Parasit Vectors. 2013;6:268.
Silva AP, Santos JM, Martins AJ. Mutations in the voltage-gated sodium channel gene of anophelines and their association with resistance to pyrethroids—a review. Parasit Vectors. 2014;7:450.
Sun H, Tong KP, Kasai S, Scott JG. Overcoming super-knock down resistance (super-kdr) mediated resistance: multi-halogenated benzyl pyrethroids are more toxic to super-kdr than kdr house flies. Insect Mol Biol. 2016;25:126–37.
The authors thank Yan Yang, Xiuye Yang, Jipu Zou, Jing Feng, Jinsong Li, Chunxue Le, Yanjun Zeng, Xiqiang Zeng, Tingting Sang, and Xiping Luo for assistance in mosquito collection and the reviewers for their helpful comments and suggestions.
This work was supported by grants from the National Major Science and Technology Project of China on Infectious Diseases (no. 2017ZX10303404) and the State Key Laboratory of Integrated Management of Pest Insects and Rodents (IPM1915). The funders had no role in the study design, data collection, analysis, decision to publish, or preparation of the manuscript.
Ethics approval and consent to participate
Consent for publication
We declare we have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
About this article
Cite this article
Ni, R., Liu, N., Li, M. et al. Identification and phylogenetic analysis of voltage-gated sodium channel haplotypes in the malaria vector Anopheles sinensis using a high-throughput amplicon sequencing approach. Parasites Vectors 14, 499 (2021). https://doi.org/10.1186/s13071-021-05009-5
- Anopheles sinensis
- Next-generation sequencing
- Voltage-gated sodium channel
- Knock-down resistance
- Phylogenetic analysis