Skip to main content
Download PDF

Molecular basis of resistance to organophosphate insecticides in the New World screw-worm fly

Parasites & Vectors volume 13, Article number: 562 (2020) Cite this article

Abstract

Background

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.

Methods

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.

Results

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.

Conclusions

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.

Background

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 [3], 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 [4]. 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 [5]. 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 [7].

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 [16]. Infestations in humans are an increasing issue [17] 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 [1]. 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 [20]). The inactivation of AChE leads to an overstimulation and blockage of its receptors (e.g. muscle tissue), resulting in paralysis and insect death [20]. 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 [24] while the replacement of Tryptophan (Trp) by Serine (Ser) at position 251 (Trp251Ser mutation) confers resistance to dimethyl-OPs (e.g. malathion) [23], and possibly to pyrethroid insecticides [25]. As a first attempt to uncover more candidate genes involved in the dimethyl-OP insecticide resistance, Carvalho et al. [22] 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.

Methods

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. [22]. The same resistant and control samples were used in this study and by Carvalho et al. [22] 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. [22]. 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 [26]. 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 [27], 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 [28]. 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 [29].

De novo transcriptome assembly and annotation

The transcriptome assembly was performed using the program Trinity [30] (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 [31] 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 [32] (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 [33] (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 [34] (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 [35] (v1.51) and RSEM [36] (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 [37] (v1.2.25) and DESeq2 [38]. 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).

Variant analysis

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 [39] (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 [40] (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 [41] (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 [42] (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 [42](version 5.2.5).

Results

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 [43].

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.

Fig. 1
figure 1

DE transcripts found by DEseq2 and EdgeR on counts generated by eXpress and RSEM. This diagram shows the overlapping results among four different strategies to identify and validate DE transcripts (control versus resistant). Most of the DE transcripts were found in all comparisons (142 transcripts)

Full size image

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).

Fig. 2.
figure 2

Heatmap of the transcripts (142) found differentially expressed between the control and resistant conditions, organized by functional categories. The expression abundances estimated by RSEM and eXpress were normalized by library size and averaged between both counting methods

Full size image

Metabolic responses

Also consistent with previous studies (see mini-review from [52]), esterase genes, known for their detoxifying role (reviewed in [53]), 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).

Cuticle-related transcripts

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).

Variant analysis

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 [54], the insecticide used in this study.

As expected and in accordance with the Esterase E3 genotyping of Carvalho et al. [23], the resistant individuals had predominantly the esterase E3 (c13624_g1_i1) 251Ser allele (Table 1), which confers resistance to dimethyl-OPs [54]. 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. [1]. It is to note that we detected a few 251Trp alleles in the resistant condition (Table 1).

Table 1 Esterase E3 (transcript c13624_g1_i1) polymorphisms (positions 137 and 251) in the control and resistant conditions. Complete information on these transcripts and SNPs (including amino acid changes) can be found in Additional file 5: Dataset S4
Full size table
Fig. 3.
figure 3

Diagram summarizing the molecular mechanisms associated with the resistance to dimethyl-OP insecticides. Main modulation of gene expression and differential allelic usage in resistant larvae (upper panel). Upward and downward arrows indicate up- and down-regulation in the resistant condition versus the control condition, respectively. Diagram of the Esterase E3 mutations at positions 137 and 251 in the control and resistant conditions (lower panel). Both mutations providing resistance to insecticides, 137Asp and 251Ser, are not found on the same haplotype. In the dimethyl-OP-resistant condition, we encountered almost invariably the 137Gly and 251Ser polymorphisms (see Table 1), which corroborate the findings of Carvalho et al. [25] that showed that 137Gly and 251Ser were on the same haplotype. It is to note that few 251Trp alleles were found in the resistant condition. We hypothesize that they come from heterozygous individuals (dashed line). The skull and bones represent the death of susceptible larvae exposed to the OP insecticide

Full size image

With the present RNA-seq data, it was also possible to explore the expression of genes presenting condition-specific polymorphisms. Using the KisSplice [39] 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.

Table 2 Annotated transcripts found DE and containing condition-specific polymorphisms
Full size table

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.

Discussion

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. [55], 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 [23] 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 [54] and L. cuprina [56], respectively. Taşkin and colleagues [54] 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. [22], 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 [61]. However, unlike the OP diazinon, which acts through its bioactive oxon metabolite, dichlorvos directly inhibits its target [62]. 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 [63] and could indicate an adaptive homeostatic response [64].

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 [49], the diamondback moth [48], and the Asian corn borer [45]. 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 [65]). 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 [68]. 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 [69]. 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 [69]. Nucleotide changes in OBP sequences can also cause variation in olfactory behavior in response to chemical stimulus in D. melanogaster [70]. 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 [73]. 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 [72]. 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 [74] or by using RNAi.

Conclusions

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.

Abbreviations

AChE:

Acetylcholinesterase

Asp:

Aspartic acid

cDNA:

Complementary DNA

CRISPR/Cas9:

Clustered Regularly Interspaced Short Palindromic Repeats associated protein-9 nuclease

CYP450s:

Cytochrome P450s

DDT:

Dichlorodiphenyltrichloroethane

DE:

Differential expression

OP:

Organophosphate (OP)-based insecticides

Dimethyl-OP:

2,2-Dichlorovinyl dimethyl phosphate C4H7Cl2O4P

FDR:

False discovery rate

Gly:

Glycine

GO:

Gene Ontology

LC90:

Lethal concentration for 90% of the population

Leu:

Leucine

SNPs:

Single nucleotide polymorphism

qPCR:

Quantitative polymerase chain reaction

RNA-seq:

RNA-sequencing

Ser:

Serine

Trp:

Tryptophan

References

  1. 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.

    Article  PubMed  CAS  Google Scholar 

  2. 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.

    Article  PubMed  CAS  Google Scholar 

  3. Zumpt F. Myiasis in Man and Animals in the Old World: A Textbook for Physicians, Veterinarians, and Zoologists. Butterworths; 1965.

  4. Hall M, Wall R. Myiasis of humans and domestic animals. Adv Parasitol. 1995;35:257–334.

    Article  CAS  PubMed  Google Scholar 

  5. 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.

    Article  PubMed  Google Scholar 

  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.

    Google Scholar 

  7. Sereno JRB, Catto JB, Sereno FTPS. Prevencao de miiases umbilicais em bezerros criados extensivamente, no Pantanal, atraves da utilizacao de ivermectin. 1996.

  8. 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.

    Google Scholar 

  9. 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.

    Article  PubMed  Google Scholar 

  10. Solomon M, Lachish T, Schwartz E. Cutaneous Myiasis. Curr Infect Dis Rep. 2016;18:28.

    Article  PubMed  Google Scholar 

  11. 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.

    Article  CAS  PubMed  Google Scholar 

  12. Chodosh J, Clarridge J. Ophthalmomyiasis: a review with special reference to Cochliomyia hominivorax. Clin Infect Dis. 1992;14:444–9.

    Article  CAS  PubMed  Google Scholar 

  13. da Silva BB, Borges US, Pimentel ICC. Human vaginal myiasis caused by Cochliomyia hominivorax. Int J Gynaecol Obstet. 2005;89:152–3.

    Article  PubMed  Google Scholar 

  14. 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.

    Article  PubMed  Google Scholar 

  15. 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.

    Article  PubMed  Google Scholar 

  16. 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.

    Article  PubMed  Google Scholar 

  17. Coronado A, Kowalski A. Current status of the New World screwworm Cochliomyia hominivorax in Venezuela. Med Vet Entomol. 2009;23(Suppl 1):106–10.

    Article  PubMed  Google Scholar 

  18. El-Azazy OM. Wound myiasis caused by Cochliomyia hominivorax in humans in Libya. Trans R Soc Trop Med Hyg. 1990;84:747–8.

    Article  CAS  PubMed  Google Scholar 

  19. Beesley WN. The New World screw-worm fly in north Africa. Ann Trop Med Parasitol. 1991;85:173–9.

    Article  CAS  PubMed  Google Scholar 

  20. Fukuto TR. Mechanism of action of organophosphorus and carbamate insecticides. Environ Health Perspect. 1990;87:245–54.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. 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.

    Article  PubMed  Google Scholar 

  22. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. 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.

    Article  CAS  Google Scholar 

  24. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. 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.

    Article  CAS  PubMed  Google Scholar 

  26. Vargas ME, Espin AM. Genetic variability in mitochondrial DNA of the screwworm, Cochliomyia hominivorax (Diptera: Calliphoridae), from Brazil. Biochem Genet. 1995;33:237–56.

    Article  CAS  PubMed  Google Scholar 

  27. 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.

  28. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114–20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Tandonnet S, Torres TT. Traditional versus 3’ RNA-seq in a non-model species. Genomics Data. 2017;11:9–16.

    Article  PubMed  Google Scholar 

  30. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. 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.

    Article  PubMed  CAS  Google Scholar 

  32. 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.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  33. 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.

    Article  CAS  PubMed  Google Scholar 

  34. 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.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  35. Roberts A, Pachter L. Streaming fragment assignment for real-time analysis of sequencing experiments. Nat Methods. 2013;10:71–3.

    Article  CAS  PubMed  Google Scholar 

  36. Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12:323.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26:139–40.

    Article  CAS  PubMed  Google Scholar 

  38. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  39. 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.

    Article  PubMed  PubMed Central  Google Scholar 

  40. Kent WJ. BLAT–the BLAST-like alignment tool. Genome Res. 2002;12:656–64.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. 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.

    Article  CAS  PubMed  Google Scholar 

  42. 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.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  43. 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.

    Article  CAS  PubMed  Google Scholar 

  44. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. 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.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  46. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. 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.

    Article  CAS  PubMed  Google Scholar 

  48. 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.

    Article  PubMed  PubMed Central  Google Scholar 

  49. 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.

    Article  CAS  Google Scholar 

  50. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  51. 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.

    Article  CAS  PubMed  Google Scholar 

  52. Hemingway J. The molecular basis of two contrasting metabolic mechanisms of insecticide resistance. Insect Biochem Mol Biol. 2000;30:1009–15.

    Article  CAS  PubMed  Google Scholar 

  53. 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.

    Article  CAS  PubMed  Google Scholar 

  54. 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.

    PubMed  Google Scholar 

  55. 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.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  56. 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.

    Article  CAS  PubMed  Google Scholar 

  57. 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.

    Article  PubMed  CAS  Google Scholar 

  58. 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.

    Article  CAS  PubMed  Google Scholar 

  59. 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.

    Article  CAS  PubMed  Google Scholar 

  60. 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.

    Article  PubMed  PubMed Central  Google Scholar 

  61. 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.

    Chapter  Google Scholar 

  62. International Programme on Chemical Safety WHO. Dichlorvos. World Health Organization; 1989.

  63. Yang T, Liu N. Genome analysis of cytochrome P450s and their expression profiles in insecticide resistant mosquitoes, Culex quinquefasciatus. PLoS ONE. 2011;6:e29418.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  64. Morgan ET. Regulation of Cytochromes P450 during inflammation and infection. Drug Metab Rev. 1997;29:1129–88.

    Article  CAS  PubMed  Google Scholar 

  65. Balabanidou V, Grigoraki L, Vontas J. Insect cuticle: a critical determinant of insecticide resistance. Curr Opin Insect Sci. 2018;27:68–74.

    Article  PubMed  Google Scholar 

  66. 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.

    Article  Google Scholar 

  67. 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.

    Article  CAS  Google Scholar 

  68. 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.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  69. 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.

    Article  CAS  PubMed  Google Scholar 

  70. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  71. Sparks TC, Lockwood JA, Byford RL, Graves JB, Leonard BR. The role of behavior in insecticide resistance. Pestic Sci. 1989;26:383–99.

    Article  CAS  Google Scholar 

  72. Zalucki MP, Furlong MJ. Behavior as a mechanism of insecticide resistance: evaluation of the evidence. Curr Opin Insect Sci. 2017;21:19–25.

    Article  CAS  PubMed  Google Scholar 

  73. 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.

    Article  CAS  PubMed  Google Scholar 

  74. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

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.

Funding

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.

Author information

Author notes

  1. Renato Assis de Carvalho

    Present address: Bayer Crop Science, Rua Domingos Jorge, 1100, São Paulo, SP, Brasil

Authors and Affiliations

  1. Departamento de Genética e Biologia Evolutiva, Instituto de Biociências, Universidade de São Paulo (USP), São Paulo, SP, Brazil

    Sophie Tandonnet, Gisele Antoniazzi Cardoso, Pedro Mariano-Martins, Raquel Dietsche Monfardini, Vanessa A. S. Cunha, Renato Assis de Carvalho & Tatiana Teixeira Torres

Authors
  1. Sophie Tandonnet
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Gisele Antoniazzi Cardoso
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Pedro Mariano-Martins
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Raquel Dietsche Monfardini
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Vanessa A. S. Cunha
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Renato Assis de Carvalho
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Tatiana Teixeira Torres
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

ST, RAC and TTT carried out larvae experiments and RNA sequencing. ST and TTT carried out all expression and SNP analysis. VASC and ST performed the Esterase E3 analyses. GAC, PMM, RDM and ST annotated the transcriptome. ST and TTT wrote the manuscript. TTT supervised the project and coordinated all activities. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Tatiana Teixeira Torres.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Additional file 1: Table S1.

Number of reads in the different samples before and after trimming and collapsing.

Additional file 2: Dataset S1.

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.

Additional file 3: Dataset S2.

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.

Additional file 4: Dataset S3.

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).

Additional file 5: Dataset S4.

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).

Additional file 6: Dataset S5.

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.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

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://0-doi-org.brum.beds.ac.uk/10.1186/s13071-020-04433-3

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://0-doi-org.brum.beds.ac.uk/10.1186/s13071-020-04433-3

Keywords

Parasites & Vectors

ISSN: 1756-3305

Contact us