|Home | About | Journals | Submit | Contact Us | Français|
Anopheles punctulatus sibling species (An. punctulatus s.s., Anopheles koliensis, and Anopheles farauti species complex [eight cryptic species]) are principal vectors of malaria and filariasis in the Southwest Pacific. Given significant effort to reduce malaria and filariasis transmission through insecticide-treated net distribution in the region, effective strategies to monitor evolution of insecticide resistance among An. punctulatus sibling species is essential. Mutations in the voltage-gated sodium channel (VGSC) gene have been associated with knock-down resistance (kdr) to pyrethroids and DDT in malarious regions. By examining VGSC sequence polymorphism we developed a multiplex assay to differentiate wild-type versus kdr alleles and query intron-based polymorphisms that enable simultaneous species identification. A survey including mosquitoes from seven Papua New Guinea Provinces detected no kdr alleles in any An. punctulatus species. Absence of VGSC sequence introgression between species and evidence of geographic separation within species suggests that kdr must be monitored in each An. punctulatus species independently.
Malaria vectors in the Southwest Pacific region belong to the Anopheles punctulatus group. The species group range extends from Moluccas to Vanuatu, including New Guinea and islands of the Bismarck Archipelago, the Solomon Islands, and northern Australia.1 Thirteen sibling species (An. punctulatus sensu stricto [s.s.], An. species near punctulatus, morphologically indistinguishable Anopheles farauti 1–8 [Farauti complex; former An. farauti 1, 2, 3, and 7 now An. farauti s.s., Anopheles hinesorum, Anopheles torresiensis, and Anopheles irenicus, respectively], Anopheles koliensis, Anopheles clowi, and Anopheles renellensis)1–5 are known to comprise the species group. Blood feeding behaviors and breeding habitats of the Punctulatus group members have been characterized as unspecialized.6 Distribution of individual sibling species is heterogeneous, displaying unique and occasional overlapping patterns of dispersal.6 In Papua New Guinea (PNG), studies have specifically implicated An. farauti s.s., An. hinesorum, An. farauti 4, An. koliensis, and An. punctulatus s.s. as the primary vectors of malaria and filariasis.7–12 As evidence accumulates to verify the reproductive independence of these sibling species, it becomes increasingly important to collect evidence regarding their species-specific susceptibility to vector control strategies and competence as malaria vectors.
The importance of malaria transmission by An. punctulatus group mosquitoes was first recognized by Heydon and others in the 1920s.13,14 Despite establishing this connection over 90 years ago, mosquito and therefore malaria control in New Guinea has been inconsistently supported and pursued. As malaria casualties during WWII consistently outnumbered those from combat, military activities in New Guinea prompted intense United States and Australian interest in malaria control from 1942 to 1945.15–17 During this time troop protection included military discipline, efforts of malaria control and survey units beginning in 1943, atabrine prophylaxis and treatment of closed environments with “bug bomb” (pyrethrum).15–17 Perhaps the most effective means for controlling mosquito populations, DDT (dichloro-diphenyl-trichloroethane), was not available to troops in the Southwest Pacific area until late 1944, when military operations in New Guinea were drawing to a close.
DDT has been of greatest relevance to mosquito population management in New Guinea because it changed World Health Organization (WHO) strategies from malaria control to eradication.18,19 During the 1940–50s DDT mixed in oil-based solvents and sprayed on larval breeding sites showed short-term effectiveness against malaria transmission.20–22 However, DDT indoor residual spraying (IRS) was associated with reduced exposure to An. punctulatus group mosquitoes in New Guinea23–26 and remained active for 6 months after application. These combined features suggested that DDT IRS would enable deployment of time-efficient “attack” phase strategies by limited numbers of disease control specialists and unskilled workers.19,25,27 These findings encouraged DDT IRS activities from the late 1950s to 1970s,27–30 to cover ~50% of the population in 1972.29 Limited WHO insecticide susceptibility tests performed in the early 1970s showed high levels of DDT susceptibility among An. punctulatus group species in surveyed parts of PNG and the Solomon Islands.28,31 In 1972 the WHO Malaria Eradication Program in PNG ceased because of operational failure,32 however DDT IRS remained an important control strategy through provincial government activities.33,34 Further attempts at National Program IRS coverage in the late 1970s appeared to reduce anopheline mosquito and Plasmodium falciparum prevalence, but by 1983 insufficient resources and strained community relations brought this effort to a close.35–37
In the 1980s, the PNG Institute of Medical Research and Ministry of Health began to assess the effect of low-cost bednets. Although the studies were conducted in two different endemic sites north of the PNG Central ranges, results showed that both untreated and permethrin-treated nets (Wosera38,39 and Madang,40 respectively) were associated with protection from malaria exposure and prevalence. Additionally, studies showed that untreated bednets were associated with reduced transmission of lymphatic filariasis on Bagabag Island north of Madang.41
Presently, the Global Fund to Fight AIDS, Tuberculosis and Malaria (GFATM) has provided support for the distribution of over 2.5 million long-lasting insecticide-treated nets (LLIN; deltamethrin) at no cost to individual families.42 Introduction of LLINs to all of PNG's 20 Provinces represents the first attempt at a comprehensive plan for mosquito-based malaria control. Of importance, our recently published study observed no phenotypic evidence of resistance to the synthetic pyrethroid insecticides among Anopheles mosquitoes collected at five different PNG malaria-endemic sites.43
In response to intensive insecticide exposure, malaria vectors in other regions of the world have developed resistance to a number of different classes of insecticide compounds. Because several insect species have shown evidence for development of specific polymorphisms in the voltage-gated sodium channel (VGSC) after prolonged exposure to DDT and pyrethroid insecticides,44,45 a number of studies have focused on sequence encoding the S4–S6 region of VGSC domain II.46,47 In Anopheles mosquitoes, knockdown resistance (kdr) has been associated with an A → T mutation (L1014F)46 or a T → C mutation (L1014S)48 in VGSC. These single nucleotide polymorphisms (SNPs) have been found in African, Indian, and Asian malaria vectors,49–54 and most recently in neighboring Indonesia,55 showing the capacity of these mutations to develop independently in different Anopheles species. Another point mutation in VGSC associated with very high levels of DDT resistance, super-kdr (M918T), has been described in many insects, however this point mutation has not yet been detected in Anopheles.44,45,56
To date, Southwest Pacific malaria vectors have not been surveyed for any mutations associated with insecticide resistance. Assuming validity of An. punctulatus sibling species definitions, and the heterogenous history of insecticide application in PNG, we hypothesize that development, prevalence, and distribution of these mutations would occur independently among members of the An. punctulatus group. Here, we have developed a multiplex DNA-based assay to detect knockdown resistance (kdr) polymorphisms and simultaneously identify individual species in the An. punctulatus group by analysis of sequence polymorphisms in the VGSC gene. Results from our study make it possible to evaluate kdr mutations in a species-specific context by a single assay and provide insight into the evolution of VGSC sequence in the An. punctulatus group. Although additional factors can contribute to insecticide resistance and avoidance, population-based studies from other malaria-endemic regions and results presented here may provide important background data on VGSC polymorphism preceding significant selective pressure of insecticides in PNG.
The Entomology Unit of Papua New Guinea Institute of Medical Research (PNGIMR) collected mosquitoes throughout PNG as part of surveillance studies associated with distribution of LLINs and insecticide susceptibility trials.43 Mosquitoes were collected through landing catch methods, larva collections, and Center for Disease Control light traps as previously described.8,57 Collection sites span 24 villages from seven provinces in PNG and all locations were defined by the global positioning system (Table 1). A subset of these mosquitoes (n = 62) were collected as larvae, and then females reared to adulthood were evaluated by WHO insecticide susceptibility assays by Keven and others.43
Morphological identification of mosquitoes was completed using methods previously described where members of the An. punctulatus group were classified as An. koliensis, An. farauti s.l. (sensu lato) or An. punctulatus s.l. by morphologically distinguishing characteristics.10,58 Collected larvae were allowed to mature in the laboratory to facilitate adult morphological identification.43 Mosquitoes morphologically identified as members of the An. punctulatus species group (n = 357) were kept in coded vials containing silica gel until DNA extraction could be completed.
Genomic DNA was extracted from single whole mosquitoes and molecular species identification was determined by analyzing polymorphisms within the internal transcribed spacer 2 (ITS2) unit of ribosomal DNA (rDNA) using a high-throughput molecular probe method previously described.59
Polymerase chain reaction (PCR) amplification of the IIS4-IIS6 region of the VGSC including introns preceding and following the coding region containing the kdr allele (referred to as intron 1 and intron 2, respectively) was modified from Martinez-Torres and others60 molecular characterization of VGSC in An. gambiae. The PCR amplification reactions (25 μL) were preformed (67 mM Tris-HCl, pH 8.8, 6.7 mM MgSO4, 16.6 mM (NH4)2SO4, 10 mM 2-mercaptoethanol, 100 μM dATP, dGTP, dCTP, and dTTP, and 2.5 units of thermostable DNA polymerase) using upstream and downstream primers (aDip1 5′-TGGCCSACRCTGAATTTACTC-3′212133 and cDip2 5′-TTKGACAAAAGCAAGGCTAAGAA-3′). The thermocycling program included: 95°C 2 min (1×), 95°C 30 sec, 60°C 30 sec, 72°C 1 min 45 sec (40×), 72°C 4 min (1×) performed using a BioRad DNA Engine Tetrad 2 Peltier Thermal Cycler (Hercules, CA). To evaluate overall amplification efficiency, PCR products were separated by electrophoresis on 2% agarose gels (1× TBE), stained with SYBR Gold (Invitrogen, Carlsbad, CA), and visualized on a Storm 860 using ImageQuant, 5.2 software (Molecular Dynamics, Sunnyvale, CA). The VGSC PCR yielded a product ranging 1,304–1,380 bp.
TOPO TA cloning (Invitrogen, Grand Island, NY) of the amplified VGSC region was preformed for 91 individual mosquitoes (An. punctulatus s.s. n = 31, An. koliensis n = 20, An. farauti s.s. n = 14, An. hinesorum n = 13, An. farauti 4 n = 12, and Anopheles longirostris, a non-An. punctulatus group member, n = 1) followed by single-pass bidirectional plasmid sequencing (Beckman Coulter Genomics, Danvers, MA) (GenBank accession nos.: HQ173903–HQ173993). DNA sequence data were analyzed using Geneious 22.214.171.124 ClustalW262 (default settings: gap open penalty = 10, gap extend penalty = 6.66) was used to generate the VGSC species consensus sequence alignment as shown in Figure 1.
A point mutation was introduced in a VGSC PCR representative creating a positive control possessing the kdr mutation (TTT), using the Stratagene QuikChange Site-Directed Mutagenesis Kit (Agilent Technologies, La Jolla, CA) in accordance with the manufacture's protocol. These samples were then sequenced confirming that the introduction of the kdr point mutation (TTT) had occurred. Additionally, 100 bp oligonucleotides were synthesized to serve as controls, harboring the TTT, TCA, or TTA genotype (Integrated DNA Technologies, Coralville, IA).
A post-PCR ligase detection reaction-fluorescent microsphere assay (LDR-FMA) was designed based upon VGSC sequence analysis. The LDR probes were designed to detect the presence (TTT or TCA) or absence (TTA) of the L1014F/S mutations associated with knockdown resistance (Table 2). The LDR probes were also designed to detect polymorphic sequence motifs within the VGSC amplified region that appeared to segregate in patterns consistent with An. punctulatus species group designations (Table 2, Figure 1). These LDR probes were predicted to enable the differentiation of An. punctulatus s.s., An. koliensis, An. farauti s.s., An. hinesorum, and An. farauti 4.
Following amplification of the VGSC target site, PCR products were added to a multiplex LDR where sequence-specific upstream classification probes, specific to the three different kdr locus genotypes (TTT, TCA, or TTA), and species-specific differentiation probes (containing anti-TAG sequences specific to Luminex microsphere sets) ligate to downstream sequence reporter probes (3′-biotinylated) when appropriate target template sequences are available (Table 2). The LDRs (15 μL) were conducted in a solution containing 20 mM Tris-HCL buffer, pH 7.6, 25 mM potassium acetate, 10 mM magnesium acetate, 1 mM NAD+, 10 mM dithiothrietol, 0.1% Triton X-100, 10 nM (200 pmol) of each LDR probe, 1 μL of each PCR product and 2 units of Taq DNA ligase (New England Biolabs, Beverly, MA). The LDR reactions were initially heated at 95°C for 1 min, followed by 32 thermal cycles of 95°C for 15 sec and 58°C for 2 min. The labeling and detection of LDR products was completed as previously described63 using Bio-Plex array reader (Bio-Rad Laboratories).
Previous analysis of a small number of mitochondrial and nuclear genes has demonstrated that significant genetic polymorphism exists within and between members of the An. punctulatus species group species considered to be the most significant vectors of malaria in PNG.59,64–67 Additionally, some polymorphisms within species, namely An. farauti s.s. and An. koliensis, have been shown to be geographically localized.10,64,68 Given these previous findings, and because the VGSC for the An. punctulatus group has not yet been described, we preformed DNA sequence analysis across a segment of this gene (~1,365 bp) observed to include insecticide resistance polymorphisms in a number of insect species.47,69
The analysis was based on comparative alignment of 90 An. punctulatus group individual's VGSC sequences where mutations associated with knockdown (kdr) and super knockdown resistance (super-kdr) have previously been identified. This includes sequence encoding the fourth to sixth transmembrane regions of domain II; codons 908 to 1026, including two introns. The occurrence of kdr mutations was evaluated following alignment of 31 alleles of An. punctulatus s.s. (98.5% similar), 14 of An. farauti s.s. (99.4% similar), 13 of An. hinesorum (98.1% similar), 12 of An. farauti 4 (99.0% similar), and 20 of An. koliensis (98.4% similar), kdr (L1014F/S), and super-kdr (M918T). Among the surveyed sequences there was no evidence of mutations associated with insecticide resistance. Additional VGSC sequence comparisons were made between An. punctulatus species group members and similar sequence fragments from An. longirostris (GenBank accession no. HQ173903) co-endemic to PNG, An. dirus (DQ026439) Southeast Asia/Malaysia, and An. gambiae (Y13592) Africa.
Further examination of the coding region (358 bp) revealed 53 variable sites where 56 unique mutations were found. Among these, 20 SNPs were observed repeatedly within and between the An. punctulatus sibling species; 36 SNPs were observed only once. Although this error rate (1.08 × 10−3; 35 of 32,310 nucleotides) is higher than standard amplitaq error rates [2.0 × 10−5 to 2.1 × 10−4/bp]70,71), it was not possible to validate these latter SNPs through repeated observation in this study.
Within the VGSC coding region sequence, 10 SNPs occurred in multiple An. punctulatus sibling species (codons 909, 910, 924, 956, 973, 987, 996, 999, 1007, 1026; all occupying the third position of each codon [Table 3]). Codons 909 (G>C) and 910 (G>A) were found to harbor mutations that occurred together and separately in each species with the exception of An. farauti s.s. where the mutations were not observed together. The predominant 909-910 haplotype was GG occurring in 41 of 90 sequenced alleles (An. punctulatus s.s., n = 17; An. koliensis, n = 12; An. farauti 4, n = 6). In An. farauti s.s. CG was the predominant 909-910 haplotype (n = 7) and for An. hinesorum, the GA and CG SNP combinations were observed most often and equally (n = 4 for each). At codon 1026 a C>A SNP was observed in 31 of 90 sequenced alleles (An. punctulatus s.s. = 10, An. koliensis = 8, An. farauti s.s. = 5, An. Hinesorum = 5, An. farauti 4 = 2). Six of these repeatedly observed synonymous mutations were shared between two species. At codon 924 a G>A SNP was observed once in both An. farauti s.s. and An. hinesorum; at 956 a C>T transition was observed in An. hinesorum (n = 2) and An. punctulatus s.s. (n = 2); at 973 a T>C transition was observed in An. farauti s.s. (n = 1) and An. koliensis (n = 1); at 987 a T>C transition was observed at fixation in both An. farauti 4 (n = 12) and An. punctulatus s.s. (n = 31); at 996 C>T was observed at fixation in both An. farauti 4 (n = 12) and An. punctulatus s.s. (n = 31); 999 C>T was observed in An. farauti 4 (n = 1) and An. koliensis (n = 9), and A>C at 1007 was observed in An. farauti s.s. (n = 11) and An. punctulatus s.s. (n = 30). The remaining six repeatedly observed synonymous mutations observed in five codons (codons 962, 964, 980, 981, and 985) were species specific. At codon 962 a C>T SNP was observed in two An. farauti 4 alleles; at 964 an A>G SNP was observed at fixation in An. koliensis (n = 20); at 980 a T>G transversion was observed at fixation in An. farauti 4 (n = 12) while a T>C transition was observed at fixation in An. koliensis (n = 20); at 981 A>G was observed at fixation in An. farauti 4 (n = 12), and 985 A>G was observed at fixation in An. punctulatus s.s. (n = 31).
Coding sequence mutations were observed in separate codons in multiple species (codons 947, 967, 999, 1017). At codon 947 a T>C SNP was observed in one An. hinesorum allele and one An. punctulatus s.s. allele (codon position 1) and resulted in a phenylalanine to leucine amino acid substitution; at 967 an A>G SNP was observed in one An. farauti s.s. allele (codon position 2) and one An. punctulatus s.s. allele and resulted in an asparagine to serine substitution; at 999 a T>C SNP was observed in two An. farauti 4 alleles and one An. farauti s.s. allele (codon position 2) and resulted in a valine to alanine substitution; at 1017 a T>C SNP observed in one An. hinesorum, one An. koliensis, and one An. punctulatus s.s. allele (codon position 2) resulted in a leucine to proline substitution.
To augment the analysis of sequence variation among the An. punctulatus complex sibling species evaluated here, we performed additional VGSC sequence comparisons with similar coding region sequence fragments for An. longirostris (358 bp), An. dirus (213 bp), and An. gambiae (348 bp). Between An. longirostris and the 90 An. punctulatus sibling species alleles we observed five synonymous mutations unique to An. longirostris (954 T>C, 979 C>T, 981 A>C, 989 G>C, and 1012 C>A). Two synonymous mutations in An. longirostris were shared with at least one An. punctulatus species group allele (995 C>T An. punctulatus s.s. [n = 1]; 996 C>T An. punctulatus s.s. [n = 31], An. farauti 4 [n = 14]). One synonymous An. longirostris mutation occurred at the An. punctulatus s.s. SNP site 981 with a unique sequence change A>C compared with A>G. Between a partial An. dirus sequence and the An. punctulatus sibling species alleles we identified eight synonymous mutations unique to An. dirus (978 G>A, 989 G>C [also observed in An. longirostris], 993 T>C, 999 C>G, 1004 T>C, 1012 C>A [also observed in AL], 1013 C>T and 1014 T>C). Six synonymous mutations were observed to be shared among one or more An. punctulatus group members (956 C>T An. hinesorum [n = 2], An. punctulatus s.s. [n = 2]; 973 T>C An. punctulatus s.s. [n = 1], An. koliensis [n = 1]; 980 T>G An. farauti 4 [n = 12]; 981 A>G An. farauti 4 [n = 12], 996 C>T An. farauti 4 [n = 12], An. punctulatus s.s. [n = 31]; and 1026 C>A An. farauti 4 [n = 2], An. hinesorum [n = 5], An. farauti s.s. [n = 6], An. koliensis [n = 8], An. punctulatus s.s. [n = 10]). Finally, between a partial An. gambiae sequence and the An. punctulatus sibling species alleles we identified 19 unique synonymous mutations in An. gambiae (940 G>A, 955 T>G, 963 A>G, 969 G>A, 974 T>C, 979 T>C [also observed in An. longirostris], 982 G>A, 989 G>A, 992 C>T, 999 C>A, 1002 T>A, 1003 T>A, 1004 C>T, 1008 C>T, 1010 A>G, 1012 C>A [also observed in An. longirostris and An. dirus], 1013 C>T [also observed in An. dirus], 1015 G>C, and 1018 T>C). Nine synonymous SNPs were observed to be shared among one or more An. punctulatus group members (956 C>T [also observed in An. dirus] An. hinesorum [n = 2], An. punctulatus s.s. [n = 2]; 962 C>T An. farauti 4 [n = 2]; 981 A>G [also observed in An. dirus] An. farauti 4 [n = 12]; 983 T>C An. punctulatus s.s. [n = 1]; 995 C>T [also observed in An. longirostris] An. punctulatus s.s. [n = 1]; 997 C>T An. farauti 4 [n = 1]; 1000 T>C An. koliensis [n = 1]; 1007 A>C An. farauti s.s. [n = 11], An. punctulatus s.s. [n = 30]; and 1026 C>A [also observed in An. dirus] An. farauti 4 [n = 2], An. hinesorum [n = 5], An. farauti s.s. [n = 6], An. koliensis [n = 8], An. punctulatus s.s. [n = 10]).
Because the potential exists for this gene to be under selection and since there has been a history of insecticide exposure among these species it becomes necessary to look at the nature of the mutations in the coding regions of the VGSC. Previously, we described 20 SNPs (16 synonymous and 4 nonsynonymous) occurring in the coding region analyzed, which were found in more than one sequence. Seven synonymous mutations (964 A>G An. koliensis n = 20, 980 T>C An. koliensis n = 20, 980 T>G An. farauti 4 n = 12, 981 A>G An. farauti 4 n = 12, 985 A>G An. punctulatus s.s. n = 31, 987 T>C An. punctulatus s.s. n = 31, and An. farauti 4 n = 12) were fixed within at least one An. punctulatus sibling species. The remaining 9 synonymous mutations were found to be polymorphic in one or more of the sibling species (909 G>C An. farauti 4 n = 3, An. farauti s.s. n = 7, An. hinesorum n = 13, An. koliensis n = 4, An. punctulatus s.s. n = 31, 910 G>A An. farauti 4 n = 5, An. farauti s.s. n = 3, An. hinesorum n = 6, An. koliensis n = 6, An. punctulatus s.s. n = 8, 924 G>A An. farauti s.s. n = 1, An. hinesorum n = 1, 956C>T An. punctulatus s.s. n = 2, An. hinesorum n = 2, 962 C>T An. farauti 4 n = 2, 973 T>C An. farauti s.s. n = 1, An. farauti 4 n = 1, 999 C>T An. koliensis n = 9, An. farauti 4 n = 1, 1007 A>C An. punctulatus s.s. n = 30, An. farauti s.s. n = 11, 1026 C>A An. punctulatus s.s. n = 10, An. koliensis n = 8, An. farauti s.s. n = 6, An. hinesorum n = 5, An. farauti 4 n = 5). All 4 nonsynonymous mutations were found to be polymorphic in one or multiple sibling species (947 T>C An. punctulatus s.s. n = 1, An. hinesorum n = 1, 967 An. punctulatus s.s. n = 1, An. farauti s.s. n = 1, 999 T>C An. farauti s.s. n = 1, An. farauti 4 n = 2, 1017 T>C An. punctulatus s.s. n = 1, An. koliensis n = 1, An. hinesorum n = 1). The assumption of neutrality was tested by the McDonald-Kreitman test72 using this coding region data. The results were non-significant suggesting that there is no evidence for selection in this region of the gene among the species in this sample. Given that the focus of this analysis was only on a portion of the VGSC gene, it is possible that selection could be occurring elsewhere.
Intron 1 was of variable length within (exception An. farauti s.s.) and between species; An. punctulatus s.s. ranged from 880 to 923 bp, An. farauti s.s. 928 bp, An. hinesorum 939 to 950 bp, An. farauti 4 925 to 938 bp, and An. koliensis 940 to 953 bp. Intron 2 varied in size between species; An. punctulatus s.s. 65 bp, An. farauti s.s. 65 bp, An. hinesorum 68 bp, An. farauti 4 66 bp, and An. koliensis 68 bp. Further unique species-specific insertions and deletions were also observed within introns 1 and 2. These differences are shown in Figure 1 where consensus sequences from each An. punctulatus sibling species of interest were aligned using ClustalW. When the An. punctulatus sibling species consensus sequences were compared with An. longirostris, they were found to be 77–79% similar.
We then analyzed the intronic regions for each species independently (Table 4). For intron 1 (ranging 880–953 bp) all An. farauti s.s. sequences had the greatest pairwise percent identity, 99.5%. High percentage identity is not surprising as all 14 sequences were the same length, leaving differences to be SNPs alone. Pairwise percent identity for the remaining within species comparisons were as follows: An. punctulatus s.s. (98.1%), An. hinesorum (97.8%), An. farauti 4 (97.8%), and An. koliensis (98.1%). Pairwise percent identity was 88.5% when comparing all An. punctulatus complex members. When comparing all An. punctulatus sibling species with An. longirostris the pairwise percent identity dropped slightly to 87.8%, however when comparing An. longirostris individually to a representative of each species, pairwise percent identities were < 69%. Intra- and inter-species complex relationships are further defined by considering the percent of identical sites within each comparison. Intra-species comparisons observed nucleotide identity over the analyzed gene fragment of 84.2% for An. punctulatus s.s., 97.2% for An. farauti s.s., 91.3% for An. hinesorum, 95.5% for An. farauti 4, and 91.3% for An. koliensis. With inter-species comparisons among An. punctulatus sibling species allelic percent identity dropped to 51.7%; 40% when An. longirostris was added to this sequence comparison. This further demonstrates that intronic sequences appear to have been conserved within each species and among members of the An. punctulatus species group. Given that this gene is possibly under selective pressure and natural selection reduces the genetic variation of organisms we would expect to see little to no variation within this gene within species.
Unlike intron 1, intron 2 was not observed to vary in size within species and is much smaller (65–68 bp). Pairwise percent identity of all the species was found to be greater than 97.7% (An. hinesorum) with four species (An. punctulatus s.s., An. farauti s.s., An. farauti 4, and An. koliensis) above 99%. Comparison of all An. punctulatus species group members showed 90.8% pairwise percent identity and when An. longirostris was added to the comparison this percentage dropped slightly to 90.2%. Conserved intron 2 sequences further demonstrated species similarities through analysis of sequence identity. An. punctulatus s.s., An. farauti s.s., An. farauti 4, and An. koliensis all have greater than 96.9% identical sites, and An. hinesorum has 89.7% identical sites. When comparing all An. punctulatus s.s. group members together, there were only 60.6% identical sites; overall sequence identity dropped to 43.7% when An. longirostris was added to these comparisons.
Given that species-specific polymorphism in both the coding and intronic regions were observed in the multiple species alignment (Figure 1) in addition to the multiple SNPs shared among more than one species (Table 3), we preformed a phylogenetic comparison of the entire sequenced portion of VGSC (1,303–1,379 bp) of the 90 An. punctulatus species group members. Anopheles longirostris was used as an out-group. Using ClustalW (at default parameters), sequences were aligned and phylogenetic relationships were inferred using the Neighbor-Joining method in MEGA573 with 1,000 bootstrap replicates (Figure 2). Branches with bootstrap values < 80 were collapsed.
We observed that each of the An. punctulatus sibling species formed its own clade with bootstrap support greater than 97, and that no individual mosquito was found in a clade that was different from its initial species identification based on detailed comparisons of ITS2 rDNA sequence polymorphisms. Phylogenetic analyses also observed significant stratification between An. hinesorum mosquitoes collected north and south of the PNG Central Ranges. Overall, our sequence comparisons suggested that the VGSC gene sequence would be a reliable candidate for molecular species identification that could be coupled with probes investigating the kdr-associated SNPs.
Conserved regions of the amplified VGSC sequence surrounding both the location of the M918T and L1014F/S mutations afforded the opportunity to build ligase-detection (LDR) probes to monitor for the presence (heterozygous or homozygous) or absence of these mutations in a post-PCR-based LDR-FMA. Given the high bootstrap support for the ability of VGSC to differentiate members of the An. punctulatus group, we targeted species-specific regions of the amplified VGSC region in development of LDR probes to differentiate members of the An. punctulatus group. These probes, coupled with the kdr detection probes, could simultaneously monitor the presence or absence of the kdr mutations and differentiate members of the An. punctulatus species group. Locations where kdr mutation detection and species-specific classification LDR reporter probes hybridized with VGSC template PCR products are shown in the VGSC consensus alignment for An. punctulatus s.s., An. farauti s.s., An. hinesorum, An. farauti 4, and An. koliensis (Figure 1).
The specificity of the assay was demonstrated using probes identified in Table 2 to differentiate PCR products amplified from previously sequenced species-specific controls including two controls, known to harbor the sensitive (TTA) or resistant (TTT or TCA) kdr mutation at amino acid position 1014. As members of the An. punctulatus species group are known to carry Plasmodium9,11,12 and Wuchereria parasites10 as well as human blood meals,3,12 six additional controls were included containing P. falciparum, Plasmodium vivax, Plasmodium malariae, Plasmodium ovale, Wuchereria bancrofti, and human genomic DNA.
Results of control LDR-FMA tests (Figure 3) showed high specificity for kdr allele differentiation (kdr-sensitive [TTA] mean fluorescence intensity [MFI] > 23,000; kdr-resistant [TTT or TCA] > 15,000 MFI; background signals < 3,000 MFI) as each LDR-FMA probe set detected only the kdr allele present in the genomic controls and displayed. Simultaneously, each species-specific LDR-FMA probe set (An. punctulatus s.s. MFI > 21,000; An. koliensis > 23,000; An. farauti s.s. > 23,000; An. hinesorum > 10,500; An. farauti 4 > 22,000), detected only the species expected; all species probe background signals were below MFIs of 2,000. Results showed that the LDR probes designed to target the species-specific polymorphisms described here detected only the species present and no others. Fluorescent signals from the parasite and human samples were below a background of 160 MFI (data not shown), showing that the presence of parasites and human blood meal would not obscure Anopheles species identification or kdr mutation detection.
The VGSC multiplex assay was then used to evaluate a sample collection of 312 mosquitoes for the kdr mutation and species identification. All 312 mosquitoes assayed were determined to be homozygous for the sensitive kdr allele (TTA/TTA). Results in Table 5 summarize the comparison of species identification using the VGSC multiplex assay to a similar molecular classification method that probes differences in An. punctulatus sibling species ITS2 rDNA sequences.59 One discordant sample was detected among the 312 mosquitoes assayed; by VGSC LDR-FMA this mosquito was identified as An. punctulatus s.s., where ITS2 LDR-FMA species were identified as An. farauti s.s.. These results show 99.7% concordance between the two species identification methods while simultaneously identifying the kdr mutation status in a large sample set of field-collected mosquitoes.
Mosquitoes that were evaluated for insecticide susceptibility using WHO insecticide bioassays (as described previously in Keven and others43) were also analyzed in this assay. Sixty-two mosquitoes, collected from five geographically distinct areas (Table 1) were previously tested and reported to be susceptible to either 0.05% lambda-cyhalothrin-treated filter paper (n = 26, locations: Dimer, Peneng, and Lorengau) or new 55 mg/m2 deltamethrin-treated LLIN (n = 36, locations: Ramu [Madang, Naru, and Dimer). These mosquitoes were evaluated using the VGSC LDR-FMA. All 62 samples gave a strong signal for the presences of the susceptible genotype (TTA; average MFI 8,182; signals for TTT and TCA probes fell within expected background [< 3,000 MFI]) showing that these samples were both phenotypically susceptible and lacked the L1014F/S mutations associated with kdr.
Here, we have characterized DNA sequence polymorphism in a segment of the VGSC gene (~1,350 bp) among sibling species of the An. punctulatus group known to be significant vectors of human malaria and filariasis in Papua New Guinea.7–12 This region of the VGSC gene encodes the fourth to sixth transmembrane regions (domain II) of the VGSC protein,60 and is known to harbor mutations associated with knockdown resistance to pyrethroids and DDT in many insect species (Anopheles,48,52,53,55,60,74–79 Culex,80,81 Musca domestica [House fly],82,83 Blattella germanica [German cockroach]).83,84
In this characterization of VGSC genetic variation we were able to determine that the well-known kdr-associated mutations (L1014F/S; M918T82) were not observed among the extant mosquitoes of the An. punctulatus species group. Although our collections were not comprehensive across PNG, individual mosquitoes were available at multiple locations along the PNG north coast where transmission of malaria and filariasis is intense. Because many of these collection sites occurred in regions that have historically been the focus of DDT larviciding applications22,27 and IRS programs,25,29 it is possible that ancestral An. punctulatus species group populations have been exposed to positive selection that would favor the kdr-associated mutations. This included mosquito collection sites (Table 1) in six mainland Provinces and Manus Island Province. The greatest distance between sample collection sites within each species ranged from 85 km for An. farauti 4 to 883 km for An. farauti s.s., with all species (except for An. farauti 4) having collections spanning at least 580 km in distance. For this study An. farauti 4 collections were limited to pockets north of the Central Ranges. Significant geographical barriers also separated collection sites for An. farauti s.s. (island locations) and An. hinesorum (locations north and south of the Central Ranges). Interestingly, although our phylogenetic analysis provided evidence for significant population stratification within An. hinesorum collected north and south of the PNG Central Ranges, significant VGSC sequence variation was not observe between An. farauti s.s. mosquitoes collected between the PNG mainland and Manus Island.
Furthermore, the overall sequence analysis revealed a wide range of synonymous and nonsynonymous SNP and insertion/deletion polymorphisms, the latter exclusively within “so-called” introns 1 and 2. The additional coding region SNPs enabled us to perform population genetic comparisons to evaluate evidence for selective pressure on this gene sequence and sequence-based relationships within and between sibling species. McDonald-Kreitman testing,72 based on the observation of 16 synonymous versus 4 non-synonymous SNPs in the VGSC coding region analyzed here indicated that there was no evidence for selection in this gene segment. However, it may be important to note that these observations have been based on mosquitoes collected over the past 5–10 years during a time when no concerted insecticide dispersal program has been under implementation in PNG. Therefore, it is not possible to rule out the potential that the kdr-associated mutations arose in the past and have now drifted out of current An. punctulatus sibling species populations.
Given results of the McDonald-Kreitman test, reservations about using VGSC sequences to perform phylogenetic analyses may be reduced. The Neighbor-Joining analysis reported in Figure 2,73,85,86 shows very strong conservation of sequence within each species cluster and no evidence for gene flow between species. Because many of the polymorphisms characterizing the different species included consistently observed insertions and deletions, our results suggested that it might be possible to develop highly specific DNA probes for species identification. Consistent with these phylogenetic results (apart from one sample), DNA probes based on VGSC sequence polymorphisms produced species identification results consistent with earlier ITS2 analyses for over 300 mosquitoes collected across our widely distributed sampling sites; expanded surveys will be necessary to verify that these results can be consistently reproduced throughout the Southwest Pacific region where the An. punctulatus sibling species are endemic. Although our resulting LDR-FMA multiplex assay did not detect any evidence of kdr-associated mutations from any of the PNG study sites surveyed, our analyses revealed an important feature of VGSC sequence evolution and an important potential challenge for mosquito control. Because our results provide evidence for independent acquisition and maintenance of VGSC sequence variation, kdr-associated mutations in this gene sequence would also be expected to emerge independently within each species. Therefore, whether VGSC kdr-associated mutations or mutations in other genes lead to insecticide resistance, assays must be developed to monitor each of the An. punctulatus sibling species with equal efficiency, specificity, and sensitivity. Therefore, molecular diagnostic assay developed here, offers advantages over other PCR-based species identification strategies by querying insecticide resistance associated and species-specific polymorphisms in a single multiplex test. Currently, countrywide planning and distribution of LLINs is underway in the Solomon Islands34,87 and Papua New Guinea,42,43 respectively. A recent report describing the discovery of the VGSC L1014F kdr mutation in populations of Anopheles from neighboring Indonesia suggests that An. punctulatus sibling species in this same region55 may be experiencing insecticide exposure that will select for kdr-associated mutations.
Integrated vector management strategies require a means to continuously monitor and assess the vector species composition in an endemic region and elements of the vector's biology. Although kdr is a concern, it is important to understand this phenotype within the context of additional gene sequence and gene expression polymorphism, physiological and behavioral heterogeneity. Because insecticide resistance can develop through many different routes88 it is becoming necessary to broaden evaluation of insecticide resistance beyond candidate gene-based investigations as demonstrated by recent genomic comparison studies of An. gambiae.89,90
We are grateful to PNGIMR Entomology Unit for aid in mosquito collections and morphological identification. We thank Wallace Peters, David Serre, and Scott Small for helpful discussions and constructive criticism leading to completion of this manuscript.
Financial support: This study was supported by grants from National Institutes of Health (AI065717) and Fogarty International Center (TW007872, TW007377, and TW007735). Mosquito collection was also supported by the Global Fund to Fight Aids, Tuberculosis and Malaria Round 3 malaria grant to PNG.
Authors' addresses: Cara N. Henry-Halldin, Kogulan Nadesakumaran, Allison M. Zimmerman, James W. Kazura, and Peter A. Zimmerman, Center for Global Health and Diseases, Case Western Reserve University School of Medicine, Cleveland, OH, E-mails: vog.cdc@5xgv, ude.esac@08nxk, ude.esac@32zma, ude.esac@41kxj, and ude.esac@zap. John Bosco Keven, Edward Thomsen, and Lisa J. Reimer, Papua New Guinea Institute of Medical Research-Madang, Madang, Papua New Guinea, E-mails: moc.liamg@nevekbj, firstname.lastname@example.org, and email@example.com. Peter Siba, Ivo Mueller, and Manuel W. Hetzel, Papua New Guinea Institute of Medical Research-Goroka, Eastern Highlands, Papua New Guinea, E-mails: firstname.lastname@example.org, mf.liamtsaf@relleumovi, and email@example.com.