Vgsc-interacting proteins are genetically associated with pyrethroid resistance in Aedes aegypti

Abstract

Association mapping of factors that condition pyrethroid resistance in Aedes aegypti has consistently identified genes in multiple functional groups. Toward better understanding of the mechanisms involved, we examined high throughput sequencing data (HTS) from two Aedes aegypti aegypti collections from Merida, Yucatan, Mexico treated with either permethrin or deltamethrin. Exome capture enrichment for coding regions and the AaegL5 annotation were used to identify genes statistically associated with resistance. The frequencies of single nucleotide polymorphisms (SNPs) were compared between resistant and susceptible mosquito pools using a contingency χ2 analysis. The -log10(χ2 p value) was calculated at each SNP site, with a weighted average determined from all sites in each gene. Genes with -log10(χ2 p value) ≥ 4.0 and present among all 3 treatment groups were subjected to gene set enrichment analysis (GSEA). We found that several functional groups were enriched compared to all coding genes. These categories were transport, signal transduction and metabolism, in order from highest to lowest statistical significance. Strikingly, 21 genes with demonstrated association to synaptic function were identified. In the high association group (n = 1,053 genes), several genes were identified that also genetically or physically interact with the voltage-gated sodium channel (VGSC). These genes were eg., CHARLATAN (CHL), a transcriptional regulator, several ankyrin-domain proteins, PUMILIO (PUM), a translational repressor, and NEDD4 (E3 ubiquitin-protein ligase). There were 13 genes that ranked among the top 10%: these included VGSC; CINGULIN, a predicted neuronal gap junction protein, and the aedine ortholog of NERVY (NVY), a transcriptional regulator. Silencing of CHL and NVY followed by standard permethrin bottle bioassays validated their association with permethrin resistance. Importantly, VGSC levels were also reduced about 50% in chl- or nvy-dsRNA treated mosquitoes. These results are consistent with the contribution of a variety of neuronal pathways to pyrethroid resistance in Ae. aegypti.

About the Authors:

Corey L. Campbell

Roles Formal analysis, Investigation, Methodology, Validation, Writing – original draft

* E-mail: [email protected]

Affiliation: Department of Microbiology, Immunology and Pathology, Colorado State University, Fort Collins, Colorado, United States of America

ORCID logo http://orcid.org/0000-0002-6384-1259

Karla Saavedra-Rodriguez

Roles Formal analysis, Investigation, Methodology

Affiliation: Department of Microbiology, Immunology and Pathology, Colorado State University, Fort Collins, Colorado, United States of America

Tristan D. Kubik

Roles Formal analysis

Affiliation: Department of Microbiology, Immunology and Pathology, Colorado State University, Fort Collins, Colorado, United States of America

Audrey Lenhart

Roles Resources, Writing – review & editing

Affiliation: Division of Parasitic Diseases and Malaria, Center for Global Health, Centers for Disease Control and Prevention, Atlanta, Georgia, United States of America

ORCID logo http://orcid.org/0000-0002-9156-4228

Saul Lozano-Fuentes

Roles Formal analysis

Affiliation: Department of Microbiology, Immunology and Pathology, Colorado State University, Fort Collins, Colorado, United States of America

ORCID logo http://orcid.org/0000-0003-1517-6853

William C. Black IV

Roles Formal analysis, Funding acquisition, Resources, Writing – review & editing

Affiliation: Department of Microbiology, Immunology and Pathology, Colorado State University, Fort Collins, Colorado, United States of America

Introduction

The major arbovirus vector Aedes aegypti continues to pose significant threats to human health in tropical and subtropical urban areas [1,2]. In southern Mexico, where dengue is hyperendemic [3], insecticides used for Ae. aegypti control include organophosphates for larval control and pyrethroids for adults [4]. The voltage-gated sodium channel (Vgsc, aaeNAV, LOC5567355) is one target of pyrethroid toxicity [5]. Increased resistance in Ae. aegypti to pyrethroid insecticides has been documented in many locations [6–9], and is associated with the presence of specific VGSC alleles. In addition, multiple lines of evidence indicate that metabolic resistance mechanisms may be equally important as target-site resistance [10–13]. Metabolic resistance to pyrethroids is generally described as a number of Vgsc-independent mechanisms derived from gene duplication, transcript overexpression, sequestration or increased reduction/oxidation (redox) activity. Cytochromes P450 (CYP), episilon class glutathione S-transferases (eGST) and esterases (EST) are major contributors to metabolic resistance, which are stimulated in response to oxidative stress [14–16].

To further increase our understanding of pyrethroid resistance mechanisms, the aim of the present study was to use association mapping to characterize allele frequency differences in genes for two Ae. aegypti aegypti field locations. One was located in the community of Viva Caucel in the city of Merida in the Mexican state of Yucatán and the other in Vergel in the same city. These locations were chosen because the association between VGSC genotype and permethrin knockdown resistance (kdr) phenotype were known for each site [9]. Specifically, Vgsc target site resistance loci V1016I and F1534C (house fly protein annotation), were present in each population. According to Vera-Maloof et al. 2015, in the Viva Caucel collection, the V1016I and F1534C mutations were present in 75% and 93% of individuals assayed, respectively. Similarly, in the Vergel collection, V1016I and F1534C were present in 80% and 98% of individuals assayed, respectively. Thus, in both locations, the F1534C allele was close to fixation. These characteristics indicated that the natural populations from which the collections were derived had already been selected for permethrin resistance.

Our hypothesis was that pyrethroid resistance-associated genes could be those with 1) chromosomal physical proximity to VGSC, resulting from a selective sweep, 2) selection due to direct activity in resistance or 3) genetic drift. For the present study, exome capture genomic DNA (gDNA) high throughput sequencing (HTS) data from both locations with two pyrethroid insecticides were analyzed to identify SNPs that differed significantly in frequency between susceptible and resistant pools in a given collection. We then highlighted those genes with significant genetic association to pyrethroid resistance. Genetic association refers to the concurrent presence of genotypic polymorphisms, as indicated by–log10(p value) scores from χ2 contingency table calculations, and the pyrethroid resistance phenotype. We used orthology information to assign functional attributes for all coding genes and used that information to identify those with statistical association to pyrethroid resistance. In addition, the gene set common to all treatment groups was subjected to GSEA [17]. Identification of genes common to both permethrin- and deltamethrin-resistant mosquitoes would give insight into overall pyrethroid resistance mechanisms. The identification of common mechanisms would also build support for increasing variability of rotating pyrethroid spraying regimens to reduce selection of resistant loci [18].

Results and discussion

Treatment and high throughput sequencing

A total of 3 pyrethroid treatment groups were evaluated for this work. Two groups, from the Viva Caucel location, were treated with either permethrin or deltamethrin. The third group, from the Vergel location, was treated with permethrin. Our rationale for evaluating combined data for type 1 and type 2 pyrethroids was 1) to identify genomic loci that are common among both treatment types and 2) reduce the identification of population-specific loci. Previously published bottle assay results are summarized in Table 1 [8,9]. Dual independent replicate pools of adult females (n = 25 per pool) from each population were prepared for a total of 12 library sets, comprising four from each location/treatment group (see Fig 1 legend). Then they were exome-capture enriched for genic regions [19,20]. From 27.25–86 million trimmed paired end reads were obtained for each replicate library (Table 2, S1 Table), and 94 to 98% of all trimmed reads aligned to the AaegL5 chromosome length reference sequence [21], which allowed for mapping of genes to their approximate physical location and eliminated the problems with unmapped genes that we experienced in previous studies [19].

[Figure omitted. See PDF.]

Fig 1. Genes with statistical association to pyrethroid resistance.

Genes (FDR, p<0.01, -log10 (p value) of χ2 contingency calculation ≥ 4.0) in each of three collections, Viva Caucel permethrin (VCP, gray circle), Viva Caucel deltamethrin (VCD, open circle), and Vergel permethrin (VeP, black circle) were plotted according to the average physical position of each gene (AaegL5 genome build). X-axis shows beginning and end locations of polymorphisms on each chromosome.

https://doi.org/10.1371/journal.pone.0211497.g001

[Figure omitted. See PDF.]

Table 1. Bottle assays.

https://doi.org/10.1371/journal.pone.0211497.t001

[Figure omitted. See PDF.]

Table 2. Polymorphisms and coverage.

https://doi.org/10.1371/journal.pone.0211497.t002

At each SNP, χ2 contingency table calculations determined statistical association of a given SNP to pyrethroid resistance. Resistance associated gene-wise scores were calculated from the weighted average–log10(p value) of the top 5% of SNP scores per gene (see Methods). Prior to further analysis, a multiple-testing adjustment (false discovery rate, FDR, alpha = 0.01) [22] and -log10(χ2 p value) ≥ 4.0 cut-offs were applied. Collection-specific details are shown in S2 Table. Gene-wise association scores above the cut-off were mapped relative to each gene’s physical position for each collection (Fig 1).

Groups with high association to pyrethroid resistance

To identify high association genes and functional categories associated with pyrethroid resistance among all three treatment groups, the intersection of genes above the thresholds were determined (Fig 2A, n = 1,054). One of these was a non-coding RNA and therefore was excluded from further analysis, leaving 1053 coding genes (S2 Table). This high association set was annotated by predicted function based on orthology to model organisms [23]. Overall functional groups for these genes are displayed in Fig 2B. Within the common gene set, the relative physical position of each gene and the corresponding genetic association values by collection are shown in Fig 3. Strikingly, though the common gene subset spanned all three chromosomes (S2 Table), most genes (1,023 of 1,053) were on chromosomes 2 and 3.

[Figure omitted. See PDF.]

Fig 2. Genes with statistical association to resistance were common among all three collections (n = 1054).

A) Venn diagram indicates common genes among all three treatment groups; Viva Caucel permethrin (VCP), Viva Caucel deltamethrin (VCD), and Vergel permethrin (VeP). B) Functional categories for genes in the high association set are listed from least abundant (top) to most abundant in the order of the diagram slices. X-axis shows beginning and end locations of polymorphisms on each chromosome.

https://doi.org/10.1371/journal.pone.0211497.g002

[Figure omitted. See PDF.]

Fig 3. Physical positions of genes in high association subset (n = 1054).

Genes were plotted according to their physical positions (AaegL5 reference). Y axis indicates weighted gene-wise average–log10(p value) of χ2 contingency table calculations. Red circle indicates approximate chromosomal location of the permethrin target site VGSC (nt positions 315,926,360–316,405,639) in VCP, VCD and VeP; red arrows show approximate location of NVY (nt 301,564,432–301,941,839) and CHL (nt 312,986,638–313,140,457).

https://doi.org/10.1371/journal.pone.0211497.g003

VGSC and associated genes.

As expected, our tests revealed high statistical association of VGSC, the target of permethrin, in all three treatment groups (S2 Table). VGSC is located on the q arm of chromosome 3 [21,24] (Fig 3). Intriguingly, some genes in the high association group interact with VGSC in other organisms. For example, CHL (LOC5578279), a neuronal-specific zinc-Cys2His2 type transcriptional regulator [25], genetically interacts with the drosophilid ortholog of VGSC (PARA) [26]. Synaptic control of motoneuron function has been characterized in dipteran models of epilepsy, wherein increased frequency and regularity of calcium signaling waves, synaptic excitation, are indicative of pathology [26–29]. For example, in a drosophilid epilepsy model, reduced expression of CHL further increased sensitivity of the PARA bang-sensitive (PARA bss) allele, which has a high basal sensitivity to paralysis. PARA bss was also sensitive to alterations in calcium levels [29], which is consistent with a functional interaction between Vgsc and calcium channels. Similarly, deltamethrin has been shown to reduce calcium channel efflux in mammalian and mosquito cell culture [30–32]. Consistent with this evidence, a calcium channel alpha subunit (caα1, LOC5564339) was also within the high association group.

The ankyrin protein domain is present in a large variety of structural proteins that link transmembrane proteins to the actin cytoskeleton [33,34]; they serve a variety of functions, including cytoskeletal anchoring and mechano-transduction of sensory responses [35,36]. A subset of these proteins physically interacts with Vgsc in mammals [28]. Five ankyrin domain-containing proteins were also identified in our high association gene set (S2 Table). Finally, NEDD4, also a high association gene, has predicted physical interactions with NAV (Vgsc) in mammals (S2 Table) [37]. Additional support for NAV-Nedd4 interactions was also observed in its regulation of cardiac NAV [38]. NEDD4 codes for E3 ubiquitin-protein ligase, which regulates degradation of plasma membrane proteins with PY motifs [39]. Aedine Vgsc has a partial PY motif (PPS, amino acids 1170–1172) rather than the canonical PY of mammalian NAV proteins (PPSYDSV) [40]. Until now, no support for VGSC-NEDD4 interactions has been reported for dipterans.

The importance of each of the associations summarized above was further substantiated by GSEA of each corresponding functional category. We tested each functional group for relative enrichment within the high association set compared to all coding genes. This was done using a hypergeometric test, which models a binomial distribution and tests the probability of a gene being represented more often than expected by chance [41]. Specific results are described below in the context of each functional group.

Selected gene silencing and phenotype validation.

To further explore the role of selected genes in pyrethroid resistance, we evaluated 2 putative transcription regulators in the high association data set. We hypothesized that when silenced, relevant components would result in increased pyrethroid susceptibility. One gene, NVY (LOC5563881), codes for a predicted central and peripheral nervous system-specific transcriptional regulator [42] with dual functions in axon guidance [43] and transcriptional repression [44]. The other gene, CHL, described above, is genetically associated with VGSC in drosophilids. We used RNA interference-induced gene silencing via dsRNA-injection to temporarily silence each gene in Vergel colony mosquitoes (F24). Silencing was confirmed relative to levels in non-specific dsRNA (β-Gal-dsRNA) treated controls (Fig 4A), and mosquitoes were subjected to a 1.5 μg permethrin discriminating dose in standard CDC bottle assays [45]. Both nvy-dsRNA and chl-dsRNA injected groups were significantly more susceptible to permethrin knockdown than (β-GAL- dsRNA)-injected controls (analysis of variance (ANOVA), p<0.0001) (Fig 4B). In contrast, PBS-injected controls showed sensitivity and viability levels that paralleled that of the β-GAL-injected group, indicating that exogenous dsRNA treatment alone had no effect (S1 Fig). The proportion of mosquitoes with lethal phenotypes at 24 hours post-treatment was also significantly greater in nvy-dsRNA and chl-dsRNA groups than that of β-GAL-dsRNA-injected mosquitoes (Fisher's exact test, p = 0.001 and p<0.0001, respectively) (Fig 4B).

[Figure omitted. See PDF.]

Fig 4. CHL- and NVY-silenced mosquitoes show increased susceptibility to permethrin and reduced VGSC transcript levels.

A) Gene silencing validation. Mosquitoes were injected with either chl-dsRNA or nvy-dsRNA and processed in pools (n = 5) to isolate RNA at 3 days-post-treatment (dpt). Target transcripts were silenced an average of 93% and 78% in chl-dsRNA and nvy-dsRNA mosquitoes, respectively, relative to βGal-dsRNA-injected controls. B) Mosquitoes were injected with βGal-dsRNA (n = 113), chl-dsRNA (n = 172) or nvy-dsRNA (n = 129). At 3 dpt, each replicate was subjected to ~1.5 ug permethrin in a CDC bottle assay; knockdown was recorded at 10 minute intervals. There was a significant difference between chl and nvy compared to β-Gal controls ***, (ANOVA, CHL, p<0.0001; NVY, p<0.0001). At 24 hours, mortality was also significantly elevated (Fisher’s Exact Test, CHL, p<0.0001, NVY, p = 0.001). Left Y axis indicates percent knockdown from 10–60 minutes. Right Y-axis indicates percent mortality at 24 hours. Error bars indicate standard error of the mean (SEM). Bottle assay data represent a compilation of 4 to 5 biological replicates. C) Pools of 5 chl-dsRNA and nvy-dsRNA treated mosquitoes were subjected qRT-PCR of VGSC relative to βGal-dsRNA injected controls. Results shown are the average of 3 biological and 3 technical replicates.

https://doi.org/10.1371/journal.pone.0211497.g004

Analysis of VGSC transcript levels showed a ~50% reduction following nvy-dsRNA or chl-dsRNA treatment compared to β-GAL-dsRNA-injected controls, indicating that each transcriptional regulator plays a role in modulating VGSC transcript levels (Fig 4C). The association mapping data was interrogated to identify CHL and NVY non-synonymous mutations in each of the collections. No non-synonymous mutations (LOD >4.0) were identified in either CHL or NVY in the VCD collection, suggestive of cis-regulatory effects. A single common mutation at A2T occurred in VCP and VeP, with a few other isolated amino acid substitutions in CHL (Fig 4D).

Detoxification and metabolic genes with high association to pyrethroid resistance.

Detoxification is a common pyrethroid resistance mechanism [14,15,46–48], and several key enzymes in this category are within the metabolic functional group. Overall, seventy-eight high association genes had predicted activity in metabolic processes, as shown in Fig 2B (S2 Table), (78/1053), hypergeometric test p value = 0.02). Though type I, or non-cyano group bearing pyrethroids, eg., permethrin, and type II, or alpha-cyano group bearing, eg., deltamethrin, have distinct toxicity mechanisms [30], individual detoxification genes may respond to both types I and II [49]. The major pyrethroid detoxification proteins are mixed function oxidases and alpha- and beta-esterases (S2 Table) [50,51].

In anophelines, an acetylcholine esterase resistant allele (ACE-1R) can partially reduce deltamethrin sensitivity, and ace-1R works synergistically with VGSC resistant alleles to substantially reduce sensitivity to permethrin [52]. However, there is little data to indicate an association between selection of acetylcholine esterase alleles and permethrin resistance in Ae. aegypti. Nevertheless, two acetylcholine esterases were included in the high association group (S2 Table, LOC5570776 and LOC5578456). Three genes coding for carboxylesterases (COEs, LOC5564561, LOC5576941, LOC5571034) were also in this group. COEs have been tied to organophosphate resistance mechanisms and, to a lesser extent, to deltamethrin resistance [53]. In addition to supporting a role for these genes in pyrethroid detoxification, their presence in the common gene set is consistent with the idea that there may be some overlap between type I and II pyrethroid detoxification mechanisms. In contrast, there were no epilson class glutathione S transferases in the high association group, which is surprising, given that they have been shown to be pertinent to resistant phenotypes [12,14,15,54–56].

Over-representation of the metabolic gene category is consistent with studies in drosophilids wherein amino acid catabolism was tied to permethrin sensitivity [57]. Specifically, Brinzer et al. found a decrease in proline levels in resistant larvae that had been exposed to permethrin, which could be explained by either increased conversion from proline to glutamate or inhibition of the proline precursor ornithine. Conversion of proline to glutamate has also been reported upon exposure to malathion or DDT [58,59]. In the current work, the aedine ortholog of 5-oxoprolinase, responsible for the degradation of glutamate, which would be required if it was in excess, is present in the high association group (S2 Table) [57].

Additional metabolic pathway genes in the high association group included lipid biosynthetic genes (S2 Table), for which transcripts are reportedly enriched in temephos-resistant mosquitoes [60]. Interestingly, there were also multiple genes that participate in lipid metabolic processes, ie., 3 lipases and 4 phospholipases (S2 Table). Increased lipase activity has been tied to insecticide resistance in Tribolium spp [61]. One possible explanation of this is that the catabolism of choline phospholipids is protective against the loss of choline, caused by the release of acetylcholine from synaptosomes [62] and inhibition of acetyl cholinesterase [63] that occurs in pyrethroid-susceptible insects. Similarly, a sphingomyelin phosphodiesterase (SMAse) activator domain protein was also in the high association group. SMAses are hydrolases that break down sphingomyelin into phosphocholine and ceramide (reviewed in [64]).

Some pyrethroid detoxification proteins are in the redox functional category. For our analyses, the redox and mitochondrial-associated gene subsets were combined, as key redox genes are localized to mitochondria. This category was not over-represented compared to all coding genes. Nevertheless, genes within this high association group included eighteen confirmed or putative CYPs, eight of which are CYP4c orthologs. CYP4c is an insect-specific member of the CYP4 family; CYP4 family proteins are best known for oxidation of fatty acids and localize to peroxisomes and the endoplasmic reticulum. Recent work has also reported their involvement in permethrin resistance [16]. Though CYPs do not directly detoxify permethrin, they may function downstream of esterases, which are responsible for much of permethrin catabolism (reviewed in [50]). Moreover, CYP6 and CYP9 (S2 Table) are also well documented for their involvement in resistance to both deltamethrin and permethrin [48,49].

Though CYPs are the best studied of the redox functional group [10,11,13,15], mitochondrial function genes also have been implicated in permethrin resistant anopheline mosquitoes, as indicated by the presence of elevated reactive oxygen species [46]. In addition, key electron transport gene cytochrome B (CYTB) transcripts are elevated in Ae.aegypti following permethrin treatment [65], consistent with a requirement for increased mitochondrial function. Here, there were five genes coding for cytochrome B and C associated genes in the high association subset (S2 Table).

Transport genes with high association to resistance.

The transport functional category was the most statistically significant of all over-represented categories evaluated (120/1053, hypergeometric test p value = 1.3 x 10−6). This is intriguing, because one would not necessarily expect the category to be over-represented among all coding genes. Because the major pyrethroid target site, Vgsc, is a sodium transporter [5,66–71], it is possible that off-target receptor/channel functions are perturbed upon pyrethroid exposure and thus selected. Indeed, 55 of 120 genes in this category code for transporter, receptor or channel proteins. An alternate possibility is that selection of transport genes represents compensatory mutations to restore overall membrane homeostasis disrupted by pyrethroid binding to Vgsc in polarized neuronal cells (reviewed in [72]).

Repair/replication/transcription/translation genes.

The largest specific functional category was for genes that affect DNA replication/ repair/transcription/translation (RRTT) functions. About 12.7% of all common genes were in this category, and though they were not over-represented relative to all coding genes (134/1053), there was a trend toward enrichment (hypergeometric test p value = 0.066). Of these genes, twenty are predicted to regulate transcription, including CHL, described above. Two orthologs of the transcription factor grauzone were also in this group. Grauzone is a Cys2His2 zinc-finger positive regulator of transcription and is required for meiosis in oogenesis [73]. Its possible role in insecticide resistance is unknown.

The effects of changes in transcriptional regulation of pyrethroid resistance are understudied in mosquitoes. Within the RRTT functional group, ribosomal protein subunits, histone proteins and histone methyltransferases were identified, hinting that insecticide resistance may require alterations from chromatin modifications to the initiation of gene expression to translational processes for establishment of the resistant phenotype.

Cytoskeletal/structural genes.

Cytoskeletal components are crucial to proper neuronal function (reviewed in [74]). Genes coding for predicted cytoskeletal/structural function proteins were not over-represented among the subset common to all three collections. However, specific cytoskeletal/structural category genes with ankyrin domains are predicted to interact with Vgsc [27,28,33] and are genetically associated with PARAbss in Drosophila spp. [27]. In addition, sickie (LOC5564933), a cytoskeletal positive regulator of neuronal axon growth, is also in the high association set [75], as was β-spectrin (β-SPEC) (S2 Table), a cytoskeletal protein important to maintaining neuronal structure in mammals and dipterans [74,76]. In mammals, β-SPEC is a neuronal structural component that physically interacts with NAV channels [33] and also affects presynaptic stability [33].

Signaling.

Genes with predicted function in cellular signaling were significantly over-represented (31/1053, hypergeometric probability p value = 0.013). G-protein coupled receptors (GPCR) are suspected of activating expression of CYPs in permethrin resistant Culex quinquefasciatus [77]. Twelve of 24 genes in this subset were among the high association set; these included GPCRs, GTP-binding proteins and GTPases. Second messenger signaling components are also important to overall control of cytoskeletal changes required for neural synapse function (reviewed in [74]).

Genes with exceptionally high association levels.

To identify the highest resistance-associated genes among all treatment groups (VCP, VCD and VeP), genes among the top 10% of–log10(p value) values in the high association set were identified. Values within this extreme group were all greater than 9.02, 7.95 and 9.85 for VCP, VCD and VeP, respectively. Thirteen genes, including VGSC, were within this subset, and eight of these were located on the p arm of chromosome 3, in the same region as VGSC. Therefore, extreme association on the 3p arm may be due to a selective sweep of VGSC, resulting in incidental selection of genes that have no functional association. However, if this was the case, we would expect a higher percentage of the 737 high association genes located on the chromosome 3p arm to be included in the extreme group. Specifically, the extreme association group represents just 1.7% of the VGSC -proximal high association cluster. This is consistent with the hypothesis that key genes within this subset are specifically selected due to functional attributes rather than by proximity to VGSC alone. Importantly, these observations do not negate the possibility of a selective sweep in this region.

In addition to VGSC, the extreme association group included two transcriptional regulators (NERVY (NVY) LOC5563881 and LOC5579659), a Gr14 gustatory receptor (LOC5575007), a putative lipid binding protein, neural LAZARILLO (LOC5572156), CINGULIN (LOC5577979), a structural component of gap junctions at neuronal synapses [78] and three hypothetical genes. LAZARILLO, a lipid-binding apolipoprotein, though not studied in mosquitoes, is a biomarker for deltamethrin liver toxicity in mammals [79]. For all of these gene products, the mechanisms of association with pyrethroid resistance in mosquitoes remain to be explored.

Genes associated with synaptic function.

Permethrin binding to the neuronal-specific transmembrane protein, Vgsc, could result in synaptic injury. The resulting toxic effects may select for mutations in genes that control pre- and post- synaptic processes. To explore this hypothesis, the high association group was interrogated for the presence of genes with known association to pre- and post-synaptic function. Twenty-one genes with demonstrated associations to this process were identified (S2 Table) [74,80,81]. Calcium signaling is requisite to motor neuron function [82,83], and the characteristic phenotype of pyrethroid treatment is the loss of motor control, due to insecticide binding to Vgsc. Therefore, the genetic association of both CAa1 and VGSC to pyrethroid resistance is consistent with complementary contributions of each protein toward motor synaptic function. In the drosophilid model, PARA expression is stimulated upon limiting synaptic activity levels [83]. Consistent with these observations, the translational repressor, PUM, helps regulate motoneuronal activity, a function that is conserved in mammals and dipterans [84,85]. It does this via control of VGSC (PARA) translation, thus regulating synaptic excitation [83]. Of note, PUM was also present in our high association set, consistent with a role in the selection of resistance (S2 Table). Synaptotagmins were another gene cluster associated with synaptic function; these genes code for synaptic vesicle calcium sensors (S2 Table) [86] that interact with NAV in mammals [87].

Neuronal cytoskeletal genes are crucial to maintenance of synaptic control (reviewed in [74]). Several cytoskeletal and associated regulatory genes important to synaptic activity were in the high association set (S2 Table), including an atypical protein kinase C (aPKC) ortholog (LOC5576493). aPKC stabilizes the microtubule cytoskeleton for proper synaptic function (reviewed in [74]). β-SPEC is an important structural feature of synaptic connections [35] and interacts with NAV channels in mammals [88]. Moreover, neuroligin-1 was also present in the high association group, also consistent with selection at excitatory synapses [89].

Conclusions

We identified gene subsets with significant enrichment in populations resistant to pyrethroids. Importantly, genes were also identified that genetically, functionally or physically interact with Vgsc in other model organisms or in cell culture. Genes which were present in the high association set and for which multiple lines of evidence support genetic or physical interactions with VGSC are highlighted in Fig 5 [25,26,28,32,33,37]. We also found that aedine transcriptional regulators CHL and NVY are associated with permethrin resistant phenotypes. This is consistent with the hypothesis that transcriptional processes are under selective pressure. Both CHL and NVY are located on the p arm of chromosome 3 (red arrows, Fig 3) near VGSC within the area of a proposed selective sweep. The functional associations described in this report are consistent with the idea that specific genes within the VGSC-proximal cluster are associated with resistant phenotypes, are under selective pressure, and are not selected merely due to chromosomal proximity.

[Figure omitted. See PDF.]

Fig 5. Putative interacting genes were inferred from 1) their presence in the high association set, and 2) functional or genetic evidence in model organisms or 3) reported here for NVY and CHL in Ae.

aegypti. ANK, ankyrin-2 [33,35]; aPKC, atypical Protein kinase C [74]; β-SPEC, beta-spectrin [35,74,76]; CAa1, calcium channel alpha subunit [30,31]; CHL, CHARLATAN [25,26]; NEDD4, E3 ubiquitin-protein ligase [37,38]; NVY, NERVY [42–44], PUM, pumilio [27], SYT1, SYNAPTOTAGMIN-1 [87,90]; VGSC, voltage-gated sodium channel [74]. Dotted line indicates indirect interaction due to shared roles in synaptic function. Blue line indicates genetic association in drosophilids. Red line indicates genetic evidence in Ae.aegypti or mosquito cell culture. Brown line indicates evidence in mammals.

https://doi.org/10.1371/journal.pone.0211497.g005

One explanation for over-representation of the “Transport” functional group is that polarized neuronal cells are more sensitive to alteration of transport functions and therefore under selective pressure in pyrethroid treatment areas. Mosquitoes may develop compensatory mutations in genes coding for transport proteins in an effort to reduce cell toxicity caused by synaptic damage. Though not studied in dipterans, patch clamp analyses of rat Nav1.8 sodium channel isoform in Xenopus oocytes indicated that membrane depolarization was not significantly impacted by permethrin or deltamethrin treatment, however there were significant decreases in calcium influx and glutamate release that were similar for both permethrin and deltamethrin [30]. Thus, taken together with the evidence above, we hypothesize that, in addition to VGSC, channels/transporters involved in glutamate and calcium homeostasis are major drivers of selection of transport category genes that contribute to resistance.

The genetic associations identified here need to be more comprehensively validated in mosquitoes to further support their roles in pyrethroid resistance.

Materials and methods

Collections and bioassays

Both populations used for this work were from the Merida region of Yucatan, Mexico and were collected from locations that were 13 kilometers apart. Population #1 was from Viva Caucel (20°59'57.83"N 89°42'15.41"W) was collected in 2011 by María Alba Loroño-Pino and Julian García-Rejón of Universidad Autónoma de Yucatán. The original larval collection (n = 60) was obtained from ~20 containers; more details are reported in Saavedrea-Rodriguez et al. 2015 [7]. Progeny from this collection was reared two generations in the laboratory. The third generation (about 3000 eggs) was reared to adulthood; a subset (Table 1) was exposed to permethrin-coated bottles using an amount that approximated the LC50 (lethal concentration, 50%, 25 μg per bottle active ingredient, Chem Service, West Chester, PA) [8]. Bioassays consisted of treatment of ~50 female mosquitoes (3–4 days old) for one hour [45]. Active mosquitoes, which were not knocked down (resistant), were transferred to cardboard cups and frozen at -80ᵒC. Mosquitoes that were knocked down were transferred to a second cardboard cup and placed into an incubator at 28ᵒC and 70% humidity. After a four hour recovery period, mosquitoes that remained inactive were scored as susceptible. DNA was isolated from individual mosquitoes using a salt extraction method [91] and resuspended in 150 μl of TE buffer (10 mM Tris-HCl, 1 mM EDTA pH 8.0). We typically isolate DNA from individual mosquitoes for archival purposes. For the current work, individual mosquito DNAs were carefully measured to allow for equimolar amounts of each mosquito in the sequencing libraries (see below).

Population #2 was also from Viva Caucel and was treated in bottles coated with 3 μg deltamethrin, also the LC50, using the same scheme as described above. Population #3 was from Vergel, Yucatán State in Southern Mexico (20°57'30.08"N, 89°35'4.33"W). It was also collected in 2011 by Universidad Autónoma de Yucatán and had been through two generations in the laboratory. The original larval collection (n = 47) was sampled from ~20 container sites [7]. Progeny (about 3000 eggs) from the 3rd generation of this collection was reared to adulthood; a subset (Table 1) were treated with permethrin, using the same conditions, as described above for population #1.

Gene silencing/bioassays

A 500 base-pair region of CHL or NVY was synthesized and inserted into a bacterial plasmid, pUC-IDTkan(Integrated DNA technologies). dsRNAs were generated using T7 primers appended to gene-specific primers, using the methods of Campbell et al. [92], wherein T7 dsDNA was prepared from synthesized gene fragments, and dsRNA was subsequently transcribed from the dsDNA. All primers are listed in S3 Table. Three to seven day old Vergel colony mosquitoes (F24) were injected with ~400 ng dsRNA and assayed at 3 days post-injection. Gene silencing was validated by qRT-PCR. Briefly, total RNA was extracted from pools of 5 mosquitoes (Qiagen RNAEasy kit), and 10 ng RNA was used as input into each 20 μl qRT-PCR reaction (NovaScript, Qiagen). Data was compiled from at least three biological replicates and three technical replicates for each treatment group. Average values for amplified ribosomal protein S7 (RPS7) and aaeACTIN were used for normalizing quantitation of relative expression using the delta-delta Ct relative gene expression calculation, relative to β-Gal-dsRNA injected control mosquitoes. At 3 days-post-injection, pools of mosquitoes were subjected to standard CDC bottle assays [45] using the discriminating dose of 1.5 μg permethrin per bottle (Sigma-Aldrich). Because the colony had been maintained for 24 generations in the lab without permethrin selection, the population had regained substantial sensitivity to insecticide treatment. Knockdown was measured in 10 minute increments for one hour, and mortality was scored at 24 hours post-treatment. Knockdown data represent a compilation of 5 experimental replicates.

Pooling and quantification of samples

For each location/treatment group, gDNA from 25 resistant females was pooled for each of two replicates; similarly gDNA from 2 pools of 25 susceptible females was prepared. Before pooling, DNA from each individual mosquito was quantified using the Quant-IT Pico Green kit (Life Technologies) and ~40 ng from each individual DNA sample (25 individuals per library) was used for a final DNA pool of 1 μg. A Covaris S2 sonicator (Covaris Ltd, Brighton UK) sheared pooled DNA to an average size of 500 bp. Sonication conditions were: duty cycle 10%, Intensity 5.0, Cycles per burst 200, Duration 40 seconds, Mode Frequency sweeping, Displayed Power 23W, Temperature 5.5° - 6°C. Each TruSeq DNA LT (v.2) library was prepared using 1 μg of sheared genomic DNA following manufacturer’s recommendations. We prepared one library for each of the twelve DNA pools following the Illumina TruSeq DNA Sample preparation guide (Illumina, San Diego CA).

Equimolar quantities of prepared libraries were pooled and enriched for coding sequences by exome capture using custom SeqCap EZ Developer probes (Nimblegen) [19]. In total, 26.7Mb of the genome (2%) was targeted for enrichment, as described elsewhere [20]. Overlapping probes covering the protein coding sequence (not including UTRs) in the AaegL1.3 gene annotations (https://www.vectorbase.org/organisms/aedes-aegypti/liverpool-lvp/AaegL1.3) were produced by Nimblegen. Enrichment followed the Nimblegen SeqCap EZ protocol. Briefly, pooled TruSeq libraries were hybridized to the probes for 64 hours at 47°C, unbound DNA was washed away, and the targeted DNA was eluted and amplified. These were then sequenced on 2 lanes of a HiSeq2500 (Illumina) for paired-end 2 x 100 nt sequencing by the Centers for Disease Control and Prevention and Control in Atlanta, GA, producing reads with quality scores > 30.

Bioinformatics

Alignments and population genetics pipeline.

The analytical methods reported here have been used previously to identify polymorphic differences between aedine subspecies and sex-specific polymorphisms [19,93]. Though the published versions of the papers reported FST values rather than contingency χ2 –log10(p values), preliminary analyses of–log10(p values) showed trends identical to those in the published reports. Importantly, the sex-specific polymorphisms we described were also replicated independently by others [94]. All raw reads were trimmed of adapters and filtered using cutadapt [95]. The AaegL5 genome build [21] of 18,081 coding and non-coding genes was used, including all 5’UTRs, exons, introns, 3’UTRs. The 5’ and 3’ non-transcribed regions in previously reported alignments were excluded [19]. Individual replicate trimmed fastq files were aligned to a custom reference file using GSNAP (version 2013-10-28), allowing 10% divergence [96]. GSNAP outputs were converted to *.mpileup files, using SAMtools [97]. The “readcounts” command in Varscan2 (v2.3.5) [98] was used to convert *.mpileup files to readcounts output, using the following options:—min-coverage 15—min-base-qual 30. FORTRAN programs used in this analysis are available upon request. The number of aligned reads was determined in .bam files, using SAMtools flagstat and multiplied by the read-length (100 nt) to achieve the total nucleotides aligned. The ratio of variant sites per chromosome was multiplied by 1000 and divided by the total aligned nucleotides (S1 Table) to obtain the “Ratio sites/aligned nucleotide * 1000.”

Because the probes were designed using genome build L1.3, and the current reference genome is L5, we checked one of the libraries (VCDA-rep1) to determine the percentage and depth of sequencing for capture target sequences that aligned to L5. Coverage depth files of .bam files were generated using SAMtools, then the analogous L5 capture coordinate intervals were located. About 15.3% (chr 1), 6.5% (chr 2), and 11.7% (chr 3) of L5 capture targets were not either not sequenced or had coverage less than the 15 read coverage cut-off. In addition, a blastn search of L5 transcripts using the capture probes indicated that they cover (13,942/18,081) 77% of the genes in L5.

All sequencing data generated under this project are available at the National Center for Biotechnology Information (NCBI) Sequence Read Archive, Bioproject accession number PRJNA393171. VCPA permethrin-treated resistant replicates 1 and 2 are SRR5805471 and SRR5805472, respectively. VCPD permethrin-treated susceptible replicates 1 and 2 are SRR5805473 and SRR5805470, respectively. VCDA deltamethrin-treated resistant replicates 1 and 2 are SRR5805467 and SRR5805466, respectively. VCDD deltamethrin-treated susceptible replicates 1 and 2 are SRR5805469 and SRR5805468, respectively. VePA permethrin-treated resistant replicates 1 and 2 are SRR5805465 and SRR5805464, respectively. VePD permethrin-treated susceptible replicates 1 and 2 are SRR5805475 and SRR5805474, respectively.

In-house FORTRAN scripts were used to calculate heterogeneity χ2 values with degrees of freedom equal to the number of alternate nucleotides minus 1. The contingency χ2 calculation tested for the presence of an alternate nucleotide at each position along the genome and computed a -log10(p value) at each SNP [19]. For the purpose of highlighting genes with high statistical association to pyrethroid resistance, calculations of weighted gene-wise average scores were made by determining the mean of the top 5% -log10(p value) values per gene. Our rationale for this calculation was as follows. Along the length of a given coding sequence, many nucleotide sites showed little change relative to the reference sequence, thus resulting in many -log10(p value)at or near zero for a given SNP. Coding changes indicative of selection toward pyrethroid resistance may be present in a small subset of sites relative to the full length gene. Therefore, a weighted average allowed us to identify those genes with localized high value SNPs. Statistical significance was determined using the Benjamini-Hochberg False Discovery Rate [22] (FDR, alpha = 0.01) and a cut-off of 4.0 was established for gene-wise averages.

Gene annotations were obtained from Vectorbase [23] using BioMart. Functional groups were assigned to the following categories: ‘apoptosis’(APOP); ‘chemosensory response‘ (CSR); cytoskeletal/structural (CYT/STR); ‘diverse’ (DIV) for genes with multiple or less clearly defined function(s); ‘mitochondrial’ (MIT) for gene products localized to mitochondrial compartment, regardless of function; (LIPID) for genes with predicted function in lipid processing or biosynthesis; (ReDox) for genes with predicted function in oxidation/reduction processes; proteolysis or proteosomal activity (PROT); replication/(DNA)repair/transcription /translation’ (RRTT); signal transduction (SigT); metabolism (MET); ‘unknown’ (UNK) for uncharacterized genes. Finally, the ‘transport’ (TRP) category included all gene products predicted to be involved in moving molecules across membranes, including receptors, exclusive of secondary messengers and signaling receptors. For GSEA, the mitochondrial and redox categories were combined, as several key redox activities occur in proteins localized to mitochondria.

Statistics.

Descriptive statistics were calculated in R (version 3.0.2). The ratio of variant sites per nucleotide of aligned reads was calculated as follows: binomial probability distributions were performed and hypergeometric tests were calculated; each functional subset was compared to the corresponding group within the set of all coding genes (n = 14,626). This analysis allowed us to identify the gene sets that were enriched compared to similar categories among all coding genes.

Supporting information

[Figure omitted. See PDF.]

S1 Table. Sequencing details for all libraries.

Number of reads aligned to reference, percent trimmed reads mapped to reference, number of variant sites per chromosome, ratio of sites per aligned nucleotide * 1000.

https://doi.org/10.1371/journal.pone.0211497.s001

(XLSX)

S2 Table. χ2 contingency values for all collections in the high association set (n = 1053).

Genes are arranged in order of their physical position. Genes in the top ten percent of all common genes are highlighted in gray, those with predicted motoneuron synaptic function are highlighted in pink, and predicted transcriptional regulators are highlighted in green. RefSeq locus, chromosome, function, aliases, description, genomic_nucleotide_accession, start_position, end_position, # exons. For each collection- “_pvalue”, weighted average–log(p value) per gene; AvgHexp_Resistant; AvgHexp_Susceptible.

https://doi.org/10.1371/journal.pone.0211497.s002

(XLSX)

S3 Table. Primers used in gene silencing and qRT-PCR.

https://doi.org/10.1371/journal.pone.0211497.s003

(XLSX)

S1 Fig. CHL- and NVY-silenced mosquitoes show increased susceptibility to permethrin.

Mosquitoes were injected with phosphate-buffered saline (PBS, n = 109), βGal-dsRNA (n = 113), chl-dsRNA (n = 172), nvy-dsRNA (n = 129) or left untreated (n = 136). At 3 dpt, each replicate was subjected to ~1.5 ug permethrin in a CDC bottle assay; knockdown was recorded at 10 minute intervals. Left Y axis indicates percent knockdown from 10–60 minutes. Right Y-axis indicates percent mortality at 24 hours. Error bars indicate SEM. Data represent a compilation of 4 to 5 biological replicates.

https://doi.org/10.1371/journal.pone.0211497.s004

(TIF)

Acknowledgments

We thank María Alba Loroño-Pino and Julian García-Rejón of Universidad Autónoma de Yucatán, México for providing mosquito material. We would also like to thank Dr. Venkatachalam Udhayakumar for his valuable assistance.

Disclaimer

The findings and conclusions in this paper are those of the authors and do not necessarily represent the official position of the Centers for Disease Control and Prevention.

Citation: Campbell CL, Saavedra-Rodriguez K, Kubik TD, Lenhart A, Lozano-Fuentes S, Black WC IV (2019) Vgsc-interacting proteins are genetically associated with pyrethroid resistance in Aedes aegypti. PLoS ONE 14(1): e0211497. https://doi.org/10.1371/journal.pone.0211497

References

1. Bhatt S, Gething PW, Brady OJ, Messina JP, Farlow AW, Moyes CL, et al. (2013) The global distribution and burden of dengue. Nature 496: 504–507. pmid:23563266

2. Boeuf P, Drummer HE, Richards JS, Scoullar MJ, Beeson JG (2016) The global threat of Zika virus to pregnancy: epidemiology, clinical perspectives, mechanisms, and impact. BMC Med 14: 112. pmid:27487767

3. Lorono-Pino MA, Farfan-Ale JA, Zapata-Peraza AL, Rosado-Paredes EP, Flores-Flores LF, Garcia-Rejon JE, et al. (2004) Introduction of the American/Asian genotype of dengue 2 virus into the Yucatan State of Mexico. Am J Trop Med Hyg 71: 485–492. pmid:15516647

4. NOM (2010) NOM-032-SSA2-2010 para la vigilancia epidemiológica,prevención y control de enfermedades transmitidas por vector. Mexico City: Norma Oficial Mexicana.

5. Du Y, Nomura Y, Satar G, Hu Z, Nauen R, He SY, et al. (2013) Molecular evidence for dual pyrethroid-receptor sites on a mosquito sodium channel. Proc Natl Acad Sci U S A 110: 11785–11790. pmid:23821746

6. Harris AF, Rajatileka S, Ranson H (2010) Pyrethroid resistance in Aedes aegypti from Grand Cayman. Am J Trop Med Hyg 83: 277–284. pmid:20682868

7. Saavedra-Rodriguez K, Beaty M, Lozano-Fuentes S, Denham S, Garcia-Rejon J, Reyes-Solis G, et al. (2015) Local evolution of pyrethroid resistance offsets gene flow among Aedes aegypti collections in Yucatan State, Mexico. Am J Trop Med Hyg 92: 201–209. pmid:25371186

8. Saavedra-Rodriguez K, Maloof FV, Campbell CL, Garcia-Rejon J, Lenhart A, Penilla P, et al. (2018) Parallel evolution of vgsc mutations at domains IS6, IIS6 and IIIS6 in pyrethroid resistant Aedes aegypti from Mexico. Sci Rep 8: 6747. pmid:29712956

9. Vera-Maloof FZ, Saavedra-Rodriguez K, Elizondo-Quiroga AE, Lozano-Fuentes S, Black IV WC (2015) Coevolution of the Ile1,016 and Cys1,534 Mutations in the Voltage Gated Sodium Channel Gene of Aedes aegypti in Mexico. PLoS Negl Trop Dis 9: e0004263. pmid:26658798

10. David JP, Faucon F, Chandor-Proust A, Poupardin R, Riaz MA, Bonin A, et al. (2014) Comparative analysis of response to selection with three insecticides in the dengue mosquito Aedes aegypti using mRNA sequencing. BMC Genomics 15: 174. pmid:24593293

11. Faucon F, Dusfour I, Gaude T, Navratil V, Boyer F, Chandre F, et al. (2015) Identifying genomic changes associated with insecticide resistance in the dengue mosquito Aedes aegypti by deep targeted sequencing. Genome Research 25: 1347–1359. pmid:26206155

12. Saavedra-Rodriguez K, Strode C, Flores Suarez A, Fernandez Salas I, Ranson H, Hemingway J, et al. (2008) Quantitative trait loci mapping of genome regions controlling permethrin resistance in the mosquito Aedes aegypti. Genetics 180: 1137–1152. pmid:18723882

13. Saavedra-Rodriguez K, Suarez AF, Salas IF, Strode C, Ranson H, Hemingway J, et al. (2012) Transcription of detoxification genes after permethrin selection in the mosquito Aedes aegypti. Insect Molecular Biology 21: 61–77. pmid:22032702

14. Riveron JM, Yunta C, Ibrahim SS, Djouaka R, Irving H, Menze BD, et al. (2014) A single mutation in the GSTe2 gene allows tracking of metabolically based insecticide resistance in a major malaria vector. Genome Biology 15: R27. pmid:24565444

15. Saavedra-Rodriguez K, Strode C, Flores AE, Garcia-Luna S, Reyes-Solis G, Ranson H, et al. (2014) Differential transcription profiles in Aedes aegypti detoxification genes after temephos selection. Insect Mol Biol 23: 199–215. pmid:24299217

16. Terhzaz S, Cabrero P, Brinzer RA, Halberg KA, Dow JA, Davies SA (2015) A novel role of Drosophila cytochrome P450-4e3 in permethrin insecticide tolerance. Insect Biochem Mol Biol 67: 38–46. pmid:26073628

17. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. (2005) Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A 102: 15545–15550. pmid:16199517

18. Hamid PH, Prastowo J, Ghiffari A, Taubert A, Hermosilla C (2017) Aedes aegypti resistance development to commonly used insecticides in Jakarta, Indonesia. PLoS One 12: e0189680. pmid:29253003

19. Dickson LB, Campbell CL, Juneja P, Jiggins FM, Sylla M, Black WCt (2017) Exon-enriched libraries reveal large genic differences between Aedes aegypti from Senegal, West Africa, and populations outside Africa. G3 (Bethesda) 7: 571–582.

20. Juneja P, Ariani CV, Ho YS, Akorli J, Palmer WJ, Pain A, et al. (2015) Exome and transcriptome sequencing of Aedes aegypti identifies a locus that confers resistance to Brugia malayi and alters the immune response. PLoS Pathog 11: e1004765. pmid:25815506

21. Matthews BJ, Dudchenko O, Kingan SB, Koren S, Antoshechkin I, Crawford JE, et al. (2018) Improved reference genome of Aedes aegypti informs arbovirus vector control. Nature 563: 501–507. pmid:30429615

22. Benjamini Y, Hochberg Y (1995) Controlling the false discovery rate—a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society Series B-Methodological 57: 289–300.

23. Lawson D, Arensburger P, Atkinson P, Besansky NJ, Bruggner RV, Butler R, et al. (2007) VectorBase: a home for invertebrate vectors of human pathogens (http://www.vectorbase.org/). Nucleic Acids Res. pp. D503–505. pmid:17145709

24. Timoshevskiy VA, Severson DW, Debruyn BS, Black WC, Sharakhov IV, Sharakhova MV (2013) An integrated linkage, chromosome, and genome map for the yellow fever mosquito Aedes aegypti. PLoS Negl Trop Dis 7: e2052. pmid:23459230

25. Tsuda L, Kaido M, Lim YM, Kato K, Aigaki T, Hayashi S (2006) An NRSF/REST-like repressor downstream of Ebi/SMRTER/Su(H) regulates eye development in Drosophila. The EMBO Journal 25: 3191–3202. pmid:16763555

26. Howlett IC, Rusan ZM, Parker L, Tanouye MA (2013) Drosophila as a model for intractable epilepsy: gilgamesh suppresses seizures in para(bss1) heterozygote flies. G3 (Bethesda) 3: 1399–1407.

27. Lin WH, Giachello CN, Baines RA (2017) Seizure control through genetic and pharmacological manipulation of Pumilio in Drosophila: a key component of neuronal homeostasis. Dis Model Mech 10: 141–150. pmid:28067623

28. Shao D, Okuse K, Djamgoz MB (2009) Protein-protein interactions involving voltage-gated sodium channels: Post-translational regulation, intracellular trafficking and functional expression. Int J Biochem Cell Biol 41: 1471–1481. pmid:19401147

29. Streit AK, Fan YN, Masullo L, Baines RA (2016) Calcium Imaging of Neuronal Activity in Drosophila Can Identify Anticonvulsive Compounds. PLoS One 11: e0148461. pmid:26863447

30. Breckenridge CB, Holden L, Sturgess N, Weiner M, Sheets L, Sargent D, et al. (2009) Evidence for a separate mechanism of toxicity for the Type I and the Type II pyrethroid insecticides. Neurotoxicology 30 Suppl 1: S17–31.

31. Jenson LJ, Anderson TD, Bloomquist JR (2016) Insecticide sensitivity of native chloride and sodium channels in a mosquito cell line. Pestic Biochem Physiol 130: 59–64. pmid:27155485

32. Soderlund DM, Bloomquist JR (1989) Neurotoxic actions of pyrethroid insecticides. Annu Rev Entomol 34: 77–96. pmid:2539040

33. Jenkins SM, Bennett V (2001) Ankyrin-G coordinates assembly of the spectrin-based membrane skeleton, voltage-gated sodium channels, and L1 CAMs at Purkinje neuron initial segments. J Cell Biol 155: 739–746. pmid:11724816

34. Sedgwick SG, Smerdon SJ (1999) The ankyrin repeat: a diversity of interactions on a common structural framework. Trends Biochem Sci 24: 311–316. pmid:10431175

35. Pielage J, Fetter RD, Davis GW (2006) A postsynaptic spectrin scaffold defines active zone size, spacing, and efficacy at the Drosophila neuromuscular junction. J Cell Biol 175: 491–503. pmid:17088429

36. Zhang W, Cheng LE, Kittelmann M, Li J, Petkovic M, Cheng T, et al. (2015) Ankyrin Repeats Convey Force to Gate the NOMPC Mechanotransduction Channel. Cell 162: 1391–1403. pmid:26359990

37. Fotia AB, Ekberg J, Adams DJ, Cook DI, Poronnik P, Kumar S (2004) Regulation of neuronal voltage-gated sodium channels by the ubiquitin-protein ligases Nedd4 and Nedd4-2. J Biol Chem 279: 28930–28935. pmid:15123669

38. van Bemmelen MX, Rougier JS, Gavillet B, Apotheloz F, Daidie D, Tateyama M, et al. (2004) Cardiac voltage-gated sodium channel Nav1.5 is regulated by Nedd4-2 mediated ubiquitination. Circ Res 95: 284–291. pmid:15217910

39. Ingham RJ, Gish G, Pawson T (2004) The Nedd4 family of E3 ubiquitin ligases: functional diversity within a common modular architecture. Oncogene 23: 1972–1984. pmid:15021885

40. Rougier JS, van Bemmelen MX, Bruce MC, Jespersen T, Gavillet B, Apotheloz F, et al. (2005) Molecular determinants of voltage-gated sodium channel regulation by the Nedd4/Nedd4-like proteins. Am J Physiol Cell Physiol 288: C692–701. pmid:15548568

41. Cho RJ, Huang M, Campbell MJ, Dong H, Steinmetz L, Sapinoso L, et al. (2001) Transcriptional regulation and function during the human cell cycle. Nat Genet 27: 48–54. pmid:11137997

42. Feinstein PG, Kornfeld K, Hogness DS, Mann RS (1995) Identification of homeotic target genes in Drosophila melanogaster including nervy, a proto-oncogene homologue. Genetics 140: 573–586. pmid:7498738

43. Terman JR, Kolodkin AL (2004) Nervy links protein kinase a to plexin-mediated semaphorin repulsion. Science 303: 1204–1207. pmid:14976319

44. Wildonger J, Mann RS (2005) Evidence that nervy, the Drosophila homolog of ETO/MTG8, promotes mechanosensory organ development by enhancing Notch signaling. Dev Biol 286: 507–520. pmid:16168983

45. Brogdon WG, McAllister JC (1998) Simplification of adult mosquito bioassays through use of time-mortality determinations in glass bottles. J Am Mosq Control Assoc 14: 159–164. pmid:9673916

46. Otali D, Novak RJ, Wan W, Bu S, Moellering DR, De Luca M (2014) Increased production of mitochondrial reactive oxygen species and reduced adult life span in an insecticide-resistant strain of Anopheles gambiae. Bull Entomol Res 104: 323–333. pmid:24555527

47. Pavlidi N, Monastirioti M, Daborn P, Livadaras I, Van Leeuwen T, Vontas J (2012) Transgenic expression of the Aedes aegypti CYP9J28 confers pyrethroid resistance in Drosophila melanogaster. Pesticide Biochemistry and Physiology 104: 132–135.

48. Riveron JM, Irving H, Ndula M, Barnes KG, Ibrahim SS, Paine MJI, et al. (2013) Directionally selected cytochrome P450 alleles are driving the spread of pyrethroid resistance in the major malaria vector Anopheles funestus. Proceedings of the National Academy of Sciences of the United States of America 110: 252–257. pmid:23248325

49. Vontas J, Kioulos E, Pavlidi N, Morou N, della Torre A, Ranson H (2012) Insecticide resistance in the major dengue vectors Aedes albopictus and Aedes aegypti. Pesticide Biochemistry and Physiology 104: 126–131.

50. Chandor-Proust A, Bibby J, Regent-Kloeckner M, Roux J, Guittard-Crilat E, Poupardin R, et al. (2013) The central role of mosquito cytochrome P450 CYP6Zs in insecticide detoxification revealed by functional expression and structural modelling. Biochemical Journal 455: 75–85. pmid:23844938

51. Zerba E (1988) Insecticidal activity of pyrethroids on insects of medical importance. Parasitol Today 4: S3–7. pmid:15463086

52. Maestre-Serrano R, Gomez-Camargo D, Ponce-Garcia G, Flores AE (2014) Susceptibility to insecticides and resistance mechanisms in Aedes aegypti from the Colombian Caribbean Region. Pestic Biochem Physiol 116: 63–73. pmid:25454522

53. Cui F, Li MX, Chang HJ, Mao Y, Zhang HY, Lu LX, et al. (2015) Carboxylesterase-mediated insecticide resistance: Quantitative increase induces broader metabolic resistance than qualitative change. Pestic Biochem Physiol 121: 88–96. pmid:26047115

54. Fonseca-Gonzalez I, Quinones ML, Lenhart A, Brogdon WG (2011) Insecticide resistance status of Aedes aegypti (L.) from Colombia. Pest Manag Sci 67: 430–437. pmid:21394876

55. Francis S, Saavedra-Rodriguez K, Perera R, Paine M, Black WCt, Delgoda R (2017) Insecticide resistance to permethrin and malathion and associated mechanisms in Aedes aegypti mosquitoes from St. Andrew Jamaica. PLoS One 12: e0179673. pmid:28650966

56. Wongtrakul J, Pongjaroenkit S, Leelapat P, Nachaiwieng W, Prapanthadara LA, Ketterman AJ (2010) Expression and Characterization of Three New Glutathione Transferases, an Epsilon (AcGSTE2-2), Omega (AcGSTO1-1), and Theta (AeGSTT1-1) From Anopheles cracens (Diptera: Culicidae), a Major Thai Malaria Vector. Journal of Medical Entomology 47: 162–171. pmid:20380296

57. Brinzer RA, Henderson L, Marchiondo AA, Woods DJ, Davies SA, Dow JA (2015) Metabolomic profiling of permethrin-treated Drosophila melanogaster identifies a role for tryptophan catabolism in insecticide survival. Insect Biochem Mol Biol 67: 74–86. pmid:26474926

58. Corrigan JJ, Kearns CW (1963) Amino acid metabolism in DDT poisoned American roaches. Journal of Insect Physiology 9: 1–12.

59. Mansingh A (1965) The effect of malathion on the metabolism of amino acids in the German cockroach Blattella germanica. J Insect Physiol 11: 1389–1400. pmid:5829268

60. Grigoraki L, Lagnel J, Kioulos I, Kampouraki A, Morou E, Labbe P, et al. (2015) Transcriptome Profiling and Genetic Study Reveal Amplified Carboxylesterase Genes Implicated in Temephos Resistance, in the Asian Tiger Mosquito Aedes albopictus. PLoS Negl Trop Dis 9: e0003771. pmid:26000638

61. Julio AH, Gigliolli AA, Cardoso KA, Drosdoski SD, Kulza RA, Seixas FA, et al. (2017) Multiple resistance to pirimiphos-methyl and bifenthrin in Tribolium castaneum involves the activity of lipases, esterases, and laccase2. Comp Biochem Physiol C Toxicol Pharmacol 195: 27–43. pmid:28163254

62. Feng G, Marion JR, Clark JM (1992) Suppression of pyrethroid-dependent neurotransmitter release from synaptosomes of knockdown-resistant house flies under pulsed-depolarization conditions during continuous perfusion. Pesticide Biochemistry and Physiology 42: 64–77.

63. Badiou A, Meled M, Belzunces LP (2008) Honeybee Apis mellifera acetylcholinesterase—a biomarker to detect deltamethrin exposure. Ecotoxicol Environ Saf 69: 246–253. pmid:17215041

64. Hannun YA, Obeid LM (2002) The Ceramide-centric universe of lipid-mediated cell regulation: stress encounters of the lipid kind. J Biol Chem 277: 25847–25850. pmid:12011103

65. Zhao L, Pridgeon JW, Becnel JJ, Clark GG, Linthicum KJ (2009) Mitochondrial gene cytochrome b developmental and environmental expression in Aedes aegypti (Diptera: Culicidae). J Med Entomol 46: 1361–1369. pmid:19960681

66. Garcia GP, Flores AE, Fernandez-Salas I, Saavedra-Rodriguez K, Reyes-Solis G, Lozano-Fuentes S, et al. (2009) Recent rapid rise of a permethrin knock down resistance allele in Aedes aegypti in Mexico. PLoS Negl Trop Dis 3: e531. pmid:19829709

67. Jones CM, Liyanapathirana M, Agossa FR, Weetman D, Ranson H, Donnelly MJ, et al. (2012) Footprints of positive selection associated with a mutation (N1575Y) in the voltage-gated sodium channel of Anopheles gambiae. Proc Natl Acad Sci U S A 109: 6614–6619. pmid:22493253

68. Moyes CL, Vontas J, Martins AJ, Ng LC, Koou SY, Dusfour I, et al. (2017) Contemporary status of insecticide resistance in the major Aedes vectors of arboviruses infecting humans. PLoS Neglected Tropical Diseases 11: e0005625. pmid:28727779

69. O'Reilly AO, Khambay BP, Williamson MS, Field LM, Wallace BA, Davies TG (2006) Modelling insecticide-binding sites in the voltage-gated sodium channel. Biochem J 396: 255–263. pmid:16475981

70. Yan Y, Yang Y, You J, Yang G, Xu Y, Huang N, et al. (2011) Permethrin modulates cholinergic mini-synaptic currents by partially blocking the calcium channel. Toxicol Lett 201: 258–263. pmid:21251955

71. Chen M, Du Y, Nomura Y, Zhu G, Zhorov BS, Dong K (2017) Mutations of two acidic residues at the cytoplasmic end of segment IIIS6 of an insect sodium channel have distinct effects on pyrethroid resistance. Insect Biochem Mol Biol 82: 1–10. pmid:28111191

72. Ganetzky B, Wu CF (1986) Neurogenetics of membrane excitability in Drosophila. Annu Rev Genet 20: 13–44. pmid:2434020

73. Harms E, Chu T, Henrion G, Strickland S (2000) The only function of Grauzone required for Drosophila oocyte meiosis is transcriptional activation of the cortex gene. Genetics 155: 1831–1839. pmid:10924478

74. Menon KP, Carrillo RA, Zinn K (2013) Development and plasticity of the Drosophila larval neuromuscular junction. Wiley Interdiscip Rev Dev Biol 2: 647–670. pmid:24014452

75. Abe T, Yamazaki D, Murakami S, Hiroi M, Nitta Y, Maeyama Y, et al. (2014) The NAV2 homolog Sickie regulates F-actin-mediated axonal growth in Drosophila mushroom body neurons via the non-canonical Rac-Cofilin pathway. Development 141: 4716–4728. pmid:25411210

76. He J, Zhou R, Wu Z, Carrasco MA, Kurshan PT, Farley JE, et al. (2016) Prevalent presence of periodic actin-spectrin-based membrane skeleton in a broad range of neuronal cell types and animal species. Proc Natl Acad Sci U S A 113: 6029–6034. pmid:27162329

77. Li T, Cao C, Yang T, Zhang L, He L, Xi Z, et al. (2015) A G-protein-coupled receptor regulation pathway in cytochrome P450-mediated permethrin-resistance in mosquitoes, Culex quinquefasciatus. Sci Rep 5: 17772. pmid:26656663

78. Lynn BD, Li X, Nagy JI (2012) Under construction: building the macromolecular superstructure and signaling components of an electrical synapse. J Membr Biol 245: 303–317. pmid:22722764

79. Arora D, Siddiqui MH, Sharma PK, Singh SP, Tripathi A, Mandal P, et al. (2016) Evaluation and physiological correlation of plasma proteomic fingerprints for deltamethrin-induced hepatotoxicity in Wistar rats. Life Sci 160: 72–83. pmid:27142831

80. Nix P, Hammarlund M, Hauth L, Lachnit M, Jorgensen EM, Bastiani M (2014) Axon regeneration genes identified by RNAi screening in C. elegans. J Neurosci 34: 629–645. pmid:24403161

81. Zhang YV, Hannan SB, Stapper ZA, Kern JV, Jahn TR, Rasse TM (2016) The Drosophila KIF1A Homolog unc-104 Is Important for Site-Specific Synapse Maturation. Front Cell Neurosci 10: 207. pmid:27656128

82. Gordon S, Dickinson MH (2006) Role of calcium in the regulation of mechanical power in insect flight. Proc Natl Acad Sci U S A 103: 4311–4315. pmid:16537527

83. Mee CJ, Pym EC, Moffat KG, Baines RA (2004) Regulation of neuronal excitability through pumilio-dependent control of a sodium channel gene. J Neurosci 24: 8695–8703. pmid:15470135

84. Driscoll HE, Muraro NI, He M, Baines RA (2013) Pumilio-2 regulates translation of Nav1.6 to mediate homeostasis of membrane excitability. J Neurosci 33: 9644–9654. pmid:23739961

85. Muraro NI, Weston AJ, Gerber AP, Luschnig S, Moffat KG, Baines RA (2008) Pumilio binds para mRNA and requires Nanos and Brat to regulate sodium current in Drosophila motoneurons. J Neurosci 28: 2099–2109. pmid:18305244

86. Fernandez-Chacon R, Konigstorfer A, Gerber SH, Garcia J, Matos MF, Stevens CF, et al. (2001) Synaptotagmin I functions as a calcium regulator of release probability. Nature 410: 41–49. pmid:11242035

87. Sampo B, Tricaud N, Leveque C, Seagar M, Couraud F, Dargent B (2000) Direct interaction between synaptotagmin and the intracellular loop I-II of neuronal voltage-sensitive sodium channels. Proc Natl Acad Sci U S A 97: 3666–3671. pmid:10737807

88. Srinivasan Y, Elmer L, Davis J, Bennett V, Angelides K (1988) Ankyrin and spectrin associate with voltage-dependent sodium channels in brain. Nature 333: 177–180. pmid:2452986

89. Song JY, Ichtchenko K, Sudhof TC, Brose N (1999) Neuroligin 1 is a postsynaptic cell-adhesion molecule of excitatory synapses. Proc Natl Acad Sci U S A 96: 1100–1105. pmid:9927700

90. Poskanzer KE, Marek KW, Sweeney ST, Davis GW (2003) Synaptotagmin I is necessary for compensatory synaptic vesicle endocytosis in vivo. Nature 426: 559–563. pmid:14634669

91. Black WC, DuTeau NM (1997) RAPD-PCR and SSCP analysis for insect population genetic studies. In: Crampton CBBaCL J., editor. The Molecular Biology of Insect Disease Vectors: A Methods Manual New York: Chapman and Hall. pp. 361–373.

92. Campbell CL, Keene KM, Brackney DE, Olson KE, Blair CD, Wilusz J, et al. (2008) Aedes aegypti uses RNA interference in defense against Sindbis virus infection. BMC Microbiol 8: 47. pmid:18366655

93. Campbell CL, Dickson LB, Lozano-Fuentes S, Juneja P, Jiggins FM, Black WCt (2017) Alternative patterns of sex chromosome differentiation in Aedes aegypti (L). BMC Genomics 18: 943. pmid:29202694

94. Fontaine A, Filipovic I, Fansiri T, Hoffmann AA, Cheng C, Kirkpatrick M, et al. (2017) Extensive Genetic Differentiation between Homomorphic Sex Chromosomes in the Mosquito Vector, Aedes aegypti. Genome Biol Evol 9: 2322–2335. pmid:28945882

95. Martin M (2011) Cutadapt Removes Adapter Sequences From High-Throughput Sequencing Reads. EMBnetjournal 17: 10–12.

96. Wu TD, Nacu S (2010) Fast and SNP-tolerant detection of complex variants and splicing in short reads. Bioinformatics 26: 873–881. pmid:20147302

97. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. (2009) The Sequence Alignment/Map format and SAMtools. Bioinformatics 25: 2078–2079. pmid:19505943

98. Koboldt DC, Zhang Q, Larson DE, Shen D, McLellan MD, Lin L, et al. (2012) VarScan 2: somatic mutation and copy number alteration discovery in cancer by exome sequencing. Genome Res 22: 568–576. pmid:22300766