Molecular basis of resistance to organophosphate insecticides in the New World screw-worm fly
Parasites & Vectors volume 13, Article number: 562 (2020)
The emergence of insecticide resistance is a fast-paced example of the evolutionary process of natural selection. In this study, we investigated the molecular basis of resistance in the myiasis-causing fly Cochliomyia hominivorax (Diptera: Calliphoridae) to dimethyl-organophosphate (OP) insecticides.
By sequencing the RNA from surviving larvae treated with dimethyl-OP (resistant condition) and non-treated larvae (control condition), we identified genes displaying condition-specific polymorphisms, as well as those differentially expressed.
Both analyses revealed that resistant individuals have altered expression and allele-specific expression of genes involved in proteolysis (specifically serine-endopeptidase), olfactory perception and cuticle metabolism, among others. We also confirmed that resistant individuals carry almost invariably the Trp251Ser mutation in the esterase E3, known to confer OP and Pyrethroid resistance. Interestingly, genes involved in metabolic and detoxifying processes (notably cytochrome P450s) were found under-expressed in resistant individuals. An exception to this were esterases, which were found up-regulated.
These observations suggest that reduced penetration and aversion to dimethyl-OP contaminated food may be important complementary strategies of resistant individuals. The specific genes and processes found are an important starting point for future functional studies. Their role in insecticide resistance merits consideration to better the current pest management strategies.
Pest organisms are defined as those which harm humans and human interests. Detrimental effects include disease, suffering and annoyance, reduction of livestock and agricultural yields, and quality loss of products among others. Pest populations have been almost exclusively controlled by pesticides. The overuse of pesticides, however, represents a strong selective pressure on natural populations leading to the selection of resistant individuals and to the emergence of resistant populations [1, 2]. The replacement of susceptible populations by resistant ones is of worldwide concern as it has greatly reduced the efficacy of pesticides, increasing the damages caused by pests.
The New World screw-worm fly, Cochliomyia hominivorax, an ectoparasite of mammals, is particularly troublesome. This species is one of the most important myiasis-causing flies of the neotropical region. Myiasis, an infestation of tissues of vertebrates by dipteran larvae , caused by C. hominivorax occurs when females oviposit in wounds or exposed tissues of an animal host. Even a small wound (e.g. the bite of a tick) is sufficient to attract a female . When the eggs hatch, larvae promptly feed on the live tissue of the host causing a series of deleterious effects.
Cochliomyia hominivorax particularly affects livestock and is responsible for severe economic losses to the cattle industry, estimated at US$336,62 million per year in Brazil alone . One of the main effects is the reduction of leather quality due to the scars caused by the larvae [5, 6], but they also include weight loss, and decrease in milk production. Infestation of calves is life-threatening .
Myiases due to C. hominivorax are also not rare in humans and generally occur in open wounds [8,9,10], mucosae [11,12,13,14,15] and particularly in people suffering from poor hygiene or neglect . Infestations in humans are an increasing issue  and introductions beyond the New World [18, 19] are a persisting concern.
Cochliomyia hominivorax is mainly controlled in Brazil by applying organophosphate (OP)-based insecticides . They are applied directly in the infested wound, targeting larvae feeding off the living tissue of the host. This compound targets the molecule acetylcholinesterase (AChE), which regulates nerve activity by breaking down the neurotransmitter acetylcholine (reviewed in ). The inactivation of AChE leads to an overstimulation and blockage of its receptors (e.g. muscle tissue), resulting in paralysis and insect death . Owing to the constant selective pressure from OP application during the last decades, resistant lineages have been strongly selected across several populations of C. hominivorax [1, 2].
Although C. hominivorax is responsible for great economic loss, few studies have investigated the molecular origin of the resistance phenotype. It is known that mutations at positions 137 and 251 in the enzyme esterase E3, present in C. hominivorax [1, 21, 22], enable the breakdown of the OP insecticide preventing the inactivation of AChE [23, 24]. The replacement of Glycine by Aspartic acid at position 137 (Gly137Asp mutation) is involved in diethyl OP (e.g. Diazinon) resistance  while the replacement of Tryptophan (Trp) by Serine (Ser) at position 251 (Trp251Ser mutation) confers resistance to dimethyl-OPs (e.g. malathion) , and possibly to pyrethroid insecticides . As a first attempt to uncover more candidate genes involved in the dimethyl-OP insecticide resistance, Carvalho et al.  sequenced larval and adult (female and male) transcriptomes of C. hominivorax and used quantitative PCR of 18 candidate genes coding for metabolic detoxification enzymes to test their association with the resistance phenotype. Of the genes analyzed, only Cyp6g1 was differentially expressed between control and resistant larvae.
In this study, we also conducted a transcriptome-wide analysis to understand the overall molecular basis of resistance in a C. hominivorax population with segregating resistant phenotype by using RNA-seq data. We measured and compared gene expression changes in a subset of individuals of this population selected after a treatment with high concentrations of dimethyl-OP (resistant condition) and non-treated individuals (control condition). We also identified condition-specific polymorphisms, detecting allele frequency shifts caused by the selection. We, therefore, tested the relative contribution of changes in gene expression and polymorphisms in the coding region to the resistant phenotype.
Cochliomyia hominivorax colony
In this study, we used a laboratory colony of C. hominivorax, which contained OP-resistant individuals segregating the Gly137Asp and Trp251Ser mutations in the esterase E3 gene [1, 21, 22], collected from Caiapônia, Goiás, Brazil on January 2005. The colony was maintained for 11 generations in laboratory conditions according to Carvalho et al. . The same resistant and control samples were used in this study and by Carvalho et al.  in their quantitative Polymerase Chain Reaction (qPCR) comparisons. Briefly, for the resistant condition, a sample from the laboratory population was treated with the dimethyl-OP insecticide dichlorvos (2,2-dichlorovinyl dimethyl phosphate; C4H7Cl2O4P; Fort Dodge, technical grade) at a lethal concentration (20 mg/L) for 90% of the population (LC90). The insecticide was directly mixed into the diet of the larvae consisting of fresh ground beef supplemented with bovine blood and water (2:1). A total of 500 second-instar larvae were maintained on the insecticide-containing diet for 24 h. Surviving individuals (resistant sample) were collected for total RNA extraction and subsequent analysis. Individuals of the control group were simply sampled from this laboratory population and maintained on the same diet without the insecticide. Individuals from the resistant and the control samples were collected at the same time for RNA extractions.
DNA and RNA extraction and sequencing
RNA extraction followed the procedure used by Carvalho et al. . Total RNA of resistant and control samples was isolated for each individual using TRIzol (InvitrogenTM, Thermo Fisher Scientific, Waltham, Massachusetts, USA) from whole bodies of 142 larvae: 44 and 58 surviving individuals respectively from the first and second replicates of the insecticide treatment, and 40 and 60 control larvae from the control group, first and second replicates, respectively. DNase I (InvitrogenTM, Thermo Fisher Scientific, Waltham, Massachusetts, USA) was used to remove genomic DNA contamination and mRNA-enriched samples were further purified by using Nucleospin RNA Clean-up columns (Macherey & Nagel, Düren, Germany). Quantification was performed using the fluorometer Qubit Quantitation Platform (InvitrogenTM, Thermo Fisher Scientific, Waltham, Massachusetts, USA).
RNA from single individuals of each replicate was pooled for sequencing. The sequencing of the first replicate (single-end reads) was outsourced to the University of Veterinary Medicine, Vienna, Austria and a second biological replicate (paired-end reads) to Laboratório Multiusuários Centralizado, ESALQ-USP (Piracicaba, Brazil). The RNA fragments resulting from the random breaking of the transcripts were converted into a complementary DNA (cDNA) library using the mRNA-Seq Sample Prep Kit (Illumina, San Diego, USA). The samples were sequenced using the HiSeq 2000 (Illumina, San Diego, USA).
For whole genome sequencing of the control population, high-quality DNA was obtained from seven-day-old pupae, as in this stage, no food remains in their gut. DNA extraction was carried out using a phenol/chlorophorm protocol . Samples were prepared with the Illumina TruSeq DNA Low Sample Protocol with the HiSeq SBS v4 High Output Kit (Illumina, San Diego, USA). The library was sequenced in the HiSeq 2500 platform to obtain paired-end (125 bp) DNA sequences. The sequencing service was provided by the Laboratório Multiusuários Centralizado, ESALQ-USP (Piracicaba, Brazil).
Preprocessing of the reads
To assess the quality of the reads, we used the program FastQC , which provided a quick overview of the C. hominivorax sequence files. To eliminate the low-quality region of the sequences, we used the program Trimmomatic . This program trimmed the sequences in order to remove read regions with an average quality lower than 15. After trimming, sequences shorter than 20 were discarded. We excluded identical reads to reduce the processing time of the sequences during assembly .
De novo transcriptome assembly and annotation
The transcriptome assembly was performed using the program Trinity  (version trinityrnaseq_r20140717, --normalize_reads --full_cleanup) in pair-end mode. The single-end reads of the first replicate were inputted as left reads.
The completeness of the assembly was assessed by using BUSCO  to search for complete single copy, complete duplicated, and fragmented orthologs within a Diptera database (“diptera_odb9”). The annotation of the transcriptome was carried out using FunctionAnnotator  (eukaryotic mode), which uses the NCBI-nr database to find the closest BLAST hits and Blast2GO to assign GO terms.
Identification of differentially expressed (DE) transcripts
For the differential expression analysis, we removed the redundancy in the transcriptome by clustering the assembled transcripts using CD-HIT-EST  (version 4.6, -c 0.95 -M 0 -T 0). This avoided discarding ambiguous read counts (multiple transcript mappings).
To estimate the expression levels of each transcript, we mapped separately the trimmed reads of each condition (Control and Resistant) to the non-redundant transcript set using Bowtie2  (version 2.2.3, ‘--local --very-sensitive-local’ for first single-end replicate and '--local --very-sensitive-local --maxins 1000 --no-mixed --no-discordant' for the second paired-end replicate).
We obtained a table of counts using eXpress  (v1.51) and RSEM  (v1.2.25) with default parameters. Both programs use an expectation maximization approach and do not require a reference genome for counting. As RSEM does not support gapped alignments, we ran Bowtie2 without the gap argument (“--sensitive --met-stderr --dpad 0 --gbar 99999999 --mp 1,1 --np 1 --score-min L,0,-0.1” and, additionally, “--no-mixed --no-discordant” for the pair-end libraries only).
To test for differential expression between the control and resistant conditions, we used two R packages of the Bioconductor repository; EdgeR  (v1.2.25) and DESeq2 . In both, counts were normalized by the library size. We considered for further examination the transcripts for which the false discovery rate (FDR) was below 0.05 and the absolute fold change higher than 2 in all analyses performed (eXpress or RSEM counting with either DESeq2 or EdgeR testing).
To test if the resistant phenotype could be a result of Single Nucleotide Polymorphisms (SNPs) in the transcripts, we conducted a variant analysis using the pipeline KisSplice  (version 2.4.0). This tool was designed to identify and test for differential allelic usage in RNA-seq data without the need of a reference genome. Briefly, kisSplice locally assembles the regions which surround variants (kissplice -s 1). The localisation of the polymorphic regions is then determined by aligning them to the transcriptome using BLAT  (version 34, -minIdentity=80). The R package KissDE was then used to identify which variants were enriched in one condition versus the other. Finally, we filtered the polymorphisms predicted to have a functional impact (non-synonymous SNPs within the coding region) by using KisSplice2reftranscriptome (version 1.2.2), which uses the longest open reading frames (ORFs) identified by Transdecoder  (version 3.0.1). For further analyses, we considered the condition specific, non-synonymous SNPs displaying (i) a magnitude of differential allelic expression of at least 20%, (ii) a coverage of at least 20 per replicate and (iii) an adjusted P-value below 0.01.
Because we found transcripts expressing predominantly one allele in the control condition (allele contributing to more than 90% of the total expression of that transcript), but expressing two alleles in the resistant condition, we used KisSplice to find SNPs in the genomic reads from the control condition, following the steps described above. From this result, we could determine if both alleles were indeed present in the control population.
Specifically for the esterase E3 (transcript c13624_g1_i1), we analysed the positions 137 and 251, which correspond to the respective positions of the two point mutations (Gly137Asp and Trp251Ser) of resistant individuals of C. hominivorax. For this, we counted the polymorphisms at those positions and compared their frequencies in each condition. We set a 1% cut-off to remove sequencing errors.
Gene Ontology (GO) analysis
We conducted a GO enrichment analysis to identify the GO terms over-represented in the gene sets impacted by the condition (either displaying a differential expression or containing condition-specific SNPs) versus those not impacted. This analysis was performed with the Fisher’s exact test, implemented in the program Blast2GO  (version 5.2.5), using a False Discovery Rate (FDR) threshold of 0.05. We tested separately the up- and down-regulated gene sets. The genes containing condition-specific SNPs were analysed together. The lists of over-represented GO terms were reduced to the most specific terms (version 5.2.5).
Transcriptome assembly and annotation
Transcriptome assembly was performed on trimmed and collapsed reads (Additional file 1: Table S1). Collapsing of identical reads significantly reduced computation time to build the transcriptome. The Trinity-assembled transcripts were clustered to reduce redundancy (see Methods). The final reduced transcriptome consisted of 33,038 transcripts (average length of 809.75bp, SD 1001.32 and N50 1493) of which 11,055 (33.5%) were assigned at least one GO term and 16,541 (50.0%) were annotated with a blast hit to the NCBI-nr database. The blast hits belonged mostly to Lucilia cuprina (12,055, 72.9%), Musca domestica (1411, 8.5%) and Stomoxys calcitrans (1209, 7.3%), which is concordant with the divergence times between C. hominivorax and each of these species .
Expression profiling comparison
To identify the differentially expressed (DE) transcripts between the control and the resistant conditions, we used two methods to calculate the abundance of each transcript (eXpress and RSEM) as well as two R packages to test the differential expression (DEseq2 and EdgeR, see Methods). All combinations of counting and testing methods resulted in a similar set of DE transcripts (Fig. 1). Only the common DE transcripts (142) found by all analyses were further examined. Of these, 77 were found up-regulated and 65 down-regulated in the resistant condition compared to the control condition.
Response to stress and immune-related transcripts
Consistent with previous studies [44,45,46,47,48,49,50,51], heat-shock proteins (responsible for the enriched term “response to stress” in the Gene Ontology enrichment analysis) and immunity effectors were found over-expressed in the resistant condition (Fig. 2, Additional file 2: Dataset S1).
Also consistent with previous studies (see mini-review from ), esterase genes, known for their detoxifying role (reviewed in ), were found up-regulated in the resistant larvae (transcripts c14015_g1_i1, c14015_g2_i3, c14015_g4_i2, c571_g1_i1; Fig. 2 and Additional file 2: Dataset S1). Additionally, we found the GO terms “transmembrane transport”, “positive regulation of sodium ion transport”, “voltage-gated sodium channel complex”, “voltage-gated sodium channel activity”, and “extracellular region” (Fig. 2 and Additional file 3: Dataset S2) over-represented in the up-regulated transcript set. One possibility is that the up-regulation of trans-membrane transporters could help remove toxins from the cells and could therefore represent a non-specific response to the insecticide. More specific to C. hominivorax’s response to dimethyl-OP, genes involved in proteolysis (mostly serine-type endopeptidases) were found up-regulated in the resistant larvae (Fig. 2, Additional file 2: Dataset S1). Consistent with this result, the GO terms “proteolysis” and “serine-type endopeptidase activity” were found over-represented in the GO enrichment analysis (Additional file 3: Dataset S2). Finally, an odorant binding protein (c9041_g1_i2) playing a role in the sensory perception of chemical stimulus was found more than 8 times more expressed in the resistant condition compared to the control (Fig. 2 and Additional file 3: Dataset S2). It is possible that the insecticide binds specifically to this receptor.
In contrast with this result, we found a number of genes involved in metabolic processes, including putative glutathione transferases (c6832_g1_i1, c6842_g2_i1), a tryptophan phenylalanine hydroxylase (henna) (c14100_g2_i5, c14100_g1_i3, c14100_g1_i1) and cytochrome P450s (c14453_g1_i2, c12908_g1_i1, c13224_g2_i1, c14195_g1_i1, c13812_g3_i2), down-regulated in the resistant condition (Fig. 2 and Additional file 2: Dataset S1). This result was in agreement with the GO analysis, in which we found the GO terms “metabolic process”, “oxidation-reduction process”, “L-phenylalanine catabolic process”, “drug catabolic process” and “small molecule metabolic process” over-represented in the down-regulated transcript set (Additional file 3: Dataset S2).
Altered expression was observed for structural constituents of the cuticle and genes involved in the chitin metabolic process with some transcripts found up-regulated (c12563_g1_i1, c8591_g1_i1, c10801_g1_i1, c954_g1_i1, c8201_g1_i1), while others were down-regulated (c12862_g1_i1, c8899_g1_i1, c10582_g1_i1, c8899_g2_i1, c7227_g2_i1), suggestive of modification to the cuticle’s composition and properties (Fig. 2 and Additional file 2: Dataset S1).
Besides expression difference, it is possible that the resistance to dimethyl-OP insecticides is due to condition-specific polymorphisms (differential allelic expression or selected point mutations). Indeed, the original population (control condition) had individuals carrying the mutations Gly137Asp and/or Trp251Ser in the esterase E3 gene known to confer resistance. In particular, we expected to find homozygous individuals for the Trp251Ser mutation in the resistant group. This mutation was shown to increase the resistance specifically to dimethyl-OP , the insecticide used in this study.
As expected and in accordance with the Esterase E3 genotyping of Carvalho et al. , the resistant individuals had predominantly the esterase E3 (c13624_g1_i1) 251Ser allele (Table 1), which confers resistance to dimethyl-OPs . We also detected a low frequency of the amino acid Leucine (Leu) at position 251 (251Leu allele) in resistant individuals (Table 1). The mutation Gly137Asp was found in the control condition (Table 1), however, the resistant individuals only presented the non-resistant 137Gly allele (Table 1). These observations led us to conclude that resistance-associated mutations at positions 137 and 251 are most likely not found on the same haplotype in this population (Fig. 3, lower panel). There was no evidence of recombination between these two sites in agreement with Carvalho et al. . It is to note that we detected a few 251Trp alleles in the resistant condition (Table 1).
With the present RNA-seq data, it was also possible to explore the expression of genes presenting condition-specific polymorphisms. Using the KisSplice  pipeline (See Methods), we identified 830 non-synonymous condition-specific polymorphisms within the coding region of 492 transcripts (1.5 SNPs/transcripts on average) (Additional file 4: Dataset S3). Of these transcripts, 256 (52%) displayed only one allele in the resistant sample. At present, it is unclear if in these cases, only one allele is expressed, or if individuals carrying the alternative allele died as a consequence of the insecticide treatment. Few transcripts (17) were found to be DE and have condition-specific SNPs (Table 2 and Additional file 5: Dataset S4). However, the GO enrichment analysis performed on the set of transcripts containing condition-specific SNPs revealed similar processes to the ones found in the DE GO analysis enrichment (Additional file 6: Dataset S5). In particular, we found several GO terms related to proteolysis/peptidase activity (“serine-type endopeptidase inhibitor activity”, “aminopeptidase activity”, “serine-type endopeptidase activity”, “negative regulation of endopeptidase activity”, “carboxypeptidase activity”) and metabolic processes (“chitin metabolic process”, “hormone metabolic process”, “cellular carbohydrate metabolic process”, “malate metabolic process”, “hexose metabolic process”) enriched in the transcripts with condition-specific SNPs. Also of interest, the GO terms “chitin binding” and “structural constituent of cuticle” were found over-represented. This result suggests that the processes related to the cuticle, metabolism and proteolysis are sensitive to dimethyl-OP insecticides in C. hominivorax and that some alleles are preferentially expressed depending on the condition or even necessary for the individual's survival when in contact with the insecticide.
Condition-specific polymorphisms were also found in transcripts potentially involved in olfactory behavior (c13758_g3_i4, c12483_g1_i2, c13091_g1_i2/4/9, c13818_g5_i3), pupation (c12606_g1_i1, c12483_g1_i2) and starvation (c12483_g1_i2, c12331_g1_i1), based on GO term assignment.
Of particular interest, the transcript c13758_g3_i4, corresponding to an odorant-binding protein (OBP), was found to contain five polymorphisms with frequency shifts between the resistant and control groups. This transcript, as well as the differentially expressed OBP c9041_g1_i2 mentioned earlier, could be involved in the perception of the insecticide.
For 227 transcripts, we found that the control condition predominantly expressed one allele (>90% of the gene expression was contributed from one allele) while the resistant condition expressed both. To validate these results, we generated 154,067,376 genomic paired sequences with an average length of 120 bp from the control population. Using these low coverage genomic reads, we could confirm that 105 of these transcripts (46,3%) were probably segregating in the original population, as both alleles were present in the DNA reads. In the RNA-seq data for these genes, one of the alleles is prevalent, while the genomic data from the same population shows both alleles. Hence, a differential allelic expression is the most plausible explanation, i.e. heterozygous individuals only expressed one of their two alleles. We hypothesize that the remaining transcripts (122) were also genetically heterozygous but had insufficient coverage to allow their detection in the genomic reads.
In this study, we compared the transcriptomes of second instar (L2) resistant and dimethyl-OP-treated larvae (resistant condition) to non-treated L2 larvae from a mixed population composed of susceptible and resistant individuals (control condition). The differential expression analysis revealed that only a few transcripts (142, 0.43%) had an altered expression between both conditions. In addition, we identified 492 transcripts with condition-specific polymorphisms. At this point, it is unclear if these polymorphisms are due to a differential allelic expression (part of the resistance would be due to phenotypic plasticity) or if individuals lacking certain polymorphisms died as a result of the OP-treatment (genotypic resistance). These transcripts were generally involved in similar processes, although only a few were found in both analyses.
From the differential allelic usage analysis, we could also confirm that the previously characterized Trp251Ser mutation in the esterase E3 gene [1, 21, 22] is necessary to confer resistance to dimethyl-OP insecticides (Fig. 3, lower panel), whereas the 137Asp mutation, that presumably confers resistance to dietyl-OP, was not present in the resistant individuals. These observations support the model that the 251Ser mutation is linked to the 137Gly wild-type polymorphism (Fig.3, lower panel). We hypothesize that the few transcripts with the wild-type 251Trp allele found in the resistant condition originated from 251Trp/251Ser heterozygous individuals (Fig. 3, lower panel). It is to note that we also found a few 251Leu alleles in resistant individuals, in agreement with Bergamo et al. , which showed this mutation in low frequencies in some populations of C. hominivorax. The Trp251Leu mutation was also identified in dimethyl-OP-resistant strains of L. cuprina  and could, therefore, have been present in the ancestor of L. cuprina and C hominivorax. The substitution of Tryptophan at position 251 by Serine or Leucine is specific for resistance to dimethyl-OPs, as demonstrated in M. domestica  and L. cuprina , respectively. Taşkin and colleagues  suggested that the substitution of Trp by smaller residues, such as Serine and Leucine, might enhance the hydrolysis of the dimethyl-OP by α-esterases, resulting in resistance. Although the esterase E3 carried the 251Ser mutation conferring resistance, no change in its expression was detected. However, four other esterases were found up-regulated in the resistant group (an exception to this was one thioesterase found down-regulated). We hypothesize that the up-regulation of these esterases (none of which appear to have polymorphisms) enhances the detoxification of the insecticide.
Other genes related to metabolic and detoxifying processes, including six cytochrome P450s (CYP450s), a family of genes known to be involved in resistance processes, were down-regulated in the larvae that survived the insecticide treatment. This result contrasts with many studies [57,58,59,60], which found that an over-expression of these genes conferred resistance to non-OP insecticides. However, these findings are in agreement with the findings of Carvalho et al. , which reported the under-expression of the Cyp6g1 Cytochrome P450 in OP resistant C. hominivorax individuals. Looking at differential allelic expression, we also found many transcripts involved in metabolic processes displaying condition-specific SNPs (7 CYP450s, 2 glutathione-S-transferases and numerous genes involved in proteolysis). Changes to the metabolism and detoxifying processes, especially the down-regulation of CYP450s can diminish the metabolization of the insecticide into its active form . However, unlike the OP diazinon, which acts through its bioactive oxon metabolite, dichlorvos directly inhibits its target . As such, dichlorvos does not require bioactivation and it is immediately toxic. Another possibility to explain our results is that the larvae display an aversion behavior towards insecticide-poisoned food, which could manifest itself by lower and/or altered metabolization rates. Down-regulation of detoxification genes, such as CYP450s, in the context of insecticide resistance is poorly understood. It was observed in mosquitoes resistant to permethrin  and could indicate an adaptive homeostatic response .
We also observed the up-regulation of heat-shock genes and immune effectors in the resistant condition. Heat-shock proteins (HSPs) allow organisms to adjust their tolerance levels and are produced by cells in response to exposure to environmental stresses. The up-regulation of HSPs are consistently reported in arthropods exposed to insecticides. It was observed in mosquitoes [44, 47, 51], flower thrips , the diamondback moth , and the Asian corn borer . In agreement with previous studies [44, 46, 50], we observed a pattern of up-regulation of genes involved in the immune response after exposure to pesticides.
The altered expression and condition-specific polymorphisms of cuticle genes and genes involved in the chitin metabolic process is suggestive of a reduced penetration of the insecticide in the haemolymph by a thickening or a modification of the permeability of the cuticle, which is a common insecticide resistance mechanism (reviewed by ). Reduction of cuticular insecticide penetration is involved in the resistance of strains of the german cockroach, Blattella germanica [66, 67] and in the malaria vector, Anopheles gambiae . Further experiments are needed to confirm the cuticle changes in C. hominivorax and their consequence on insecticide penetration.
Interestingly, we found an odorant binding protein (OBP) displaying condition-specific polymorphisms and another one was found highly up-regulated in the resistant larvae compared to the control ones. This result was concordant with the identification of condition-specific SNPs in genes putatively involved in olfactory behavior and starvation. In Drosophila melanogaster, OBPs modulate the ingestion of bitter tastants and can contribute directly to taste perception . Lines expressing RNAi corresponding to putative orthologs of the OBP genes identified here, display altered feeding responses. Silencing of these genes in adult flies resulted in a reduction in feeding from sources with specific tastants . Nucleotide changes in OBP sequences can also cause variation in olfactory behavior in response to chemical stimulus in D. melanogaster . Changes in the expression of olfactory genes, such as the odorant-binding protein genes identified in this study, may indicate that resistant larvae are able to perceive the insecticide, which could trigger physiological responses and behavior avoidance of the toxic food. The behavioral avoidance, in turn, could lead to a reduced intake of the insecticide-contaminated food (starvation). Aversion to specific diets has been observed in screwworm larvae: they stop feeding, and this directly impacts the expression of genes involved in metabolic processes (unpublished results). Behavioral avoidance has been reported on multiple occasions and may be a pervasive resistance strategy in many cases [71, 72]. Resistant strains of B. germanica, for instance, avoid harborages treated with cypermethrin, while a susceptible strain apparently is unable to distinguish between treated and untreated harborages . Resistant strains absorb a sublethal amount sufficient for detection, leading to successful avoidance of treated harborages, as they are not overwhelmed by the toxic effects of an insecticide. However, the molecular basis of behavioral resistance remains enigmatic and poorly characterized . Here, we identified clear candidates for behavioral resistance (genes linked to olfactory, starvation and locomotion processes). Avoidance assays and food intake measurements using control and resistant larvae would have to be performed to test the magnitude of the behavioral resistance. If behavior resistance is significant, the specific OBP identified in this study could be targeted using the genetic engineering technique CRISPR/Cas9, which has been successfully carried out in C. hominivorax  or by using RNAi.
In conclusion, it is clear that the Trp251Ser mutation in the esterase E3 is necessary to confer resistance to dimethyl-OP insecticides, however it is uncertain if it is sufficient to provide full dimethyl-OP resistance. The up-regulation of other esterases, as well as the avoidance of insecticide intake, may be important complementary strategies and ought to be considered in our current pest management strategies. Finally, our study demonstrates that the resistance and response mechanisms to dimethyl-OP insecticides involve not only modulations of the expression of genes but also the presence of condition-specific polymorphisms, which is a generally overlooked mechanism in RNA-seq analyses.
Availability of data and materials
The raw sequencing data as well as transcript abundance measurements were deposited in the GEO (Gene Expression Omnibus) repository under the accession number GSE145822.
Clustered Regularly Interspaced Short Palindromic Repeats associated protein-9 nuclease
Organophosphate (OP)-based insecticides
2,2-Dichlorovinyl dimethyl phosphate C4H7Cl2O4P
False discovery rate
Lethal concentration for 90% of the population
Single nucleotide polymorphism
Quantitative polymerase chain reaction
Carvalho RA, Torres TT, de Azeredo-Espin AML. A survey of mutations in the Cochliomyia hominivorax (Diptera: Calliphoridae) esterase E3 gene associated with organophosphate resistance and the molecular identification of mutant alleles. Vet Parasitol. 2006;140:344–51.
Carvalho RA, Limia CEG, Bass C, de Azeredo-Espin AML. Changes in the frequency of the G137D and W251S mutations in the carboxylesterase E3 gene of Cochliomyia hominivorax (Diptera: Calliphoridae) populations from Uruguay. Vet Parasitol. 2010;170:297–301.
Zumpt F. Myiasis in Man and Animals in the Old World: A Textbook for Physicians, Veterinarians, and Zoologists. Butterworths; 1965.
Hall M, Wall R. Myiasis of humans and domestic animals. Adv Parasitol. 1995;35:257–334.
Grisi L, Leite RC, Martins JR, de Barros ATM, Andreotti R, Cançado PHD, et al. Reassessment of the potential economic impact of cattle parasites in Brazil. Rev Bras Parasitol Vet. 2014;23:150–6.
Grisi L, Massard CL, Borja G, Pereira JB. Impacto econômico das principais ectoparasitoses em bovinos no Brasil. Hora Veterinária. 2002;21:8–10.
Sereno JRB, Catto JB, Sereno FTPS. Prevencao de miiases umbilicais em bezerros criados extensivamente, no Pantanal, atraves da utilizacao de ivermectin. 1996.
Batista-da-Silva JA, Moya-Borja GE, Queiroz MMC. Factors of susceptibility of human myiasis caused by the New World screw-worm, Cochliomyia hominivorax in São Gonçalo, Rio de Janeiro. Brazil J Insect Sci Online. 2011;11:14.
Trombetta L, Oliva A, Galache V, Bava J, Troncoso A. Cutaneous myiasis due to Cochliomyia hominivorax in a drug user. J Infect Dev Ctries. 2009;3:873–6.
Solomon M, Lachish T, Schwartz E. Cutaneous Myiasis. Curr Infect Dis Rep. 2016;18:28.
Gomez RS, Perdigão PF, Pimenta FJGS, Rios Leite AC, de Lacerda JC, Custódio Neto AL. Oral myiasis by screwworm Cochliomyia hominivorax. Br J Oral Maxillofac Surg. 2003;41:115–6.
Chodosh J, Clarridge J. Ophthalmomyiasis: a review with special reference to Cochliomyia hominivorax. Clin Infect Dis. 1992;14:444–9.
da Silva BB, Borges US, Pimentel ICC. Human vaginal myiasis caused by Cochliomyia hominivorax. Int J Gynaecol Obstet. 2005;89:152–3.
Osorio J, Moncada L, Molano A, Valderrama S, Gualtero S, Franco-Paredes C. Role of ivermectin in the treatment of severe orbital myiasis due to Cochliomyia hominivorax. Clin Infect Dis. 2006;43:e57-59.
Gealh WC, Ferreira GM, Farah GJ, Teodoro U, Camarini ET. Treatment of oral myiasis caused by Cochliomyia hominivorax: two cases treated with ivermectin. Br J Oral Maxillofac Surg. 2009;47:23–6.
Thyssen PJ, Nassu MP, Costella AMU, Costella ML. Record of oral myiasis by Cochliomyia hominivorax (Diptera: Calliphoridae): case evidencing negligence in the treatment of incapable. Parasitol Res. 2012;111:957–9.
Coronado A, Kowalski A. Current status of the New World screwworm Cochliomyia hominivorax in Venezuela. Med Vet Entomol. 2009;23(Suppl 1):106–10.
El-Azazy OM. Wound myiasis caused by Cochliomyia hominivorax in humans in Libya. Trans R Soc Trop Med Hyg. 1990;84:747–8.
Beesley WN. The New World screw-worm fly in north Africa. Ann Trop Med Parasitol. 1991;85:173–9.
Fukuto TR. Mechanism of action of organophosphorus and carbamate insecticides. Environ Health Perspect. 1990;87:245–54.
Carvalho RA, Torres TT, Paniago MG, Azeredo-Espin AML. Molecular characterization of esterase E3 gene associated with organophosphorus insecticide resistance in the New World screwworm fly, Cochliomyia hominivorax. Med Vet Entomol. 2009;23(Suppl 1):86–91.
Carvalho RA, Azeredo-Espin AML, Torres TT. Deep sequencing of New World screw-worm transcripts to discover genes involved in insecticide resistance. BMC Genomics. 2010;11:695.
Campbell PM, Newcomb RD, Russell RJ, Oakeshott JG. Two different amino acid substitutions in the ali-esterase, E3, confer alternative types of organophosphorus insecticide resistance in the sheep blowfly, Lucilia cuprina. Insect Biochem Mol Biol. 1998;28:139–50.
Newcomb RD, Campbell PM, Ollis DL, Cheah E, Russell RJ, Oakeshott JG. A single amino acid substitution converts a carboxylesterase to an organophosphorus hydrolase and confers insecticide resistance on a blowfly. Proc Natl Acad Sci USA. 1997;94:7464–8.
Heidari R, Devonshire AL, Campbell BE, Dorrian SJ, Oakeshott JG, Russell RJ. Hydrolysis of pyrethroids by carboxylesterases from Lucilia cuprina and Drosophila melanogaster with active sites modified by in vitro mutagenesis. Insect Biochem Mol Biol. 2005;35:597–609.
Vargas ME, Espin AM. Genetic variability in mitochondrial DNA of the screwworm, Cochliomyia hominivorax (Diptera: Calliphoridae), from Brazil. Biochem Genet. 1995;33:237–56.
Andrews S. FastQC A quality control tool for high throughput sequence data. 2010. https://www.bioinformatics.babraham.ac.uk/projects/fastqc/ Accessed 7 Oct 2020.
Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114–20.
Tandonnet S, Torres TT. Traditional versus 3’ RNA-seq in a non-model species. Genomics Data. 2017;11:9–16.
Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29:644–52.
Simão FA, Waterhouse RM, Ioannidis P, Kriventseva EV, Zdobnov EM. BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. Bioinformatics. 2015;31:3210–2.
Chen T-W, Gan R-C, Fang Y-K, Chien K-Y, Liao W-C, Chen C-C, et al. FunctionAnnotator, a versatile and efficient web tool for non-model organism annotation. Sci Rep. 2017;7:10430.
Li W, Godzik A. Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences. Bioinformatics. 2006;22:1658–9.
Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009;10:R25.
Roberts A, Pachter L. Streaming fragment assignment for real-time analysis of sequencing experiments. Nat Methods. 2013;10:71–3.
Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12:323.
Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26:139–40.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.
Sacomoto GAT, Kielbassa J, Chikhi R, Uricaru R, Antoniou P, Sagot M-F, et al. KISSPLICE: de-novo calling alternative splicing events from RNA-seq data. BMC Bioinformatics. 2012;13(Suppl 6):S5.
Kent WJ. BLAT–the BLAST-like alignment tool. Genome Res. 2002;12:656–64.
Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J, et al. De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. Nat Protoc. 2013;8:1494–512.
Götz S, García-Gómez JM, Terol J, Williams TD, Nagaraj SH, Nueda MJ, et al. High-throughput functional annotation and data mining with the Blast2GO suite. Nucleic Acids Res. 2008;36:3420–35.
Marinho MT, Junqueira ACM, Paulo DF, Esposito MC, Villet MH, Azeredo-Espin AML. Molecular phylogenetics of Oestroidea (Diptera: Calyptratae) with emphasis on Calliphoridae: insights into the inter-familial relationships and additional evidence for paraphyly among blowflies. Mol Phylogenet Evol. 2012;65:840–54.
Ingham VA, Wagstaff S, Ranson H. Transcriptomic meta-signatures identified in Anopheles gambiae populations reveal previously undetected insecticide resistance mechanisms. Nat Commun. 2018;9:5282.
Cui L, Rui C, Yang D, Wang Z, Yuan H. De novo transcriptome and expression profile analyses of the Asian corn borer (Ostrinia furnacalis) reveals relevant flubendiamide response genes. BMC Genomics. 2017;18:20.
Ishak IH, Riveron JM, Ibrahim SS, Stott R, Longbottom J, Irving H, et al. The Cytochrome P450 gene CYP6P12 confers pyrethroid resistance in kdr-free Malaysian populations of the dengue vector Aedes albopictus. Sci Rep. 2016;6:24707.
Lv Y, Wang W, Hong S, Lei Z, Fang F, Guo Q, et al. Comparative transcriptome analyses of deltamethrin-susceptible and -resistant Culex pipiens pallens by RNA-seq. Mol Genet Genomics MGG. 2016;291:309–21.
Zhang LJ, Wu ZL, Wang KF, Liu Q, Zhuang HM, Wu G. Trade-off between thermal tolerance and insecticide resistance in Plutella xylostella. Ecol Evol. 2015;5:515–30.
Wang H, Reitz SR, Wang L, Wang S, Li X, Lei Z. The mRNA expression profiles of five heat shock protein genes from Frankliniella occidentalis at different stages and their responses to temperatures and insecticides. J Integr Agric. 2014;13:2196–210.
Lin Q, Jin F, Hu Z, Chen H, Yin F, Li Z, et al. Transcriptome analysis of chlorantraniliprole resistance development in the diamondback moth Plutella xylostella. PLoS ONE. 2013;8:e72314.
Vontas J, Blass C, Koutsos AC, David J-P, Kafatos FC, Louis C, et al. Gene expression in insecticide resistant and susceptible Anopheles gambiae strains constitutively or after insecticide exposure. Insect Mol Biol. 2005;14:509–21.
Hemingway J. The molecular basis of two contrasting metabolic mechanisms of insecticide resistance. Insect Biochem Mol Biol. 2000;30:1009–15.
Montella IR, Schama R, Valle D. The classification of esterases: an important gene family involved in insecticide resistance - a review. Mem Inst Oswaldo Cruz. 2012;107:437–49.
Taşkin V, Kence M, Göçmen B. Determination of malathion and diazinon resistance by sequencing the Md alpha E7 gene from Guatemala, Colombia, Manhattan, and Thailand housefly (Musca domestica L.) strains. Genetika. 2004;40:478–81.
Bergamo LW, Fresia P, Azeredo-Espin AML. Incongruent nuclear and mitochondrial genetic structure of new world screwworm fly populations due to positive selection of mutations associated with dimethyl- and diethyl-organophosphates resistance. PLoS ONE. 2015;10:e0128441.
Smyth KA, Boyce TM, Russell RJ, Oakeshott JG. MCE activities and malathion resistances in field populations of the australian sheep blowfly (Lucilia cuprina). Heredity. 2000;84(Pt 1):63–72.
Le Goff G, Boundy S, Daborn PJ, Yen JL, Sofer L, Lind R, et al. Microarray analysis of cytochrome P450 mediated insecticide resistance in Drosophila. Insect Biochem Mol Biol. 2003;33:701–8.
Ishizuka M, Tanikawa T, Tanaka KD, Heewon M, Okajima F, Sakamoto KQ, et al. Pesticide resistance in wild mammals—mechanisms of anticoagulant resistance in wild rodents. J Toxicol Sci. 2008;33:283–91.
Nikou D, Ranson H, Hemingway J. An adult-specific CYP6 P450 gene is overexpressed in a pyrethroid-resistant strain of the malaria vector, Anopheles gambiae. Gene. 2003;318:91–102.
Bergé JB, Feyereisen R, Amichot M. Cytochrome P450 monooxygenases and insecticide resistance in insects. Philos Trans R Soc Lond B Biol Sci. 1998;353:1701–5.
Feyereisen R. 8 - Insect CYP Genes and P450 Enzymes. In: Gilbert LI, editor. Insect Mol Biol Biochem. San Diego: Academic Press; 2012. p. 236–316.
International Programme on Chemical Safety WHO. Dichlorvos. World Health Organization; 1989.
Yang T, Liu N. Genome analysis of cytochrome P450s and their expression profiles in insecticide resistant mosquitoes, Culex quinquefasciatus. PLoS ONE. 2011;6:e29418.
Morgan ET. Regulation of Cytochromes P450 during inflammation and infection. Drug Metab Rev. 1997;29:1129–88.
Balabanidou V, Grigoraki L, Vontas J. Insect cuticle: a critical determinant of insecticide resistance. Curr Opin Insect Sci. 2018;27:68–74.
Anspaugh DD, Rose RL, Koehler PG, Hodgson E, Roe RM. Multiple mechanisms of pyrethroid resistance in the German cockroach, Blattella germanica (L.). Pestic Biochem Physiol. 1994;50:138–48.
Valles SM, Dong K, Brenner RJ. Mechanisms responsible for cypermethrin resistance in a strain of German cockroach, Blattella germanica. Pestic Biochem Physiol. 2000;66:195–205.
Yahouédo GA, Chandre F, Rossignol M, Ginibre C, Balabanidou V, Mendez NGA, et al. Contributions of cuticle permeability and enzyme detoxification to pyrethroid resistance in the major malaria vector Anopheles gambiae. Sci Rep. 2017;7:11091.
Swarup S, Morozova TV, Sridhar S, Nokes M, Anholt RRH. Modulation of feeding behavior by odorant-binding proteins in Drosophila melanogaster. Chem Senses. 2014;39:125–32.
Wang P, Lyman RF, Mackay TFC, Anholt RRH. Natural variation in odorant recognition among odorant-binding proteins in Drosophila melanogaster. Genetics. 2010;184:759–67.
Sparks TC, Lockwood JA, Byford RL, Graves JB, Leonard BR. The role of behavior in insecticide resistance. Pestic Sci. 1989;26:383–99.
Zalucki MP, Furlong MJ. Behavior as a mechanism of insecticide resistance: evaluation of the evidence. Curr Opin Insect Sci. 2017;21:19–25.
Hostetler ME, Brenner RJ. Behavioral and physiological resistance to insecticides in the German cockroach (Dictyoptera: Blattellidae): an experimental reevaluation. J Econ Entomol. 1994;87:885–93.
Paulo DF, Williamson ME, Arp AP, Li F, Sagel A, Skoda SR, et al. Specific gene disruption in the major livestock pests Cochliomyia hominivorax and Lucilia cuprina using CRISPR/Cas9. G3. 2019;9:3045–55.
The authors would like to thank members of the Torres lab for comments on the manuscript, and two anonymous reviewers for the very constructive comments and criticisms on an earlier version of the manuscript.
This work was supported by grants to TTT from Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP, grants 12/06819-9 and 2014/13933-8). ST was supported by a training grant from Fundação de Amparo a Pesquisa de São Paulo (FAPESP, grant 2013/00243-0). This study was also financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES) - Finance Code 001.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Number of reads in the different samples before and after trimming and collapsing.
Differentially Expressed (DE) transcripts between the control and resistant conditions and their annotations. DE transcripts were those displaying a log2 of the Fold Change (FC) greater than 1 and an adjusted P-value (FDR) lower than 0.05 in all analyses performed (combinations between the counting methods eXpress or RSEM and the testing packages EdgeR and DEseq2). In blue and red are the transcripts found down- and up-regulated in the resistant condition versus the control condition. Negative and positive log2FC indicate whether the transcripts were found down- or up-regulated in the resistant condition versus the control one.
Gene Ontology (GO) terms found over-represented in the transcripts found differentially expressed using the enrichment analysis implemented in BLAST2GO. Sheet 1 and sheet 2 correspond to the complete and reduced sets of GO terms found enriched in the DE transcripts sets, respectively. The reduced list derives from the reduction of the complete list of enriched GO terms to its most specific terms, using BLAST2GO’s functionality. ‘Nr Test’ and ‘Nr Reference’ refer to the number of transcripts annotated with the corresponding term in the test set and the Reference set, respectively. ‘Non Annot’ and ‘Non Annot Reference’ refer to the number of transcripts not annotated with the corresponding GO terms in the test set and reference set, respectively.
Condition-specific polymorphisms between the resistant and control conditions with transcript annotation information. The KisSplice pipeline (see Methods) was used to identify condition-specific polymorphisms (SNPs). We only considered SNPs displaying an adjusted P-value below 0.01, a minimum coverage of 20 per replicate and condition and a magnitude of the differential allelic usage (Absolute Deltaf/DeltaPSI) above 0.2 (20%). Only polymorphisms within the coding region of the transcript and found non-synonymous were further analyzed. Columns 3 to 10 are the normalized read counts of each variant for each replicate and condition. Annotation was performed using FunctionAnnotator (see Methods).
Transcripts found differentially expressed (DE) and displaying condition-specific polymorphisms with their annotations. Annotations were obtained using FunctionAnnotator. For gene expression abundances, we used Express and RSEM. For differential expression testing and differential allelic usage, we used EdgeR, DEseq2 and KisSplice (See Methods).
Gene Ontology (GO) terms found over-represented in the transcripts containing condition-specific SNPs. Sheet 1 and sheet 2 correspond to the complete and reduced sets of GO terms found enriched in the DE transcripts sets, respectively. The reduced list derives from the reduction of the complete list of enriched GO terms to its most specific terms, using BLAST2GO’s functionality. ‘Nr Test’ and ‘Nr Reference’ refer to the number of transcripts annotated with the corresponding term in the test set and the Reference set, respectively. ‘Non Annot’ and ‘Non Annot Reference’ refer to the number of transcripts not annotated with the corresponding GO terms in the test set and reference set, respectively.
About this article
Cite this article
Tandonnet, S., Cardoso, G.A., Mariano-Martins, P. et al. Molecular basis of resistance to organophosphate insecticides in the New World screw-worm fly. Parasites Vectors 13, 562 (2020). https://doi.org/10.1186/s13071-020-04433-3