RDL mutations predict multiple insecticide resistance in Anopheles sinensis in Guangxi, China

Background Anopheles sinensis is a major vector of malaria in China. The gamma-aminobutyric acid (GABA)-gated chloride channel, encoded by the RDL (Resistant to dieldrin) gene, is the important target for insecticides of widely varied structures. The use of various insecticides in agriculture and vector control has inevitably led to the development of insecticide resistance, which may reduce the control effectiveness. Therefore, it is important to investigate the presence and distribution frequency of the resistance related mutation(s) in An. sinensis RDL to predict resistance to both the withdrawn cyclodienes (e.g. dieldrin) and currently used insecticides, such as fipronil. Methods Two hundred and forty adults of An. sinensis collected from nine locations across Guangxi Zhuang Autonomous Region were used. Two fragments of An. sinensis RDL (AsRDL) gene, covering the putative insecticide resistance related sites, were sequenced respectively. The haplotypes of each individual were reconstructed by the PHASE2.1 software, and confirmed by clone sequencing. The phylogenetic tree was built using maximum-likelihood and Bayesian inference methods. Genealogical relations among different haplotypes were also analysed using Network 5.0. Results The coding region of AsRDL gene was 1674 bp long, encoding a protein of 557 amino acids. AsRDL had 98.0% amino acid identity to that from Anopheles funestus, and shared common structural features of Cys-loop ligand-gated ion channels. Three resistance-related amino acid substitutions (A296S, V327I and T345S) were detected in all the nine populations of An. sinensis in Guangxi, with the 296S mutation being the most abundant (77–100%), followed by 345S (22–47%) and 327I (8–60%). 38 AsRDL haplotypes were identified from 240 individuals at frequencies ranging from 0.2 to 34.8%. Genealogical analysis suggested multiple origins of the 345S mutation in AsRDL. Conclusions The near fixation of the 296S mutation and the occurrence of the 327I and 345S mutations in addition to 296S, in all the nine tested An. sinensis populations in Guangxi, strongly indicate a risk of multiple insecticide resistance. The haplotype diversity plus genetic heterogeneities in the geographical distribution, and multiple origins of AsRDL alleles call for a location-customized strategy for monitoring and management of insecticide resistance.


Background
Guangxi Zhuang Autonomous Region was once a malaria-endemic region with more than 5 million malaria patients per year being recorded before 1949 [1]. After continued efforts, particularly the "National Malaria Control Programme", the "Basically Eliminating Malaria" and the "Action Plan of Malaria Elimination (2010-2020)" [2,3], the malaria burden has been substantially reduced in Guangxi [4]. However, the risk of malaria re-emergence remains, because many imported malaria cases have been reported due to ever increasing cross-border population migration, and the natural environment in Guangxi, consisting of rice fields, provides many mosquito breeding sites [3][4][5][6].
Chemical control of vectors has played an important role in malaria control and elimination [7]. In China, Anopheles sinensis is a major vector of malaria. The use of insecticides in the control of An. sinensis itself, and in agriculture, has inevitably led to increasing insecticide resistance in Chinese An. sinensis [8][9][10]. The development of insecticide resistance can reduce the effectiveness of vector control, therefore, understanding the levels of resistance and mechanisms responsible for this resistance is a core concept to effectively manage An. sinensis. However, the insecticide resistance issue in An. sinensis has so far received little attention. Although the distribution and frequency of insecticide resistance-conferring mutations in the acetylcholinesterase and voltage gated sodium channel were determined in An. sinensis collected extensively across Guangxi in previous studies [11,12], other questions related to insecticide resistance remain to be answered.
The insect gamma-aminobutyric acid (GABA) receptor subunit encoded by the RDL (Resistant to dieldrin) gene plays a central role in neuronal signaling and is involved in various processes [13]. RDL has been the target for insecticides of various chemical structures such as cyclodienes and fipronil [14]. RDL is also a potential secondary target for neonicotinoids and pyrethroids [13]. Previous studies have documented that mutations in RDL are associated with insecticide resistance in multiple insect species. For example, the replacement of a single alanine (Ala301 in Drosophila) at the start of transmembrane segment M2 with either serine or glycine makes the Drosophila strains resistant to dieldrin [15]. Similarly, the A296S substitution (analogous to position 301 in Drosophila melanogaster) was associated with dieldrin resistance in mosquitoes such as Anopheles arabiensis, Anopheles funestus, Anopheles stephensi, Aedes aegypti, Aedes albopictus and Culex pipiens quinquefasciatus [16][17][18][19]. Another replacement of the equivalent residue (A301N) was observed in the fipronil-resistant planthopper Laodelphax striatellus [20,21]. Besides mutation at residue 296, a substitution (V327I) was reported on a background of 296S in dieldrin-resistant An. funestus [18]. In addition, two different substitutions at residue 350 (T350M and T350S) were observed in Drosophila linked with 301G [22] and 301S [23], respectively. Recently, a second RDL mutation (R300Q) in combination with A302S was found to be associated with much higher levels of fipronil resistance (237-fold) in Nilaparvata lugens [24].
Until now, there is no published report about RDL in the malaria vector An. sinensis. In this study, the RDL gene in An. sinensis (AsRDL) was identified, and the naturally occurring genetic mutations in AsRDL and their geographical distribution in Guangxi were investigated, in the hope of predicting their both historical role in cyclodiene (e.g. dieldrin) resistance and potential influence in conferring resistance to currently used insecticides such as fipronil.

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 previously [12].
Anopheles sinensis individual mosquitoes were morphologically identified [25], and put into 100-μl Eppendorf tubes containing 100% ethanol solution at 4 °C. The accuracy of morphological classification was confirmed by molecular detection of ten randomly selected adults from each population using the rDNA-ITS2 method [26].

Sequencing of the full length RDL coding sequence from Anopheles sinensis
Total RNA was extracted from ten adults of a laboratory strain of An. sinensis [27] by using TRIzol (Invitrogen, CA, USA) according to the manufacturer's protocol. First-strand cDNA was synthesized from total RNA (1 µg) using PrimeScript First Strand Synthesis Kit according to the manufacturer's instructions (Takara, Dalian, China). The full length open reading frame sequence of the GABA-gated chloride channel gene was amplified by PCR using the primers (AsRDLfullcd-F:ATGTCGCTAACCATCGAAGTTCCGC; AsRDLfullcd-R: TTACTTATCCTCACCGAGCAGCA) commercially synthesized by Invitrogen (China). The PCR mixture (50 μl) consisted of 2 μl cDNA template, 1 μl PrimeStar ® GXL (Takara), 10 μl 5× Buffer, 4 μl 2.5 μM dNTPs, 1 μl 10 mM each primer, and 31 μl ddH 2 O. The thermal cycling profile consisted of an initial step of denaturation at 95 °C for 3 min; followed by 36 cycles of 98 °C for 10 s, 55 °C for 15 s, 68 °C for 2 min; and a final extension step at 68 °C for 10 min. The PCR products were gel purified (Takara, Dalian, China) and then subcloned into pEasy-T1 and transformed the Escherichia coli Trans 5α strain (Transgen, Beijing, China). The nucleotide sequences of As-RDL gene were identified by direct sequencing of PCR products and/or clone sequencing from two directions.

Haplotype reconstruction and confirmation
Genotype data were carefully inspected by using Chromas 2.01 (Technelysium Pty Ltd, Australia). AsRDL haplotypes were reconstructed by software PHASE2.1 [29]. The presence (and accuracy) of the haplotypes that were reconstructed by PHASE in heterozygotes was confirmed by sequencing of three to five clones.

Genealogical analysis
The AsRDL haplotypes were used to build the phylogenetic topologies by maximum-likelihood (ML) and Bayesian inference (BI) methods. The best-fitting substitution model was determined under the Bayesian information criteria (BIC) using jModelTest 2.1.1 [30]. The best-fitting model of nucleotide substitution for AsRDL fragment was TVM+I. PhyML 3.0 was used to perform ML analyses [31]. The non-parametric bootstrap with 1000 replicates was used to estimate the support for each internal branch of the ML phylogeny. The program MrBayes 3.2.1 was used for Bayesian phylogenetic inference [32]. Markov chains were run for 10 7 generations with two independent runs, and trees were sampled every 1000 generations. The initial 25% of trees were discarded as burn-in, and the remaining trees (15,002) were used to estimate the majority-rule consensus tree and the Bayesian posterior probabilities (BPP). Genealogical relations among different haplotypes were also analysed using Network 5.0 [33].

Identification of AsRDL gene
The coding region of AsRDL gene was 1674 bp long, encoding a protein of 557 amino acids. Four non-synonymous variations that result in amino acid substitutions (R119G, I162V, I176V, I278V) were identified in cDNA sequences (Fig. 1). Mapping the coding region to the genomic sequence (KE525297.1) showed that the AsRDL gene consists of nine exons and eight introns (Fig. 2).

Identification of naturally occurring genetic mutations in AsRDL
A 184 bp-in-length sequence of partial exon 7 (namely 7F), and a 494 bp region containing 81 bp in intron 7 and 413 bp in exon 8 (namely 8R), were obtained from 240 individuals, respectively. Three polymorphic sites were observed in 7F (Fig. 4). From the 8R sequences, ten and eighteen polymorphic sites were identified in intron 7 and exon 8, respectively (Fig. 4).
From these sequencing data, 12 non-synonymous polymorphic sites were recorded. Two non-synonymous mutations within exon 7 were identified in the first nucleotide (T/G) of codon 296 and in the first nucleotide (A/G) of codon 327, which lead to A296S and V327I substitutions, respectively (Fig. 5). Ten non-synonymous mutations were observed within exon 8, which result in 345T/S, 349L/M, 356R/G, 357K/L/M/N, 360F/L, 408G/H and 473T/A substitutions, respectively (Fig. 1).
The distribution frequencies of these substitutions were presented in Table 1. Among these substitutions, 356G, 357L, 357N and 408H were detected in only one population at a very low frequency. In contrast, three putative resistance-related substitutions (A296S, V327I and T345S) were widely distributed with relatively higher frequencies (Table 1 and Fig. 6).
Based on the genotypes of the three resistance-related sites, eleven combinations were observed (Table 2). Most individuals carried at least one resistant mutation (combination 1-10). A few individuals from Baise (6.8%) and Liuzhou (12.5%) were susceptible homozygotes (combination 11). Two combinations (2 and 3) contained mutations in all the three positions. Notably, either 327I or 345S was present only in individuals carrying 296S, and combinations 3, 6 and 7 were distributed in all the nine populations.
The distribution frequencies of the AsRDL genotypes were showed in Fig. 6. Overall, the distribution showed obvious geographical heterogeneity. At position 296, the heterozygotes were present in Baise, Liuzhou and Guilin, while the susceptible homozygotes were found only in Baise and Liuzhou. Notably, the resistant homzygotes (296SS) was fixed in six of the nine populations. In contrast, the frequencies of resistant homozygotes (327II or 345SS) were much lower. Neither individual homozygous for 327II was detected in Hezhou and Wuzhou, nor mosquito homozygous for 345SS was observed in Guigang (Fig. 6).

Haplotype diversity of AsRDL gene
The sequence data of the 8R fragment that covers both partial intron 7 and exon 8 with a length of 494 bp were used in haplotype analysis in order to obtain more informative result than data from the 7F fragment. From the sequence data of 240 individuals, 38 AsRDL haplotypes Exon/intron organization of the AsRDL gene. E represents exon and dash line represents intron. The numbers above the box or below the dash line represent the length (base pair) of different parts of AsRDL gene. The intron phase is indicated using braced numbers: 0 means intron splicing between codons, 1 means intron splicing between the first and second nucleotides of a codon, and 2 means intron splicing between the second and third nucleotides of a codon were identified ( Table 3). The 38 haplotypes varied in frequency from 0.2 to 34.8%. Overall five common haplotypes (H1, H2, H3, H17 and H25) accounted for 67% of the haplotype diversity, while sixteen haplotypes were observed only once. The two most common haplotypes (H3 and H17) occurred in all the nine populations.

Genealogical analysis of AsRDL haplotypes
The network analysis revealed a complex reticulate genealogical relation of the 38 AsRDL haplotypes, and suggested independent mutation events leading to 345S-containing haplotypes (Fig. 7). For example, the resistant H17 was possibly derived from H3, while H32 from H25. Although the phylogenetic tree had lower resolution, the 38 haplotypes were roughly grouped into two major clades (63% bootstrap value) separated by 3 nucleotide substitutions in intron 7 (Fig. 8).

Discussion
In insects, GABA receptors are ligand-gated chloride channels consisting of presumably the homopentameric RDL subunit encoded by the RDL gene. RDLs mediate synapse inhibition in the central nervous system and are important targets for insecticides of widely varied structures such as cyclodienes and fipronil [14]. In this study, Fig. 6 The sampling location (a), and distribution frequencies of individual genotypes of AsRDL in Guangxi (b-d). BS Baise, HC Hechi, HZ Hezhou, GG Guigang, GL Guilin, LZ Liuzhou, NN Nanning, WZ Wuzhou, YL Yulin. n represents the sample size the RDL gene was identified in An. sinensis. The putative AsRDL protein has 557 amino acid residues and shares the common structural features with known Cys-loop ligand-gated chloride channels (Fig. 1). Similar to the finding obtained in An. funestus [18], only one GABAreceptor isoform is found in An. sinensis in this study.
The RDL subunit has four transmembrane segments (M1, M2, M3 and M4). Point mutations in the M2 and M3 have been documented to confer low sensitivity of the binding site and thus insecticide resistance. In this study, the insecticide resistance associated mutation 296S, and two other mutations (327I and 345S) that have been documented to be linked to 296S [18,23], were identified and widely distributed in An. sinensis populations (Table 1 and Fig. 5). The prediction of the three mutations to be resistance-related is made by reference to other species, thus further functional study is required to confirm it. Whether other mutations (356G, 357L, 357N and 408H) in An. sinensis, which are narrowly distributed at a low frequency, contribute to insecticide resistance is unclear.
Point mutations (Ala to Ser/Gly/Asn) in RDL homologues at the site equivalent to 301 in Drosophila have been documented to provide resistance to cyclodienes and phenylpyrazoles in several species [17][18][19][20][21][22][23][24]34]. From the 240 individuals from nine geographical locations, the 296S mutation was detected at a very high frequency. Notably, 100% of individuals in six of the nine populations carried this mutation in the homozygous form. These data indicate a risk of insecticide to cyclodienes and phenylpyrazoles in Guangxi An. sinensis populations. However, the GCT to GGT mutation (inducing the A296G change) found in An. gambiae [17], and the replacement of Ala with Asn observed in plant hoppers [20,21], were not detected in An. sinensis samples. Given that cyclodienes have now been withdrawn from control programmes, the high frequency of the dieldrin resistance-associated A296S mutation in An. sinensis in Guangxi may be the result of the increasing use of insecticides targeting the GABA receptor (e.g. ethiprole) in the agriculture [35].
The 327I mutation which is located in the intracellular loop between TM2 and TM3 has been reported to coexist with 296S in dieldrin-resistant An. funestus [18]. Interestingly, a point mutation identical to An. funestus (GTA to ATA) was identified in An. sinensis in this study (Fig. 5). The impact of the 327I mutation is unknown and worthy of further investigation.
Two mutations at the site 350 in the third membrane spanning domain of RDL have been documented in Drosophila [22,23]. T350 M was observed in fipronil resistant Drosophila simulans with Gly301 [22], while T350S was detected in a Drosophila melanogaster line (Ral-491) that is homozygous for Ser301 [23]. Interestingly, conserved mutations (A296G + T345 M) were also found in An. gambiae [13]. In this study, parallel sinensis. This observation lends support to the notion that there is a functional role for mutations at the 345 site [13,23]. Sequence analyses revealed that AsRDL gene had diverse polymorphisms. From the fragment sequences of AsRDL from 240 mosquitoes, 38 haplotypes were identified ( Table 3). The number and frequency of AsRDL haplotypes were different among the nine populations from Guangxi. The most abundant haplotype was H3 (34.8%), followed by H17 (13%). Both H3 and H17 were distributed in all the nine populations tested, and shared the type-2 intron. Similar to the status of AS-VGSC [12], the high polymorphism of AsRDL may be partly explained by their large population size, wide distribution range of An. sinensis [36], diverse natural landscapes and different local insecticide selective pressure in Guangxi.
The genealogical relation and phylogenetic tree analyses support two independent origins of the resistance related 345S mutation in AsRDL (Figs. 7, 8). The two most common susceptible haplotypes (H3 and H25) could be considered as possible ancestors of resistant haplotypes respectively. The six resistant haplotypes (H17-22) were derived from their common susceptible ancestor H3 with one to four mutational steps, whereas the other two resistant haplotypes (H32 and H33) were the result of one and two mutational steps respectively from H25.