Skip to main content

Advertisement

Transcriptome profiling of a spirodiclofen susceptible and resistant strain of the European red mite Panonychus ulmi using strand-specific RNA-seq

Abstract

Background

The European red mite, Panonychus ulmi, is among the most important mite pests in fruit orchards, where it is controlled primarily by acaricide application. However, the species rapidly develops pesticide resistance, and the elucidation of resistance mechanisms for P. ulmi has not kept pace with insects or with the closely related spider mite Tetranychus urticae. The main reason for this lack of knowledge has been the absence of genomic resources needed to investigate the molecular biology of resistance mechanisms.

Results

Here, we provide a comprehensive strand-specific RNA-seq based transcriptome resource for P. ulmi derived from strains susceptible and resistant to the widely used acaricide spirodiclofen. From a de novo assembly of the P. ulmi transcriptome, we manually annotated detoxification enzyme families, target-sites of commonly used acaricides, and horizontally transferred genes implicated in plant-mite interactions and pesticide resistance. In a comparative analysis that incorporated sequences available for Panonychus citri, T. urticae, and insects, we identified radiations for detoxification gene families following the divergence of Panonychus and Tetranychus genera. Finally, we used the replicated RNA-seq data from the spirodiclofen susceptible and resistant strains to describe gene expression changes associated with resistance. A cytochrome P450 monooxygenase, as well as multiple carboxylcholinesterases, were differentially expressed between the susceptible and resistant strains, and provide a molecular entry point for understanding resistance to spirodiclofen, widely used to control P. ulmi populations.

Conclusions

The new genomic resources and data that we present in this study for P. ulmi will substantially facilitate molecular studies of underlying mechanisms involved in acaricide resistance.

Background

The European red mite, Panonychus ulmi, is a member of the family Tetranychidae (Arthropoda: Chelicerata: Arachnida: Acari). Tetranychidae, or spider mites, use their stylet- like chelicerae to puncture the leaf mesophyll-cells to suck out cell contents. This results in chlorotic spots and necrosis and ultimately, leaf abscission [1]. As a result, among mite pests, spider mites are economically important with 80 % of the acaricide market used for their control [2]. Nonetheless, a major problem associated with spider mite control is the rapid development of acaricide resistance. Factors accelerating resistance evolution include not only the frequent use of acaricides, but also high fecundity rates, arrhenotokous reproduction (on haploid males selection is highly effective) and a short life cycle [3]. As a result, the two-spotted spider mite (Tetranychus urticae) and P. ulmi have developed resistance to nearly all acaricide classes [4] and rank among the top 10 of most resistant arthropod species based on the number of unique active ingredients for which resistance has been reported [2].

The acaricides spirodiclofen and spiromesifen belong, together with the insecticide spirotetramat, to the family of spirocyclic tetronic acids (keto-enols) and interrupt lipid biosynthesis by direct inhibition of acetyl coenzyme A carboxylase (ACCase) [5, 6]. Introduced relatively recently to the market [7], spirodiclofen has been shown to be an effective tool in resistance management of tetranychid mites. It shows good to excellent residual control [8, 9] and is effective against spider mites resistant to different acaricide classes [1016]. Although high levels of spirodiclofen resistance have not yet been reported in field populations, artificial selection studies have shown that spider mites have the potential to develop resistance to spirodiclofen comparatively rapid in the laboratory [17, 18]. Selection for resistance was undertaken with T. urticae, P. ulmi and P. citri and resistance in all species was synergized by piperonylbutoxide (PBO), suggesting the involvement of cytochrome P450s [1719]. Strikingly, for all three tetranychid species spirodiclofen resistance was age dependent. The concentration of spirodiclofen causing 50 % lethality (LC50) in the resistant strains was much lower in eggs then in mobile stages, a finding Demaeght et al. [20] suggested could be due to lower expression of detoxification genes in eggs.

In 2011, the high quality genome sequence of the spider mite Tetranychus urticae became available [21], and afforded new resources and tools to investigate the molecular mechanisms underlying acaricide resistance [20, 2225]. For example, although biochemical studies had shown that elevated P450 and esterase levels were associated with resistance to spirodiclofen in T. urticae (see above), studies with a genome-wide gene expression microarray allowed the identification of the genes coding for the detoxification enzymes associated with these activities (CYP392E10 and TuCCE-04) [20]. It is clear that our understanding of resistance to acaricides in phytophagous mites improved with the advent of the T. urticae genome, however, such genomic resources are scarce for other key tetranychid mites, particularly P. ulmi. With only 49 nucleotide sequences available in the NCBI public database (accessed May 2015), P. ulmi resistance research at the molecular level has hardly been possible.

For organisms without a sequenced genome, de novo transcriptome assemblies provide an important starting point for genomic analyses, and recent advances in high-throughput transcript sequencing have greatly contributed to our understanding of gene expression and structure [2629]. This is especially true in agricultural pest management, where the majority of pest species do not have a reference genome and where the elucidation of resistance mechanisms of insecticides is vital [3034].

In this study we used the Illumina HiSeq2000 technology to generate deep paired-end, strand-specific RNA-seq reads for two strains of P. ulmi. The first strain (HS) is susceptible to most currently used acaricides, while the second strain (PSR-TK) was field collected and shown to be resistant to a number of acaricides, including hexythiazox and clofentezine. The strain also showed decreased susceptibility to spirodiclofen, and was further selected in the laboratory until high levels of spirodiclofen resistance (Resistance ratio of 7000) were reached [17]. Paired-end RNAseq reads were used to assemble a P. ulmi transcriptome, and with the same assembly methodology we constructed two P. citri transcriptomes based on previously published RNA-seq data available as single end reads [33, 34]. We then performed a phylogenetic classification of Panonychus gene families involved in detoxification, such as cytochrome P450 monooxygenases (CYPs), carboxylcholinesterases (CCEs), glutathione-S-transferases (GSTs), UDP-glycosyl transferases (UGTs) and ABC transporters (ABCs), and searched for horizontally transferred genes that were previously uncovered in the genome sequence of T. urticae [21, 3539]. Moreover, a differential expression analysis between the acaricide susceptible and resistant P. ulmi strain allowed us to identify genes encoding proteins that are candidates for involvement in acaricide resistance. In addition, we identified and compared P. ulmi expressed sequences for known acaricide target sites. The assembled annotated transcriptome of P. ulmi provides an important resource that should substantially facilitate molecular studies in this species, including the elucidation of detoxification mechanisms of xenobiotics.

Results and discussion

De novo sequence assembly

We performed high-throughput Illumina sequencing of RNA extracted from adult females from acaricide susceptible (HS) and resistant P. ulmi strains (PSR-TK). With the Illumina HiSeq platform we generated 66,308,047 and 95,412,934 strand-specific paired-end reads for the HS and PSR-TK P. ulmi strains, respectively. Two approaches were used to assemble the P. ulmi transcriptome (see Methods for details). The assembly created with CLC Genomics Workbench (CLC) contains 27,777 unique P. ulmi contigs longer than 200 bp, totalling 27.6 Mb of sequence with an overall GC content of 33.6 %. The average contig size was 993 bp and the N50 was 2087 bp. The Velvet/Oases assembly, on the other hand, consists of 30,044 transcripts totalling 18.6 Mb of sequence with a GC content of 34.5 % and an average contig size/N50 of 619 bp/871 bp. Mapping all RNA-seq reads against both transcriptomes (CLC or Velvet/Oases) revealed that the overall mapping success rate, as a measure for assembly quality, was significantly lower for the Velvet/Oases assembly compared to the CLC assembly (see Additional file 1), with an average overall alignment rate of 91.7 % and 73.4 % for the CLC and Velvet/Oases assembly, respectively. Hence, unless otherwise stated, the CLC assembly was used for all further analyses described in this study. All raw reads and the CLC assembly were submitted to NCBI under BioProject PRJNA271858.

As the transcriptome assemblies of two previously published P. citri transcriptomes [33, 34] were not publicly available, and to generate assemblies with similar methodology, we assembled these transcriptomes in a similar way (CLC) as for P. ulmi (Additional files 2 and 3). P. citri assembly statistics were in line with those of P. ulmi, with 25,529 contigs and an average contig size/N50 of 708 bp/1057 bp for the Liu et al. [33] P. citri transcriptome and 32,362 contigs and an average contig size/N50 of 646 bp/969 bp for the Niu et al. [34] P. citri transcriptome.

Homology searches, species distribution and Gene Ontology (GO) analysis

A total of 9190 sequences (33.1 % of all contigs) had a BLAST hit in the NCBI non-redundant (nr) database, with 9109 sequences having a BLASTx hit against the nr protein database and 81 sequences having a BLASTn hit against the nr nucleotide database. As T. urticae gene annotations have not yet been uploaded to the NCBI database, we also performed a local BLASTx search against the T. urticae proteome. We found that 11,250 P. ulmi contigs showed a BLASTx hit with the predicted T. urticae proteome, including 2483 that had no BLASTx-hit in the NCBI nr database. In total, 11,673 P. ulmi contigs (42,0 % of all contigs) had a BLAST-hit in at least one analysis (Additional file 4). Although 42.0 % of contigs with a BLAST-hit may appear to be a relatively small number, we found that this value does not significantly differ from other de novo transcriptome studies with non-model arthropod species (e.g. in [40, 41]). In the Niu et al. transcriptome study [34] of the related P. citri the percentage of contigs was also of similar magnitude (50.7 %). Such a low percentage can be caused by multiple reasons. The majority of de novo transcriptomes contains a considerable amount of short contigs (<240 bp), containing incomplete (<80 AA) ORFs that in many cases will not BLAST at a given E-value threshold. On the other hand some contigs may be non-coding RNAs which do not BLAST with the non-redundant protein/nucleotide database (or just UTR fragments of otherwise coding genes).

Out of the contigs having a BLAST-hit, 94.1 % (10,989 contigs) of the top BLAST hits belonged to metazoan species, 2.9 % to Bacteria (342 contigs), 1.6 % (181 contigs) to Plants, 0.6 % (75 contigs) to Fungi and 0.7 % (86 contigs) to other organisms. The highest homologies for the majority of the metazoan BLAST hits belonged to the Chelicerata (55.7 %, 6498 contigs) and Arthropoda (20.6 %, 2404 contigs) (Additional file 5). A detailed analysis of all chelicerate BLAST hits revealed a high number of top BLAST hits with T. urticae (21.9 %, 2553 contigs), followed by the social velvet spider Stegodyphus mimosarum (17.8 %, 2080 contigs), the tick Ixodes scapularis (8.5 %, 987 contigs) and the predatory mite Metaseiulus occidentalis (4.7 %, 545 contigs).

We performed a Gene Ontology (GO) analysis of the P. ulmi assembled sequences after excluding contigs (227) that were considered as foreign contamination and consisted mainly of sequences belonging to Wolbachia sp. (111 contigs) and the host plant (Prunus sp.) of the P. ulmi HS and PSR-TK strains (79 contigs) (see Methods). GO terms were assigned to 5833 unique contigs (21.2 % of 27,550 contigs) and in most cases multiple GO terms were assigned to the same P. ulmi contig. We found that 17,161 (54.8 %), 6623 (21.2 %) and 7534 (24.1 %) GO terms emerged for biological process, cellular component and molecular function categories, respectively. Those were further categorized into 16 biological process, 9 molecular function and 7 cellular component sub-categories (Fig. 1). Among the 16 sub-categories of biological process, the cellular process occupied the highest number (3535, 20.6 %), followed by the metabolic process (2858, 16.7 %), single organism process (2853, 16.6 %) and response to stimulus (1308, 7.6 %) (Fig. 1). The major sub-categories of the molecular function category included binding (2705, 40.8 %), followed by catalytic activity (2409, 36.4 %) and transporter activity (445, 6.7 %) (Fig. 1). In the cellular component domain, the majority of the GO terms were shown to be specific for the cell (3047, 40.4 %), followed by organelle (2101, 27.9 %) and macromolecular complex (1355, 18,0 %) (Fig. 1). To further increase the annotation functionality, all contigs were mapped to the reference canonical pathways in the Kyoto Encyclopedia of Genes and Genomes (KEGG). KEGG mapping resulted in 861 contigs that were assigned to 117 KEGG pathways, with 35.4 % of P. ulmi KEGG enzymes identified as transferases, 25.2 % as hydrolases, 22.4 % as oxidoreductases, 8.1 % as ligases, 6.2 % as lyases, while isomerases constituted for 2.7 % of annotated enzymes (Additional file 6). Our GO and KEGG analyses of the P. ulmi transcriptome are in line with previous GO analyses of P. citri transcriptomes where cellular process, binding and cell were the three major sub-categories [33, 34] and indicates that our approach using the Illumina HiSeq platform provided an extensive representation of the P. ulmi transcriptome.

Fig. 1
figure1

Gene ontology annotation and classification (level 2) of the P. ulmi transcriptome. Results are summarized into three main categories: molecular function, biological process and cellular component. In total GO terms were assigned to 5833 P. ulmi contigs

Analysis of Panonychus genes involved in xenobiotic resistance

In general, arthropods have developed two types of mechanisms to cope with xenobiotic compounds, both of which can contribute to the development of resistance: mechanisms that decrease exposure due to quantitative or qualitative changes in major detoxification enzymes and transporters (pharmacokinetic mechanisms) and mechanisms that decrease sensitivity due to changes in target site sensitivity caused by point mutations (pharmacodynamic mechanisms) [2, 42, 43]. The pharmacokinetic mechanism can be subdivided into three phases (I-III). In phase I detoxification enzymes such as cytochrome P450 monooxygenases (CYPs) and carboxylcholinesterases (CCEs) incorporate a nucleophilic functional group (a hydroxyl, carboxyl or amine group) into the toxic compound, resulting in a more reactive and water soluble compound. During phase II, enzymes such as UDP-glycosyltransferases (UGTs) and glutathione-S-transferases (GSTs) further increase the water solubility of the phase I metabolite by conjugation with endogenous molecules like sugars and glutathione, respectively. In Phase III, conjugates are transported out of the cell by cellular transporters, e.g., ABC transporters (ABCs).

In this study we mined the P. ulmi and P. citri transcriptomes for genes encoding known target-sites and major detoxification enzymes and transporters (CYPs, CCEs, UGTs, GSTs and ABCs) and compared them with those of T. urticae, for which a high-quality genomic assembly and annotation is available [21]. In addition, we also performed a phylogenetic analysis for all major detoxification gene families. It should be noted, however, that as transcriptomic data for Panonychus sp. is compared to genomic data for T. urticae, care should be taken when comparing numbers and ortho/homologues of Panonychus detoxification genes as recent gene duplications and genes with very low expression levels might be missed. In addition, expression of detoxification genes might also be stage-dependent. Hence, gene number differences between both Panonychus species should also be carefully interpreted as both P. citri transcriptomes were assembled using RNA-seq data from mixed stages [33, 34] while the P. ulmi transcriptome was based solely on RNA-seq data of adult females (see Methods).

Cytochrome P450 monooxygenases

Cytochrome P450 monooxygenases (P450s) are heme-containing enzymes with a diverse range of functions. Many P450s are important phase I detoxification enzymes with a crucial role in detoxification of plant secondary metabolites and in metabolizing insecticides/acaricides to less toxic compounds. Four major clans can be distinguished in the CYP gene family, namely the CYP2, CYP3, CYP4 and M (mitochondrial CYP genes) [44]. The advent of the T. urticae genome allowed a first insight into the CYP gene family of a phytophagous mite, revealing 81 full length CYPs. This number is similar to what is found in insects, but with an expansion of T. urticae specific intronless genes of the CYP2 clan [21]. A total of 63 P. ulmi and 118 P. citri CYP non-allelic ORFs could be identified (Table 1, Additional file 7). Only those CYPs that did not misalign in the final alignment and had a minimum ORF length of 450 nt were included in a phylogenetic analysis (41 P. ulmi and 49 P. Citri ORFs) (Fig. 2). Based on this analysis or based on its best BLASTx hit with a T. urticae CYP, all P. ulmi and P. citri CYPs could be assigned to one of the 4 CYP clans (Fig. 2, Table 1, Additional file 7) [21].

Table 1 Comparison of the CYP gene number in T. urticae, P. ulmi and P. citri
Fig. 2
figure2

Phylogenetic analysis of P. ulmi putative CYPs. Maximum likelihood tree of P. ulmi, P. citri and T. urticae CYP protein sequences. The tree was constructed using MUSCLE [87] and Treefinder [89]. The tree was midpoint rooted and numbers at the branch point of each node represent the bootstrap value resulting from 1000 pseudoreplicates (LR-ELW). The scale bar represents 0.5 substitutions per site. Tetranychid CYPs clustered into the four main CYP clades: mitochondrial CYPs, CYP 2, CYP 3 and CYP 4. Colour and shape codes are as follows: P. ulmi, green triangle, P. citri, yellow triangle and T. urticae, red square. An arrow and an asterisk indicate CYPs associated with high levels of spirodiclofen resistance in P. ulmi (this study) and T. urticae [20], respectively. Tetranychid CYP protein sequences and accession numbers can be found in Additional file 7

In vertebrates, mitochondrial CYPs (clan M) are associated with steroid or vitamin D metabolism while in insects this clan comprises two groups: one group involved in the ecdysteroid metabolism pathway and another group that metabolize a variety of xenobiotic compounds [44]. We identified 9 P. ulmi and 11 P. citri ORFs in this clan and detected three clear cases of 1:1:1 orthology between P. ulmi, P. citri and T. urticae: CYP314A1, CYP315A1 and CYP302A1 (Additional file 7, Additional file 8, Table 1). Together with T. urticae CYP307A1 of the CYP2 clan (see below), these CYPs are the only T. urticae CYPs that could be assigned as clear orthologs of insect and crustacean CYPs. They are known to be involved in the ecdysteroid pathway in insects and their tetranychid counterparts likely fulfill an analogous role [21]. Similar to T. urticae, orthologs of insect CYP2 clan members CYP306A1 and CYP18A1 could also not be identified for both Panonychus species (Additional file 7, Fig. 2). In insects, CYP306A1s are involved in the synthesis of 20 hydroxyecdysone (20E) by hydroxylating carbon 25, while CYP18A1 is a hydroxylase/oxidase involved in ecdysteroid degradation [21]. Their absence in Panonychus sp. further corroborates that this gene cluster was lost as a whole and is an ancient trait, leading to the use of ponasterone A as the moulting hormone instead of 20E.

In insects, the number of CYP2 genes is relatively well conserved and several insect CYP2 genes are known to be involved in essential physiological functions. Compared to insects the number of CYP2 genes is more diverse in T. urticae and the cladoceran Daphnia sp. [21, 44] and at present the 392 family within the CYP2 clan is best characterized in T. urticae. A whole genome microarray gene expression analysis revealed that several members of the CYP2 clan were upregulated in acaricide resistant T. urticae strains or in a T. urticae strain adapted to a new challenging host plant [20, 36]. Furthermore, functional expression of two T. urticae CYP2 clan enzymes, CYP392A16 and CYP392E10, confirmed that these could metabolise the acaricides spirodiclofen and abamectin, respectively [20, 45]. We identified 16 and 31 CYP2 ORFs in P. ulmi and P. citri respectively. In our phylogenetic analysis, there are two clear cases of 1:1:1 orthology within the CYP2 clan (T. urticae CYP307A1, and CYP392A5). The phylogenetic analysis also suggests that the potential expansion/loss events of certain CYP2 subfamilies, like CYP392A, CYP392D and CYP392E, occurred after the split between the Tetranychus and Panonychus genus. Consequently, it was not surprising that we could not identify a clear P. ulmi/P. citri ortholog of T. urticae CYP392E10 and CYP392A16, suggesting that spirodiclofen or abamectin metabolism must depend on different P450s or other factors in Panonychus species (see differential expression analysis).

Insect CYP3 and CYP4 clans comprise the insect specific families CYP4, CYP6, CYP9 and CYP325 that are well known for their involvement in environmental response/detoxification functions against xenobiotics and phytotoxins [44]. We identified 8 P. ulmi and 26 P. citri ORFs in the CYP3 clan, while 30 P. ulmi and 50 P. citri ORFs were found in the CYP4 clan. Within the CYP3 and CYP4 clan we found 5 (T. urticae CYP382A1, CYP384A1, CYP385A1, CYP385B1, CYP385C1) and 3 (T. urticae CYP407A1, CYP391A1, and CYP4CL1) clear cases of 1:1:1 orthology between T. urticae, P. ulmi and P. citri, respectively. Within the CYP4 clan, the CYP389A subfamily consists of a single gene in T. urticae while 4 paralogs are present in each, P. ulmi and P. citri. In contrast, the CYP389C subfamily appears to be larger in T. urticae compared to the two Panonychus species. Interestingly, Ran et al. [46] showed that P. citri paralogs of T. urticae CYP389A1 were among the most up- and down regulated CYP genes in an amitraz resistant strain. Ding et al. [47], on the other hand, showed that a P. citri ortholog of T. urticae CYP4CF2 was highly induced upon exposure to pyridaben while the expression of a P. citri ortholog of T. urticae CYP4CL1 increased after induction by abamectin, azocyclotin, pyridaben, and spirodiclofen. In contrast to the CYP2 clan, none of these tetranychid CYP3 and CYP4 clan members have been functionally characterized.

Carboxylcholinesterases

Carboxylcholinesterases (CCEs) catalyse the hydrolysis of carboxylesters and have a whole gamut of physiological functions including degradation of neurotransmitters, metabolizing hormones and pheromones, regulation of behaviour and detoxification of xenobiotics [48, 49]. CCEs are phase I detoxification enzymes and their role in xenobiotic metabolism and the development of resistance in arthropods has been thoroughly studied (recently reviewed in [42]). A classification based on the CCE phylogeny, as proposed by [48] groups insect CCEs into 13 clades spread over 3 classes: the dietary/detoxification enzymes (clades A–C), the generally secreted enzymes (hormone/semiochemical processing class, clades D–G) and the neuro/developmental CCEs (clades I–M, mainly non-catalytic esterases) [48, 49]. While total numbers of CCEs in insects and T. urticae are similar, the neuro/developmental class is expanded in T. urticae with two new subclades J’ and J”. In addition, in contrast to insects no T. urticae CCEs have been identified within the currently recognized classification of insect dietary CCEs [21].

We found a total of 61 and 91 CCE non-allelic ORFs in the P. ulmi and P. citri transcriptome, of which 41 and 40 ORFs of P. ulmi and P. citri, respectively, were included into a phylogenetic analysis. Based on this analysis or the best BLASTx hit with T. urticae CCEs, P. ulmi and P. citri CCEs could be assigned to one of the 13 CCE clades (Fig. 3, Table 2, Additional file 9). Similar to T. urticae CCEs, the majority of Panonychus CCEs belong to the neuro/developmental clade, with 18 P. ulmi and 25 P. citri CCEs in clade J” and 17 P. ulmi and 23 P. citri CCEs in clade J’ (Table 2, Fig. 3, Additional file 9). In sixteen cases, 1:1:1 orthologous groups could be identified between T. urticae and Panonychus species: 1 in clade J (AChE), 1 in clade K (gliotactin), 1 in clade L (neuroligins), 1 in clade J’, 2 in clade F’ (Acari/Crustacean Juvenile hormone CCEs), 3 in undetermined clades (NC) and 7 in clade J”. Remarkably, within the subclade J’, a putative expansion of CCEs in T. urticae (17 CCEs) was found, suggesting these T. urticae J’ CCEs arose after the Tetranychus/Panonychus divergence within the Tetranychidae. Finally, in clade J’ two clear P. ulmi orthologues of P. citri PcCCE16 and PcCCE39 could be identified. The latter have their best BLASTn (E-value of 0.0) hit with P. citri PCE1 and PCE2, respectively, which have been shown to be induced by acaricide exposure [50].

Fig. 3
figure3

Phylogenetic analysis of P. ulmi putative CCEs. A set of P. ulmi, P. citri, T. urticae, D. melanogaster, Anopheles gambiae and Apis mellifera CCE protein sequences were aligned using MUSCLE, subsequently trimmed according to [49] and subjected to a maximum-likelihood analysis using Treefinder [89]. The tree was midpoint rooted and numbers at the branch point of each node represent the bootstrap value resulting from 1000 pseudoreplicates (LR-ELW). The scale bar represents 0.5 substitutions per site. Tetranychid CCEs clustered into clades [21,49]: F’ (Crustacean/Acari JhE), I, J (AChEs), J, J”, K (gliotactins), L (neuroligins), M (neurotactins) and three undetermined clades (NC). Colour and shape codes are as follows: P. ulmi, green triangle, P. citri, yellow triangle, T. urticae, red square, D. melanogaster, black dot, Anopheles gambiae, purple dot and Apis mellifera, blue dot. An arrow and an asterisk indicate CCEs associated with high levels of spirodiclofen resistance in P. ulmi (this study) and T. urticae [20], respectively. CCE protein sequences and accession numbers can be found in Additional file 9

Table 2 Comparison of the CCE gene number in T. urticae, P. ulmi, P. citri and D. melanogaster

Glutathione-S-transferases

Glutathione-S-transferases (GSTs) belong to a multifunctional enzyme family, known to play an important role in phase II of xenobiotic detoxification. In addition to pesticide detoxification by conjugation, insect GSTs have also been reported to play a role in attenuation of oxidative stress caused by pesticide exposure and sporadically by metabolism of pesticides directly [51]. Insects cytosolic GSTs can be divided into seven classes [delta (δ), epsilon (ε), omega (ω), sigma (σ), theta (θ) and zeta (ζ)] with delta and epsilon GSTs well known for their role in detoxification of organophosphates and organochlorines [42, 52]. In contrast to insects, Acari also have mu-class GSTs, which were until recently thought to be verterbrate specific [21, 53].

We found a total of 19 P. ulmi and 37 P. citri GST non-allelic ORFs of which 13 P. ulmi and 23 P. citri GSTs were included in a phylogenetic analysis (Table 3, Fig. 4 Additional file 10). Based on this analysis or based on the best BLASTp hit with T. urticae GSTs, P. ulmi and P. citri GSTs could be assigned to 4 GST-classes: delta, mu, omega and zeta. The number of P. ulmi GSTs is lower than those observed in P. citri, T. urticae and Ixodes scapularis (Table 3). Although BLASTp searches against the non-redundant database performed by [34] identified representatives of sigma and theta GST we were unable to find neither of those two classes in the transcriptome of P. ulmi. The lack of theta and sigma GSTs seems to be a common feature of Acari species analysed to date [53, 54]. Similar to the phylogenetic analysis of tetranychid CYP genes (see above) only limited orthologous relationships between tetranychid GSTs could be established, as illustrated by the presence of T. urticae specific clusters of delta and mu GSTs. In both the omega and zeta class we identified a 1:1:1 orthology between Panonychus contigs and T. urticae GSTs (TuGSTo01: PcGST6: PuGST8 and TuGSTz01: PcGST14: PuGST7). Within the mu GSTs we found two 1:1:1 orthologs (TuGSTm11/TuGSTm05: PuGST12: PcGST9 and TuGSTm02: PcGST24: PuGST11) while no 1:1:1 orthologs could be identified in the delta GSTs (Fig. 4).

Table 3 Comparison of the cytosolic GST gene number in T. urticae, P. ulmi, P. citri and I. scapularis
Fig. 4
figure4

Phylogenetic analysis of P. ulmi putative GSTs. A set of P. ulmi, P. citri, T. urticae, Ixodes scapularis, Sarcoptes scabiei, Dermatophagoides pteronyssinus GST protein sequences were aligned using MUSCLE and subjected to a maximum-likelihood analysis using Treefinder. The tree was midpoint rooted and numbers at the branch point of each node represent the bootstrap value resulting from 1000 pseudoreplicates (LR-ELW). Tetranychid GSTs clustered within classes: δ, delta class, ζ, zeta class, ω, omega class and μ, mu class. Colour and shape codes are as follows: P. ulmi, green triangle, P. citri, yellow triangle, T. urticae, red square, purple dot, I. scapularis, blue rhombus, S. scabiei and pink rhombus, Dermatophagoides pteronyssinus. GST protein sequences and accession numbers can be found in Additional file 10

Using a whole genome gene-expression microarray, members of the same GST classes were also shown to be highly upregulated in multi-resistant T. urticae strains [36]. Three GST genes (TuGSTd10, TuGSTd14, TuGSTm09) that were upregulated in a strain MAR-AB highly resistant against abamectin and other important acaricides were also functionally expressed. Out of these, TuGSTd14 showed the highest affinity toward abamectin and had a competitive type of inhibition, suggesting the acaricide may bind to the H-site of the enzyme [55]. However, as mentioned earlier no clear Panonychus orthologs of TuGSTd14 could be identified in our study.

ABC transporters

The ATP Binding Cassettte (ABC) protein family is a large and ubiquitous family of proteins. Most members of this family use ATP to transport substrates across lipid membranes, and hence are referred to as ABC transporters. Based on the sequence similarity of their nucleotide binding domain (NBD, the domain that binds ATP) the ABC protein family can be divided into eight subfamilies (A to H). Although extensively characterized for their involvement in drug resistance in vertebrates, the role of ABC transporters in arthropod xenobiotic resistance is less well known. Nevertheless, ABC transporters have been associated with resistance to as many as 27 insecticides/acaricides from nine different chemical classes [56]. However, the major source of evidence linking these arthropod ABC transporters to resistance has been ABC transporter expression quantification or synergism studies [56]. A total of 103 ABC genes were identified in the T. urticae genome, the highest number discovered in a metazoan species to date. This high number of ABC genes in T. urticae is due predominantly to lineage-specific expansions of the ABCC, ABCG and ABCH subfamilies [57].

We found 145 P. ulmi and 230 P. citri ORFs with a tBLASTn hit with T. urticae ABC protein sequences (Additional file 11). NBDs from the ORFs of these contigs, when present, were extracted and used in a maximum likelihood phylogenetic analysis. Similar to previous studies [57] C-terminal NBDs of ABCC transporters clustered together with NBDs of ABCB transporters (Fig. 5, Additional file 11). Based on this analysis or, in case the ORF did not encode an NBD, based on its best BLASTx hit with T. urticae ABC proteins, we assigned P. ulmi and P. citri ABC ORFs to one of the 8 ABC gene subfamilies (Table 4, Fig. 5, Additional file 11). As ABC genes are large genes (around 2 to 5 kb in size), they can be easily fragmented into multiple contigs assembled from RNA-seq reads. This could result in an overestimation of the number of putative ABC genes. Therefore, we used the number of NBDs to estimate quantitative differences among the ABC gene families in the Tetranychidae (Table 4). The number of T. urticae NBDs is much higher than the number of P. citri and P. ulmi NBDs, and is mainly due to a higher number of T. urticae NBDs in the ABCA, ABCC, ABCG and ABCH gene subfamily. In several cases T. urticae ABCH NBDs do not have a counterpart in Panonychus species while in both the ABCA, ABCC as ABCG gene family a putative expansion of T. urticae NBDs compared to Panonychus NBDs can be suspected. Interestingly, for exactly two of these T. urticae ABC gene subfamilies it was shown that several genes were differentially expressed in multi-pesticide resistant strains and/or in mites transferred to challenging host plants [36]. Finally, NBDs of ABC transporters that are considered as conserved in arthropods (ABCB half transporters (D. melanogaster CCG7955, CG1824), ABCC (D. melanogaster CG7806, sur), ABCD, ABCE, ABCF, ABCG (D. melanogaster CG3327, CG11069, CG31121)) were also found in both Panonychus species (Table 4, Fig. 5) [56].

Fig. 5
figure5

Phylogenetic analysis of P. ulmi ABC protein NBDs. Maximum likelihood midpoint rooted tree of N- and/or C-terminal NBDs extracted from P. ulmi and P. citri ABC protein coding ORFs and T. urticae and D. melanogaster ABC protein sequences. The scale bar represents 0.5 amino-acid substitutions per site. P. ulmi and P. citri NBDs clustered within the eight currently described ABC subfamilies (a to h). The tree was midpoint rooted and numbers at the branch point of each node represent the bootstrap value resulting from 1000 pseudoreplicates (LR-ELW). C-terminal NBDs of ABCC transporters clustered together with NBDs of ABCB transporters [57]. Colour and shape codes are as follows: P. ulmi, green triangle, P. citri, yellow triangle, T. urticae, red square and D. melanogaster, black dot. ABC protein sequences and accession numbers can be found in Additional file 11

Table 4 Comparison of the number of P. ulmi, P. citri and T. urticae ABC protein NBDs

Horizontally transferred genes

Horizontal gene transfer (HGT) refers to the asexual transfer of genetic information between non-related species. In bacteria and fungi HGT is a strong evolutionary force contributing to adaptation and colonization of challenging environments. Due to the unique characteristics of animal biology, the prevalence and impact of HGT on animal evolution was considered to be minor. However, recent studies have uncovered a vast array of horizontally transferred genes (HTGs) in the T. urticae genome with plausible roles in xenobiotics detoxification or other responses to the environment. These include 80 UDP-glycosyltransferases [39], 16 intradiol ring-cleavage dioxygenases [36], 3 carotenoid desaturases, 2 carotenoid cyclase/synthases [21, 35], 2 levanases [21], a cyanate lyase [37], β-cyanoalanine synthase [38] and a methionine synthase [21]. Orthologs of these T. urticae HTGs were identified in both Panonychus transcriptomes (Table 5 and 6, Additional file 12). The presence of these genes in the Panonychus genus strongly suggests that these horizontal gene transfers occurred before the split of the Tetranychus and Panonychus genus. As the majority of the microbial xenologues of these unique spider mite genes code for enzymes able to detoxify and break down plant secondary metabolites, it has been argued that HGT facilitated spider mite adaptations to a phytophagous lifestyle [38]. In the sections below we will discuss some of these tetranychid HGT families more into detail.

Table 5 Comparison of the UGT gene number in T. urticae, P. ulmi, P. citri and I. scapularis
Table 6 Putative horizontally transferred genes identified in the T. urticae genome and Panonychus spp. transcriptomes

UDP-glycosyltransferases

UDP-glycosyltransferases (UGTs) are common in the majority of living organisms including viruses, bacteria, plants and animals. They catalyze the conjugation of small lipophilic molecules with uridine diphosphate (UDP) sugars, increasing their water solubility. As such, UGTs play an important role in the synthesis, storage and transport of secondary metabolites. In vertebrates, UGTs are also well studied because of their role in phase II drug metabolism [39, 58]. Although a role for UGTs in arthropod xenobiotic resistance was suggested more than twenty years ago [59], only recently has functional evidence for their role in xenobiotic resistance been presented [6063]. It has been suggested that the UGT gene family might have been lost early in the Chelicerata lineage and subsequently re-gained in the tetranychid mites by means of the horizontal gene transfer from bacteria. This discovery provides important clues to UGTs functions in relation to detoxification and therefore host adaptation in the phytophagous mites.

We identified a total of 33 and 52 P. ulmi and P. citri UGT non-allelic ORFs, respectively. A subset of 24 P. ulmi and 32 P. citri ORFs were included in a phylogenetic analysis (Fig. 6). Based on this analysis or based on the T. urticae best BLASTx hit, P. ulmi and P. citri UGTs could be assigned to one of the UGT subfamilies (Fig. 6, Table 5, Additional file 12). Within Tetranychidae the UGT201 subfamily is the largest UGT subfamily with 10, 23 and 36 UGT ORFs/genes in P. ulmi, P. citri and T. urticae, respectively (Table 5, Additional file 12). The UGT201A, UGT201B, UGT202A and UGT204A subfamily are clearly more numerous in T. urticae compared to both Panonychus species while the opposite could be observed for the UGT201G subfamily, suggesting the UGTs in these subfamilies probably arose in T. urticae or were lost in P. ulmi after diversification within the Tetranychidae. Finally, within different UGT subfamilies clear 1:1:1 orthologous relationships could be observed between Panonychus UGTs and T. urticae UGT202B1, UGT203G1, UGT203F1, UGT204C1, UGT205A3, UGT206A1 and UGT207A1 (Fig. 6).

Fig. 6
figure6

Phylogenetic analysis of P. ulmi putative UGTs. Maximum likelihood midpoint rooted tree of P. ulmi, P. citri and T. urticae UGT protein sequences. The scale bar represents 0.5 amino-acid substitutions per site. Numbers at the branch point of each node represent the bootstrap value resulting from 1000 pseudoreplicates (LR-ELW). Color and shape codes are as follows: P. ulmi, green triangle, P. citri, yellow triangle and T. urticae, red square. UGT protein sequences can be found in Additional file 12

Intradiol ring-cleavage dioxygenases (IDRCDs)

Among Metazoa, tetranychid mites are the only species known to harbor genes encoding IDRCDs. In bacteria and fungi these enzymes catalyze the oxygenolytic fission of catecholic substances, allowing them to degrade aromatic rings, an essential step in the carbon cycle [21, 36]. However, their role in tetranychid mites has not yet been characterized. About half of the number of genes in this family were differentially expressed in mites upon host plant change and in multi-resistant T. urticae strains, and their expression patterns were highly correlated, suggestive of a role for this gene family in xenobiotic resistance [36]. Twelve P. ulmi and 23 P. citri contigs showed a tBLASTn hit with T. urticae IDRCD proteins. A phylogenetic analysis revealed several orthologous relationships between tetranychid IDRCDs (tetur01g00490, tetur04g00150, tetur04g08620, tetur07g02040, tetur10g00490, tetur19g02300, tetur20g01160), suggesting these IDRCD genes arose before the split between the Tetranychus and Panonychus genus (Fig. 7). Furthermore, an expansion of IDRCDs can be observed in both T. urticae (17 IDRCD genes) as in P. ulmi and P. citri (12 and 23 genes respectively). Remarkably, exactly those IDRCDs that are absent in both Panonychus species are highly upregulated in mites upon host plant change and in multi-resistant T. urticae strains [36].

Fig. 7
figure7

Phylogenetic analysis of P. ulmi putative IDRCDs. Maximum likelihood midpoint rooted tree of P. ulmi, P. citri, T. urticae and T. evansi IDRCD protein sequences. The scale bar represents 0.2 amino-acid substitutions per site. Numbers at the branch point of each node represent the bootstrap value resulting from 1000 pseudoreplicates (LR-ELW). Panonychus sp. and Tetranychus sp. specific clades are shaded in blue and green, respectively. IDRCD protein sequences and accession numbers can be found in Additional file 12

Other horizontal gene transfers

Besides UGTs and IDRCDs we also identified Panonychus orthologs of T. urticae carotenoid desaturases, carotenoid cyclases, levanases, cyanate lyase, β-cyanoalanine synthase and methionine synthase (Table 6, Additional file 12). Wybouw et al. [38] showed that a T. urticae strain overexpresses β-cyanoalanine synthase when adapted to cyanogenic plants. Functional expression also revealed that this enzyme is able to detoxify cyanide which is the main defensive phytochemical of these plants. In contrast to T. urticae, Panonychus sp. almost exclusively feed on Rosaceae [64], a cyanogenic plant family producing cyanide in various tissues as plant defense [65]. As such, β-cyanoalanine synthase might have been crucial in conferring resistance in Panonychus sp. to this continuous exposure to dietary cyanide.

Target-sites of insecticides/acaricides

Since the PSR-TK strain originates from the field and has evolved resistance to at least hexythiazox, clofentezine and spirodiclofen [17], the P. ulmi transcriptome was mined for contigs encoding known target sites of acaricides, in order to search for SNPs associated with target-site resistance. We successfully annotated acetylcholinesterase (AChE), voltage gated sodium channel (VGSC), GABA- and glutamate-gated chloride channels (Rdl and GluCl), chitin synthase 1 (CHS1), acetyl-CoA carboxylase (ACCase) and cytochrome B (cytB) [6, 23, 24, 42, 6672] (see Table 7 for an overview and the Insecticide Resistance Action Committee (IRAC) mode of action classes that target these proteins [41]). To obtain P. ulmi target-site sequences as full-length as possible, we also mined an alternative P. ulmi assembly, which was constructed using the Velvet/Oases package (Additional file 13), and combined contigs from the latter assembly with the one described in this study (CLC Genomics Workbench).

Table 7 Known target sites of acaricides, their T. urticae gene ID and presence in P. ulmi and the IRAC acaricide classes that target these proteins

Compared to insects T. urticae has a higher number of Rdl and GluCl genes, with the majority of insects having only one Rdl and GluCl gene while T. urticae has 3 and 6, respectively [23]. Similar to T. urticae we identified 5 P. ulmi GluCl and 3 Rdl genes (Additional file 14), suggesting that the Rdl and GluCl gene family diversified before the radiation of the Tetranychidae. Next, we compared all known target-site sequences of the spirodiclofen susceptible (HS) and resistant (PSR-TK) P. ulmi strains to identify non-synonymous single nucleotide polymorphisms (SNPs) that are unique for the resistant strain. For AChE, we identified the F331W substitution (Torpedo californica numbering) that was fixed in the PSR-TK strain but segregating in the HS strain. In the past, it has been shown that a substitution at this position causes resistance to organophosphate (OP) and carbamate (CB) insecticides/acaricides [42] and the presence of this substitution in both P. ulmi strains emphasizes the scope of selection exerted by OPs and CBs during the second half of the 20th century.

Comparing ACCase sequences of both strains revealed several non-synonymous fixed SNPs in the PSR-TK strain that were not present in the susceptible HS-strain. At present, only one ACCase substitution has been associated with resistance against cyclic keto-enols: an E645K substitution in the whitefly Trialeurodes vaporariorum [73]. However, the residue is not conserved in mites and at the corresponding position in the P. ulmi ACCase protein, the sequence is identical in both strains (a glutamine; at position 566, T. urticae numbering tetur21g02170). Recently, it has been shown that the activated enol derivative of the insecticide spirotetramate interacts with the carboxyltransferase domain of ACCase [6], suggesting this region as a prime site to expect resistance mutations. Within this domain region, we also found several residues that were different between the ACCase protein sequence of the PSR-TK and the HS-strain (Additional file 14, Additional file 15). However, it is unclear whether these substitutions play a role in resistance. Similarly to ACCase, other target site sequences (VGSC, cytB, GluCl, Rdl and CHS1) contained non-synonymous SNPs in the PSR-TK strain that were not present in the HS strain, but none were located at positions previously reported to be involved in acaricide resistance [42].

Differential expression analysis between an acaricide resistant (PSR-TK) and susceptible strain (HS) of P. ulmi using RNAseq data

RNA was extracted from 200 1–3 day old adult female P. ulmi mites from the PSR-TK and HS strains with four-fold biological replication (see Methods). Using the replicated RNAseq data and the DESeq2 software, we performed a differential expression analysis between the acaricide resistant (PSR-TK) and a susceptible (HS) strain of P. ulmi. Only those contigs that (1) were not considered as contamination by the NCBI contamination screen (2) have strand-specific reads and (3) have a BLASTx hit were included in our differential expression analysis (see Methods for more details). For all P. ulmi contigs the number of mapped reads per contig can be found in Additional file 16. For further analyses described below, we excluded P. ulmi contigs without a BLASTx hit as a biological role of the proteins coded by these contigs cannot be assigned. We observed that 123 of 8722 P. ulmi contigs were significantly upregulated (Benjamini-Hochberg adjusted p-value ≤ 0.05 and │FC│ > 2) in the PSR-TK strain compared to the HS strain, while 122 were significantly downregulated (Fig. 8, Tables 8 and 9, Additional file 17). For a subset of contigs (genes), the differential expression analysis based on RNAseq was validated by qPCR on cDNAs of the same populations (Fig. 9).

Fig. 8
figure8

Volcano plot of differentially expressed contigs in the acaricide resistant PSR-TK strain relative to the susceptible HS strain. The negative log10 of Benjamini-Hochberg (BH) adjusted p-values was plotted against the log2FC in gene expression. P. ulmi contigs up- and downregulated, by twofold or more (123 and 122 genes, respectively) are depicted as light grey circular dots. Circular dots of various colours depict contigs with a best BLASTx hit with T. urticae proteins putatively involved in xenobiotic detoxification ([36], see also Additional file 17). Contigs encoding detoxification genes putatively involved in spirodiclofen resistance are indicated in the plot by ID PuP450_18 (contig_01016, log2(FC) = 5.1, BH adjusted p-value = 1,1E−193) and ID PuCCE_7 (contig_00577, log2(FC) = 5.1, BH adjusted p-value =1,4E−257)

Table 8 Top 20 of upregulated contigs in the acaricide resistant P. ulmi PSR-TK strain relative to the susceptible HS strain. See Additional file 17 for a full list of differentially expressed contigs
Table 9 Top 20 of downpregulated contigs in the acaricide resistant P. ulmi PSR-TK strain relative to the susceptible HS strain. See Additional file 17 for a full list of differentially expressed contigs
Fig. 9
figure9

Validation of differentially expressed P. ulmi contigs by qPCR. Five upregulated and three downregulated contigs determined by differential gene expression analysis of RNA-seq data (see Fig. 8, Additional file 17) were selected for qPCR analysis. The data from qPCR were presented as mean of three replicates. Error bars represent the standard error of the calculated mean. Asterisks indicate significantly different expression values compared to the reference condition (HS)

Several up- and downregulated contigs showed high homology with sequences of arthropod bisegmented double stranded RNA viruses (biRNA viruses) (Tables 8 and 9, Additional file 17). There are several reports documenting viruses infecting Panonychus species [7477]. As such, differential expression of biRNA viruses between the two resistant strains might reflect the different virus composition and infection status.

Among the differentially expressed contigs, 29 coded for major detoxification enzyme families or potential new players in detoxification (see above and [36]). Among the upregulated contigs, we identified 4 CCEs, 3 CYPs, 3 GSTs, 3 ABCs, 4 UGTs, 2 lipocalins and 2 Major Facilitator Superfamily genes while 1 CYP, 2 CCEs, 1 GST, 2 ABCs and 2 IDRCDs were downregulated (Fig. 8, Additional file 17). Such a broad response in detoxifying enzyme families is typically also encountered in the spider mite T. urticae, where numerous gene-expression studies of susceptible and resistant strains have shown an overlap in the gene-families involved [20, 36, 78], but also specificity in term of the members of these families.

Two differentially expressed contigs also coded for the same small (~ 200 bp) fragment of the gene coding for ACCase, the target-site of cyclic keto-enols [6] (see above), but with a six nucleotide difference between the two contig sequences. One of these contigs was highly upregulated (contig_21107) in the PSR-TK strain according to our analysis while the other was highly downregulated (contig_08025). A close inspection of mapped reads revealed that sequences polymorphism could explain this result (see Additional file 16). Therefore, we mapped all reads from both P. ulmi strains against the manually assembled P. ulmi ACCase gene (see above) and found no significant difference in expression (data not shown). Two contigs encoding CCEs (contig_19626/PuCCE13 and contig_00577/PuCCE7) and a CYP (contig_01016/PuCYP_18) were among the most highly upregulated contigs while another CCE contig (contig_00445/PuCCE2) was one of the most downregulated (Fig. 9, Tables 8 and 9, Additional file 17). RT-qPCR data confirmed the upregulation of contig_00577/PuCCE7 and contig_01016/PuCYP18 and the downregulation of contig_00445/PuCCE2, while upregulation of contig_19626/PuCCE13 was only partially confirmed (Fig. 9). qPCR expression analysis also confirmed that the up- and downregulation of contig_00577/PuCCE7 and contig_00445/PuCCE2 in PSR-TK might be the result of spirodiclofen selection, with both contigs being slightly downregulated in Ge 16/09 (parental strain of PSR-TK) and highly up- and downregulated in PSR-TK relative to HS, respectively (Fig. 10). In addition, the same pattern was found for a CYP (contig_01016/PuP450_18) which is highly upregulated in PSR-TK-while and slightly downregulated in GE 16/09 compared to the HS strain (Fig. 10).

Fig. 10
figure10

Effect of spirodiclofen selection on expression of genes putatively involved in spirodiclofen resistance. qPCR quantification of expression levels of contig_01016 (PuP450_18), contig_00577 (PuCCE_7) and contig_00445(PuCCE_2) in P. ulmi strains HS (spirodiclofen susceptible strain), PSR-TK (spirodiclofen resistance ratio > 7000 relative to HS) and Ge 16/09, the parental strain of PSR-TK (spirodiclofen resistance ratio of 59 relative to HS) [17]. The data from qPCR were presented as mean of three replicates. Error bars represent the standard error of the calculated mean. Asterisks indicate significant different expression values compared to the reference condition (HS)

Previously, it was shown that the P450 inhibitor PBO and the esterase inhibitor DEF could enhance the toxicity of spirodiclofen in a highly spirodiclofen resistant T. urticae strain [18] and in a spirodiclofen resistant P. citri strain [19]. Similarly, Kramer and Nauen [17] showed that PBO enhanced the toxicity of spirodiclofen in PSR-TK P. ulmi strain, suggesting the involvement of CYPs in spirodiclofen resistance in both tetranychid species. The upregulation of a CYP (CYP392E10, tetur27g01030) and a CCE (TuCCE-04, tetur01g10750) was also strongly associated with spirodiclofen resistance in two unrelated T. urticae strains (SR-VP and SR-TK). Furthermore, functional expression confirmed that CYP392E10 metabolizes spirodiclofen [20]. Remarkably, the upregulated P. ulmi CCE and CYP in the PSR-TK strain cluster in different clades in our phylogenetic analyses of tetranychid CCEs and CYPs than those that are upregulated in the spirodiclofen-resistant T. urticae strains (Figs. 2 and 3). Thus, within the Tetranychidae, two different species are able to develop a similar resistance mechanism to spirodiclofen, but the enzymes recruited for this purpose do not need to have followed the same evolutionary path. This independent development of resistance in related species is further reinforced by the fact that only few of the differentially expressed contigs in the PSR-TK P. ulmi strain have a BLASTx hit with amino acid translations of differentially expressed genes in two spirodiclofen resistant T. urticae strains (SR-VP and SR-TK, [20]) (see Additional file 18). Overall, the analysis has pointed to a number of genes and enzymes that might be causal for the resistant phenotype, and future work including functional expression should confirm that the enzymes can metabolize or sequester spirodiclofen leading to resistance, as was recently undertaken for the model species T. urticae [20].

Conclusions

Using paired-end strand specific reads from RNA of spirodiclofen susceptible and resistant P. ulmi strains we generated a P. ulmi transcriptome dataset containing 27,777 contigs. As such, our transcriptome data represents a significant increase in the genomic resources that are available for this species. Phylogenetic analyses revealed that in the majority of cases detoxification gene (CCEs, CYPs, UGTs, ABCs) classes/clades/subfamilies are more numerate in T. urticae compared to P. ulmi, suggestive of a link between detoxification gene proliferation and the polyphagous nature of T. urticae. Specific radiations of subfamilies in P. ulmi were also observed. Annotation of all major target-sites in P. ulmi revealed the presence of mutations in AchE that are likely to confer carbamate and organophosphate resistance. Finally, we used the replicated RNAseq data to assess differences in gene expression between a spirodiclofen resistant and susceptible P. ulmi strain and found that a CYP and CCEs are likely to be associated with spirodiclofen resistance in P. ulmi. These results are in line with a previous report on molecular mechanisms of spirodiclofen resistance in T. urticae. However, the upregulated CYP and CCE in the resistant strains from both species seem not to have evolved from the same common ancestor, indicating both species developed spirodiclofen resistance in a similar but nevertheless independent evolutionary manner. To conclude, the new genomic resources and data that we present in this study for P. ulmi will substantially facilitate molecular studies of underlying mechanisms involved in acaricide resistance.

Methods

Mite strains

Mites from the three strains (HS, PSR-TK and Ge 16/09) used in this study were all identified morphologically as P. ulmi by Johan Witters (ILVO, Belgium) but also on the basis of cytochrome b (cytB) sequence similarity (cytB sequences of the HS and PSR-TK P. ulmi strain (Additional file 14) had their best BLASTx hit (E-value < 1E−100) with P. ulmi cytochrome b (YP_002808558.1) in the NCBI database. The HS strain was originally collected in 1990 from an apple orchard in Burscheid, Germany and is susceptible to most currently used acaricides, including spirodiclofen [17]. The Ge 16/09 is a field strain, collected in 2009 from an orchard in Heidenfahrt-Heidesheim, Germany, and characterised by high resistance levels to mite growth inhibitors hexythiazox and clofentezine (extrapolated resistance ratios in comparison with the susceptible HS strain were 20,000 and 3900 respectively) and moderate resistant to spirodiclofen with a reported spirodiclofen resistance ratio of 59 compared to the HS strain. The PSR-TK strain is a spirodiclofen-selected laboratory strain derived from Ge 16/09 and with a resistance ratio of more than 7000 compared to the HS strain [17]. All mite strains were reared under identical conditions in climatically controlled chambers at 24 +/− 1 °C, 60 % relative humidity and a 16:8 light:dark period for more than 10 generations on domestic plum trees (Prunus domestica L. var. Brompton). The PSR-TK strain was maintained under constant selection pressure by rearing on plum trees sprayed with 1000 mg a.i./L spirodiclofen until runoff.

RNA extraction, library construction, sequencing and assembly of the P. ulmi transcriptome

Total RNA was extracted from 200 1–3 day old adult female P. ulmi mites from the PSR-TK and HS strains using the RNEasy mini kit (Qiagen, Belgium) with four-fold biological replication (i.e., four replicates each for PSR-TK and HS). Each RNA sample was processed for sequencing by the sequencing service company Fasteris SA (Switzerland) according to the “HiSeq Service Stranded Standard Protocol”. The resulting library was sequenced by Fasteris SA (Switzerland) using the HiSeq 2000 Illumina technology generating strand-specific paired-reads of 2x100 bp. Adapter sequences were removed from the obtained 2x 100 bp reads by Fasteris SA and adapter trimmed reads were used for further analysis. Previously it was shown that a representative transcriptome assembly can be generated from a sub-sample of reads [40]. Hence, to reduce computation time, a subset (the first 12 million paired-reads/sample, 48 million reads/strain, 96 million paired-reads in total) of the total number of reads was used for P. ulmi transcriptome assembly using the CLC Genomics Workbench (CLC) software version 6.5.1 and default settings (mapping mode = fast, automatic bubble size = yes, minimum contig length = 200 nt, automatic word size = yes, perform scaffolding = yes, auto-detect paired distances = yes). Using the same subset of reads we also used an alternative approach to assemble the P. ulmi transcriptome. Reads were first trimmed using Sickle [79] with quality-cutoff set at 30 and length cutoff set at 90. Trimmed reads were assembled using the Velvet/Oases package [80] and with the following settings for Oases: −-kmin = 59 --kmax = 61 --kstep = 2 --merge = 61 –d “-shortPaired -strand_specific” -p “ins_length = 200”. Transcripts with maximum sequence length (30,044 transcripts) were filtered from the resulting Velvet/Oases assembly (44,903 transcripts) and used for further analysis. To compare the quality of both assemblies all RNA-seq reads (4 replicates/strain) were mapped using Bowtie version 2.1.0 [81] against either the CLC or Velvet/Oases assembly (see Differential expression analysis section below for mapping procedures).

Raw reads have been submitted to the NCBI Short Read Archive (SRA; experiment accession number SRX833872 (HS) and SRX833917 (PSR-TK)). The CLC assembly was submitted to the NCBI Transcriptome Shotgun Assembly (TSA) Sequence Database and deposited at DDBJ/EMBL/GenBank under the accession GCAC00000000. The version described in this paper is the first version, GCAC01000000. P. ulmi contigs containing more than 10 % of unassigned nucleotides (N) or more than 14 Ns in a row (1630 sequences) and contigs that were shorter than 200 nt after removal of vector contamination sequences (2 contigs) were not uploaded to the TSA Database following database submission procedures. After uploading, an NCBI contamination screen indicated that 227 out of the 26,145 uploaded sequences should be excluded as they showed very high identity with either rRNA of Bacteria, genes of Prunus sp. (host plant of P. ulmi strains) or Wolbachia sp. sequences. These sequences were removed from the TSA assembly and, finally, 25,918 contigs were uploaded. The complete collection of the CLC assembly (27,777 contigs) was added to this manuscript in Additional file 19. The (alternative) Velvet/Oases P. ulmi assembly was added as Additional file 13 to this manuscript.

Blast homology searches and sequence annotation

All P. ulmi contigs (27,777) were used for BLASTx searches against the NCBI non-redundant (nr) protein database (version of September 3rd 2014; E-value cut-off of 1E−5) using Standalone BLAST 2.2.30+ [82]. Sequences that did not yield BLASTx hits were subsequently searched against the non-redundant nucleotide database using BLASTn and with an e-value cut-off of 1E−5. As T. urticae is not yet included in the NCBI nr protein database, P. ulmi contigs were also used as query in BLASTx against the T. urticae proteome (version March 2014, http://bioinformatics.psb.ugent.be/orcae/overview/Tetur). All BLAST results were imported in Blast2GO (version 2.8) to assign gene ontology (GO) terms to sequences retrieved by BLAST search [83]. Further annotation was done with InterPro, where protein motifs were directly queried at the InterProScan web service and consequently merged to the existing annotation [84]. The annotation results were further fine-tuned with (1) the Annex function in order to augment the annotation through the inference of biological process and cellular component terms from molecular function annotations [85] and with (2) the generic GO slim reduction to summarize the functional information of the transcriptome dataset. Finally, to enable better visualisation of the results, Gene Ontology relationships and enzyme codes (ECs) were highlighted on the Kyoto Encyclopedia of Genes and Genomes (KEGG) maps.

Analysis of genes related to xenobiotic detoxification

Contigs encoding putative P. ulmi CYPs, CCEs, GSTs, UGTs and ABC proteins were retrieved from the P. ulmi transcriptome (CLC assembly) by a tBLASTn search and using protein sequences of T. urticae CYPs, CCEs, GSTs, UGTs and ABCs as query [21, 39, 57]. tBLASTn searches were performed with NCBI Standalone BLAST 2.2.30+ [82] and a cut-off E-value < 1E−5. Open reading frames (ORFs) of P. ulmi contigs encoding detoxification enzymes and ABC proteins were identified using “EMBOSS 6.3.1 getorf” integrated in the Mobyle portal framework (http://mobyle.pasteur.fr/). P. ulmi ORFs were aligned against each other using BLASTn. In cases where two ORFs showed more than 94 % identity at the nucleotide level, they were considered as allelic variants and the longest ORF was retained for phylogenetic analysis. ORFs with identical overlapping sequences were considered to be part of the same gene and merged using BioEdit version [86], except if the merged ORF misaligned when blasted to the T. urticae protein database.

Two previously published P. citri transcriptomes were also mined for transcripts encoding detoxification enzymes and ABC transporters. For both published P. citri transcriptomes, the cDNA library was constructed using a mixture of RNA extracted from different developmental stages and subsequently sequenced using Illumina HiSeq 2000 technology with a read-length of 90 bp [33, 34]. However, as the P. citri transcriptome sequence assemblies (TSAs) have not been made publicly available, we assembled these transcriptomes de novo using a similar approach as for P. ulmi and using all reads in the sequence read archives (SRAs) ERR044692-ERR044695 and SRR341928 provided by the studies of Liu et al. [33] and Niu et al. [34], respectively. P. citri ORFs encoding CYPs, CCEs, GSTs, UGTs, and ABC proteins were identified using a similar approach as for P. ulmi. Likewise, these P. citri ORFs were aligned against each other using BLASTn. In case two P. citri contigs showed at least 94 % identity they were considered as allelic variant and the longest ORF was retained for phylogenetic analysis.

A final selection of Panonychus ORFs encoding CYPs, GSTs UGTs and ABCs was translated into protein sequences and aligned with their homologues in T. urticae using MUSCLE version 3.8.31 [87]. Alignments were inspected by eye and for each alignment only protein sequences showing no misalignment and having a minimum ORF length (CYPs: 450 nt, CCEs: 450 nt, GSTs: 300 nt, UGTs: 450 nt) were retained for the final alignment. Two phylogenetic analyses were performed for the CYP gene family: one with all tetranychid CYPs and one restricted to tetranychid and D. melanogaster mitochondrial CYPs (see [88]). Panonychus mt CYPs were identified based on the phylogenetic analysis of all tetranychid CYPs. For the phylogenetic analysis of CCEs, a reference set of arthropod CCE protein sequences was also included in the alignment (Additional file 9). As N- and C termini of CCEs are highly variable, the alignment of CCE protein sequences was trimmed as previously described [49]. Only the nucleotide binding domain (NBD) of ABC proteins encoded by Panonychus ORFs was used for phylogenetic analysis of ABC proteins. NBDs were extracted using ScanProsite (http://prosite.expasy.org/scanprosite/) and the PROSITE profile PS5089. Next to T. urticae ABC transporter NBDs, C- and N terminal NBDs of D. melanogaster ABC transporters were also included in phylogenetic analysis of tetranychid ABC transporter NBDs [57]. Model selection was done with ProtTest 2.4 and according to the Akaike information criterion WAG + G + F, LG + I + G + F, LG + I + G + F, LG + I + G, WAG + G + F and LG + G + F were optimal for the phylogenetic reconstruction of CYPs, mitochondrial CYPs, UGTs, GSTs, CCEs and ABC proteins, respectively. Finally, for each alignment a maximum likelihood analysis was performed using Treefinder (version 2011) [89] bootstrapping with 1000 pseudoreplicates (LR-ELW). The resulting trees were midpoint rooted and edited with MEGA 6.0 software [90].

Annotation of horizontally transferred genes (excluding UGTs)

Annotation of P. ulmi and P. citri contigs encoding horizontally transferred genes (HTGs) was done in a similar way as for genes involved in xenobiotic detoxification. Model selection and phylogenetic analysis of IDRCD genes was done in a similar way as for genes involved in xenobiotic detoxification. According to the Akaike information criterion WAG + I + G + F was optimal for phylogenetic reconstruction of tetranychid IDRCDs.

Identification of P. ulmi target-site sequences

To obtain P. ulmi target site sequences as full length as possible we mined both P. ulmi transcriptomes (CLC and Velvet/Oases, see above) for contigs encoding target-sites of acaricides using tBLASTn (E-value cutoff < 1E−5) and T. urticae target-site protein sequences as query. P. ulmi contigs encoding target-sites (from both transcriptomes (CLC and Velvet/Oases) were aligned against each other and those contigs having identical overlapping regions were merged. The resulting contigs were manually assembled to obtain P. ulmi target-site sequences as full length as possible. All reads were then mapped against the P. ulmi target site sequences as for the differential expression analysis (see below) with this difference that resulting BAM files were merged for each strain using the SAMtools package [91]. For both the spirodiclofen susceptible and the resistant P. ulmi strain, a consensus sequence was then derived using a perl script (gene_extractor.pl) written by Rutger Vos and available at https://github.com/naturalis/fastq-simple-tools/tree/master/script. The obtained P. ulmi target site sequences were screened for all target site mutations previously described in the literature [42].

Differential expression analysis

To assess differential gene expression between the HS and PSR-TK strains, we first took the reverse complement of each P.ulmi contig where the majority of its BLASTx hits mapped on the reverse strand. P. ulmi contigs were subsequently concatenated into a single sequence with spacers (“Ns”) of length 100 added between contigs using a custom python script. For each RNA-seq replicate (4 per strain, 8 in total), reads were mapped against this concatenated sequence using Bowtie version 2.1.0 [81] with the preset option “--very-sensitive-local” and with the maximum fragment length for valid paired-end alignments set to 1000 (−X 1000). Mapped reads were counted using htseq-count that is included in the HTSeq package [92]. Forward strand reads were counted using the “--stranded = reverse” option while the “--stranded = yes” option was used to count reverse strand reads. Those P. ulmi contigs of which the proportion between the number of reads mapping to one strand and the total number of reads mapping to either the forward or the reverse strand was more than or equal to 0.95 were considered to have strand-specific reads (20,069 contigs, 72,3 % of all contigs). Read counts of P. ulmi contigs with (1) a BLAST hit, (2) strand-specific reads and (3) not considered contamination (see above) were used for differential expression analysis. Differentially expressed contigs between the resistant (PSR-TK) and susceptible (HS) P. ulmi strain were determined using the DESeq2 (version 1.6.3;[92]) and Bioconductor (http://bioconductor.org/) R-packages. The “unfiltered DESeq2 results” settings (dds < − DESeq(dds, minReplicatesForReplace = Inf) and res < − results(dds, cooksCutoff = FALSE, independent Filtering = FALSE)) were used for differential expression analysis. P. ulmi contigs with a fold change higher than two and a Benjamini-Hochberg adjusted p-value less than or equal to 0.05 were considered differentially expressed.

RT-qPCR

A set of eight differentially expressed contigs were evaluated by RT-qPCR. This set included 3 down- and 5 upregulated P. ulmi contigs. P. ulmi contig_01010 and contig_00741 showed best BLASTx hits with T. urticae RP49 (tetur18g03590) and ubiquitin (tetur06g00900), E-value of 4E−70 and 5E−59, respectively, and were used as reference genes for RT-qPCR. Primers were designed using Primer3 v.0.4.0 [93] and are listed in Additional file 20. RNA was extracted in triplicate from each strain as described above and 2 μg was reverse transcribed using the Maxima First Strand cDNA synthesis kit for RT-qPCR (Fermentas Life Sciences). All qPCR reactions were conducted with a thermal cycler Mx3005P (Stratagene). Reactions were prepared with Maxima SYBR Green qPCR/Master Mix following the manufacturer’s instructions (Fermentas Life Sciences) and run in two technical and three biological replicates. Non-template-controls were also included to eliminate potential contamination of the samples. The qPCR run protocol was as follows: initial denaturation at 95 °C for 10s, 35 cycles of 95 °C for 15 s, 55 °C for 30s, 72 °C for 30s. To exclude the possibility of nonspecific amplification, each PCR was followed by a melting curve step (ramping from 95 to 55 °C by 1 °C every 2 s) to confirm a single amplicon. Amplification efficiency of each primer pair was calculated from the standard curve, prepared of cDNA of all samples tested, with a dilution range from 0.4 to 50 ng RNA. The Ct values for P. ulmi contig_01010 (RP49) and contig_00741 (ubiquitin) were used for normalization. Significant differences in the relative expression values of the target genes were tested with pairwise fixed reallocation randomization [94, 95].

References

  1. 1.

    Jeppson LR, Keifer HH, Baker EW. Mites injurious to economic plants. Berkeley, USA: University of California Press; 1975.

  2. 2.

    Van Leeuwen T, Tirry L, Yamamoto A, Nauen R, Dermauw W. The economic importance of acaricides in the control of phytophagous mites and an update on recent acaricide mode of action research. Pestic Biochem Physiol. 2015;121:12–21.

  3. 3.

    Van Leeuwen T, Vontas J, Tsagkarakou A, Dermauw W, Tirry L. Acaricide resistance mechanisms in the two-spotted spider mite Tetranychus urticae and other important Acari: A review. Insect Biochem Mol Biol. 2010;40(8):563–72.

  4. 4.

    Whalon M, E., Mota-Sanchez R, M., Hollingworth RM, Duynslager L: Arthropod Pesticide Resistance Database. http://www.pesticideresistance.org/. Accessed May 2015. In. http://www.pesticideresistance.org/.

  5. 5.

    Bretschneider T, Benet-Buchholz J, Fischer R, Nauen R. Spirodiclofen and spiromesifen - Novel acaricidal and insecticidal tetronic acid derivatives with a new mode of action. Chimia. 2003;57(11):697–701.

  6. 6.

    Lummen P, Khajehali J, Luther K, Van Leeuwen T. The cyclic keto-enol insecticide spirotetramat inhibits insect and spider mite acetyl-CoA carboxylases by interfering with the carboxyltransferase partial reaction. Insect Biochem Mol Biol. 2014;55:1–8.

  7. 7.

    Elbert A, Brück E, Sone S, Toledo A. Worldwide uses of the new acaricide Envidor® in perennial crops. Pflanzenschutz-Nachr Bayer. 2002;55:287–304.

  8. 8.

    Wachendorff U, Nauen R, Schnorbach HJ, Rauch N, Elbert A. The biological profile of spirodiclofen (Envidor®) – a new selective tetronic acid acaricide. Pflanzenschutz-Nachr Bayer. 2002;55:149–76.

  9. 9.

    De Maeyer L, Geerinck R. The multiple target use of spirodiclofen (Envidor 240 SC) in IPM pomefruit in Belgium. Comm Agric Appl Biol Sci. 2009;74(1):225–32.

  10. 10.

    Nauen R, Stumpf N, Elbert A. Efficacy of BAJ 2740, a new acaricidal tetronic acid derivative, against tetranychid spider mite species resistant to conventional acaricides. Pests and diseases, Vol. 1. Brighton, UK: The BCPC Conference; 2000. p. 453-458

  11. 11.

    Rauch N, Nauen R. Spirodiclofen resistance risk assessment in Tetranychus urticae (Acari : Tetranychidae): a biochemical approach. Pestic Biochem Physiol. 2002;74(2):91–101.

  12. 12.

    Pree DJ, Whitty KJ, Van Driel L. Baseline susceptibility and cross resistances of some new acaricides in the European red mite, Panonychus ulmi. Exp Appl Acarol. 2005;37(3–4):165–71.

  13. 13.

    Van Leeuwen T, Stillatus V, Tirry L. Genetic analysis and cross-resistance spectrum of a laboratory-selected chlorfenapyr resistant strain of two-spotted spider mite (Acari : Tetranychidae). Exp Appl Acarol. 2004;32(4):249–61.

  14. 14.

    Van Leeuwen T, Van Pottelberge S, Tirry L. Comparative acaricide susceptibility and detoxifying enzyme activities in field-collected resistant and susceptible strains of Tetranychus urticae. Pest Manag Sci. 2005;61(5):499–507.

  15. 15.

    Thiel M, Nauen R. Investigations on acaricide resistance in populations of the European red mite, Panonychus ulmi Koch, collected at Lake Constance. Gesunde Pflanze. 2006;58(4):239–45.

  16. 16.

    Ilias A, Roditakis E, Grispou M, Nauen R, Vontas J, Tsagkarakou A. Efficacy of ketoenols on insecticide resistant field populations of two-spotted spider mite Tetranychus urticae and sweet potato whitefly Bemisia tabaci from Greece. Crop Prot. 2012;42:305–11.

  17. 17.

    Kramer T, Nauen R. Monitoring of spirodiclofen susceptibility in field populations of European red mites, Panonychus ulmi (Koch) (Acari: Tetranychidae), and the cross-resistance pattern of a laboratory-selected strain. Pest Manag Sci. 2011;67(10):1285–93.

  18. 18.

    Van Pottelberge S, Van Leeuwen T, Khajehali J, Tirry L. Genetic and biochemical analysis of a laboratory-selected spirodiclofen-resistant strain of Tetranychus urticae Koch (Acari: Tetranychidae). Pest Manag Sci. 2009;65(4):358–66.

  19. 19.

    Yu D, Wang CF, Yu CY, Huang Y, Yao J, Hu J. Laboratory selection for spirodiclofen resistance and cross-resistance in Panonychus citri. Afr J Biotechnol. 2011;10(17):3424–9.

  20. 20.

    Demaeght P, Dermauw W, Tsakireli D, Khajehali J, Nauen R, Tirry L, et al. Molecular analysis of resistance to acaricidal spirocyclic tetronic acids in Tetranychus urticae: CYP392E10 metabolizes spirodiclofen, but not its corresponding enol. Insect Biochem Mol Biol. 2013;43(6):544–54.

  21. 21.

    Grbic M, Van Leeuwen T, Clark RM, Rombauts S, Rouze P, Grbic V, et al. The genome of Tetranychus urticae reveals herbivorous pest adaptations. Nature. 2011;479(7374):487–92.

  22. 22.

    Van Leeuwen T, Demaeght P, Osborne EJ, Dermauw W, Gohlke S, Nauen R, et al. Population bulk segregant mapping uncovers resistance mutations and the mode of action of a chitin synthesis inhibitor in arthropods. Proc Natl Acad Sci U S A. 2012;109(12):4407–12.

  23. 23.

    Dermauw W, Ilias A, Riga M, Tsagkarakou A, Grbić M, Tirry L, et al. The cys-loop ligand-gated ion channel gene family of Tetranychus urticae: Implications for acaricide toxicology and a novel mutation associated with abamectin resistance. Insect Biochem Mol Biol. 2012;42(7):455–65.

  24. 24.

    Demaeght P, Osborne EJ, Odman-Naresh J, Grbic M, Nauen R, Merzendorfer H, et al. High resolution genetic mapping uncovers chitin synthase-1 as the target-site of the structurally diverse mite growth inhibitors clofentezine, hexythiazox and etoxazole in Tetranychus urticae. Insect Biochem Mol Biol. 2014;51:52–61.

  25. 25.

    Ilias A, Vontas J, Tsagkarakou A. Global distribution and origin of target site insecticide resistance mutations in Tetranychus urticae. Insect Biochem Mol Biol. 2014;48:17–28.

  26. 26.

    Wilhelm BT, Marguerat S, Watt S, Schubert F, Wood V, Goodhead I, et al. Dynamic repertoire of a eukaryotic transcriptome surveyed at single-nucleotide resolution. Nature. 2008;453(7199):1239–U1239.

  27. 27.

    Nagalakshmi U, Wang Z, Waern K, Shou C, Raha D, Gerstein M, et al. The transcriptional landscape of the yeast genome defined by RNA sequencing. Science. 2008;320(5881):1344–9.

  28. 28.

    Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Meth. 2008;5(7):621–8.

  29. 29.

    Vivancos AP, Gueell M, Dohm JC, Serrano L, Himmelbauer H. Strand-specific deep sequencing of the transcriptome. Genome Res. 2010;20(7):989–99.

  30. 30.

    Cabrera AR, Donohue KV, Khalil SMS, Scholl E, Opperman C, Sonenshine DE, et al. New approach for the study of mite reproduction: The first transcriptome analysis of a mite, Phytoseiulus persimilis (Acari: Phytoseiidae). J Insect Physiol. 2011;57(1):52–61.

  31. 31.

    Karatolos N, Pauchet Y, Wilkinson P, Chauhan R, Denholm I, Gorman K, et al. Pyrosequencing the transcriptome of the greenhouse whitefly, Trialeurodes vaporariorum reveals multiple transcripts encoding insecticide targets and detoxifying enzymes. BMC Genomics. 2011;12:56.

  32. 32.

    Pereira Firmino AA, de Assis Fonseca FC, Pepino de Macedo LL, Coelho RR, Antonino de Souza Jr JD, Togawa RC, et al. Transcriptome analysis in cotton boll weevil (Anthonomus grandis) and RNA interference in insect pests. Plos One. 2013;8(12):e85079.

  33. 33.

    Liu B, Jiang G, Zhang Y, Li J, Li X, Yue J, et al. Analysis of transcriptome differences between resistant and susceptible strains of the citrus red mite Panonychus citri (Acari: Tetranychidae). Plos One. 2011;6(12):e28516.

  34. 34.

    Niu JZ, Dou W, Ding TB, Shen GM, Zhang K, Smagghe G, et al. Transcriptome analysis of the citrus red mite, Panonychus citri, and its gene expression by exposure to insecticide/acaricide. Insect Mol Biol. 2012;21(4):422–36.

  35. 35.

    Bryon A, Wybouw N, Dermauw W, Tirry L, Van Leeuwen T. Genome wide gene-expression analysis of facultative reproductive diapause in the two-spotted spider mite Tetranychus urticae. BMC Genomics. 2013;14:815.

  36. 36.

    Dermauw W, Wybouw N, Rombauts S, Menten B, Vontas J, Grbic M, et al. A link between host plant adaptation and pesticide resistance in the polyphagous spider mite Tetranychus urticae. Proc Natl Acad Sci U S A. 2013;110(2):E113–22.

  37. 37.

    Wybouw N, Balabanidou V, Ballhorn DJ, Dermauw W, Grbic M, Vontas J, et al. A horizontally transferred cyanase gene in the spider mite Tetranychus urticae is involved in cyanate metabolism and is differentially expressed upon host plant change. Insect Biochem Mol Biol. 2012;42(12):881–9.

  38. 38.

    Wybouw N, Dermauw W, Tirry L, Stevens C, Grbic M, Feyereisen R, et al. A gene horizontally transferred from bacteria protects arthropods from host plant cyanide poisoning. eLife. 2014;3:e02365–5.

  39. 39.

    Ahn S-J, Dermauw W, Wybouw N, Heckel DG, Van Leeuwen T. Bacterial origin of a diverse family of UDP-glycosyltransferase genes in the Tetranychus urticae genome. Insect Biochem Mol Biol. 2014;50:43–57.

  40. 40.

    Francis W, Christianson L, Kiko R, Powers M, Shaner N, Haddock S. A comparison across non-model animals suggests an optimal sequencing depth for de novo transcriptome assembly. BMC Genomics. 2013;14(1):167.

  41. 41.

    Sparks TC, Nauen R. IRAC: Mode of action classification and insecticide resistance management. Pestic Biochem Physiol. 2015;121:122–8.

  42. 42.

    Feyereisen R, Dermauw W, Van Leeuwen T. Genotype to phenotype, the molecular and physiological dimensions of resistance in arthropods. Pest Biochemistry and Physiology 2015;121:61-77.

  43. 43.

    Despres L, David J-P, Gallet C. The evolutionary ecology of insect resistance to plant chemicals. Trends Ecol Evol. 2007;22(6):298–307.

  44. 44.

    Feyereisen R. Insect CYP Genes and P450 Enzymes. In: Gilbert LI, editor. Insect molecular biology and biochemistry. London: Academic Press, Elsevier; 2012. p. 236–316.

  45. 45.

    Riga M, Tsakireli D, Ilias A, Morou E, Myridakis A, Stephanou EG, et al. Abamectin is metabolized by CYP392A16, a cytochrome P450 associated with high levels of acaricide resistance in Tetranychus urticae. Insect Biochem Mol Biol. 2014;46:43–53.

  46. 46.

    Ran C, Jiang G-F, Liu B, Liu H-Q, Li H-J, Wang J-J. Selection of resistance to amitraz and analysis of expression difference of cytochrome P450 genes in Panonychus citri (Acari: Tetranychidae). Acta Entomol Sin. 2012;55(6):703–9.

  47. 47.

    Ding T-B, Niu J-Z, Yang L-H, Zhang K, Dou W, Wang J-J. Transcription profiling of two cytochrome P450 genes potentially involved in acaricide metabolism in citrus red mite Panonychus citri. Pestic Biochem Physiol. 2013;106(1–2):28–37.

  48. 48.

    Oakeshott JG. Biochemical genetics and genomics of insect esterases. In: Gilbert LI, editor. Comprehensive Molecular Insect Science – Pharmacology. 5th ed. Oxford: Elsevier; 2005. p. 309–81.

  49. 49.

    Claudianos C, Ranson H, Johnson RM, Biswas S, Schuler MA, Berenbaum MR, et al. A deficit of detoxification enzymes: pesticide sensitivity and environmental response in the honeybee. Insect Mol Biol. 2006;15(5):615–36.

  50. 50.

    Zhang K, Niu J-Z, Ding T-B, Dou W, Wang J-J. Molecular characterization of two carboxylesterase genes of the citrus red mite, Panonychus citri (Acari: Tetranychidae). Arch Insect Biochem Physiol. 2013;82(4):213–26.

  51. 51.

    Ortelli F, Rossiter LC, Vontas J, Ranson H, Hemingway J. Heterologous expression of four glutathione transferase genes genetically linked to a major insecticide-resistance locus from the malaria vector Anopheles gambiae. Biochem J. 2003;373:957–63.

  52. 52.

    Enayati AA, Ranson H, Hemingway J. Insect glutathione transferases and insecticide resistance. Insect Mol Biol. 2005;14(1):3–8.

  53. 53.

    Reddy BPN, Prasad GBKS, Raghavendra K. In silico analysis of glutathione S-transferase supergene family revealed hitherto unreported insect specific delta- and epsilon-GSTs and mammalian specific mu-GSTs in Ixodes scapularis (Acari: Ixodidae). Comput Biol Chem. 2011;35(2):114–20.

  54. 54.

    Xu Z, Zhu W, Liu Y, Liu X, Chen Q, Peng M, et al. Analysis of insecticide resistance-related genes of the carmine spider mite Tetranychus cinnabarinus based on a de novo assembled transcriptome. Plos One. 2014;9(5):e94779.

  55. 55.

    Pavlidi N, Tseliou V, Riga M, Nauen R, Van Leeuwen T, Labrou NE, Vontas J. Functional characterization of glutathione S-transferases associated with insecticide resistance in Tetranychus urticae. Pestic Biochem Physiol. 2015;121:53–60.

  56. 56.

    Dermauw W, Van Leeuwen T. The ABC gene family in arthropods: Comparative genomics and role in insecticide transport and resistance. Insect Biochem Mol Biol. 2014;45:89–110.

  57. 57.

    Dermauw W, Osborne EJ, Clark RM, Grbic M, Tirry L, Van Leeuwen T. A burst of ABC genes in the genome of the polyphagous spider mite Tetranychus urticae. BMC Genomics. 2013;14:317.

  58. 58.

    Rowland A, Miners JO, Mackenzie PI. The UDP-glucuronosyltransferases: Their role in drug metabolism and detoxification. Int J Biochem Cell Biol. 2013;45(6):1121–32.

  59. 59.

    Brattsten LB. Metabolic defenses against plant allelochemicals. In: Rosenthal GA, Berenbaum MR, editors. Herbivores Their interactions with secondary plant metabolites. 2nd ed. London: Academic; 1992. p. 175–242.

  60. 60.

    Ahn S-J, Badenes-Perez FR, Reichelt M, Svatos A, Schneider B, Gershenzon J, et al. Metabolic detoxification of capsaicin by UDP-glycosyltransferases in three Helicoverpa species. Arch Insect Biochem Physiol. 2011;78(2):104–18.

  61. 61.

    Daimon T, Hirayama C, Kanai M, Ruike Y, Meng Y, Kosegawa E, et al. The silkworm green b locus encodes a quercetin 5-O-glucosyltransferase that produces green cocoons with UV-shielding properties. Proc Natl Acad Sci U S A. 2010;107(25):11471–6.

  62. 62.

    Lee SW, Ohta K, Tashiro S, Shono T. Metabolic resistance mechanisms of the housefly (Musca domestica) resistant to pyraclofos. Pestic Biochem Physiol. 2006;85(2):76–83.

  63. 63.

    Sasai H, Ishida M, Murakami K, Tadokoro N, Ishihara A, Nishida R, et al. Species-specific glucosylation of DIMBOA in larvae of the rice armyworm. Biosci Biotechnol Biochem. 2009;73(6):1333–8.

  64. 64.

    Migeon A, Dorkeld F. Spider Mites Web: a comprehensive database for the Tetranychidae. Available at: http://www1.montpellier.inra.fr/CBGP/spmweb/. 2015.

  65. 65.

    Vetter J. Plant cyanogenic glycosides. Toxicon. 2000;38(1):11–36.

  66. 66.

    Tsagkarakou A, Van Leeuwen T, Khajehali J, Ilias A, Grispou M, Williamson MS, et al. Identification of pyrethroid resistance associated mutations in the para sodium channel of the two-spotted spider mite Tetranychus urticae (Acari: Tetranychidae). Insect Mol Biol. 2009;18(5):583–93.

  67. 67.

    Van Nieuwenhuyse P, Van Leeuwen T, Khajehali J, Vanholme B, Tirry L. Mutations in the mitochondrial cytochrome b of Tetranychus urticae Koch (Acari: Tetranychidae) confer cross-resistance between bifenazate and acequinocyl. Pest Manag Sci. 2009;65(4):404–12.

  68. 68.

    Van Leeuwen T, Vanholme B, Van Pottelberge S, Van Nieuwenhuyse P, Nauen R, Tirry L, et al. Mitochondrial heteroplasmy and the evolution of insecticide resistance: Non-Mendelian inheritance in action. Proc Natl Acad Sci U S A. 2008;105(16):5980–5.

  69. 69.

    Kwon DH, Clark JM, Lee SH. Cloning of a sodium channel gene and identification of mutations putatively associated with fenpropathrin resistance in Tetranychus urticae. Pestic Biochem Physiol. 2010;97(2):93–100.

  70. 70.

    Kwon DH, Im JS, Ahn JJ, Lee J-H, Marshall Clark J, Lee SH. Acetylcholinesterase point mutations putatively associated with monocrotophos resistance in the two-spotted spider mite. Pestic Biochem Physiol. 2010;96(1):36–42.

  71. 71.

    Kwon DH, Yoon KS, Clark JM, Lee SH. A point mutation in a glutamate-gated chloride channel confers abamectin resistance in the two-spotted spider mite, Tetranychus urticae Koch. Insect Mol Biol. 2010;19(4):583–91.

  72. 72.

    Khajehali J, Van Leeuwen T, Grispou M, Morou E, Alout H, Weill M, et al. Acetylcholinesterase point mutations in European strains of Tetranychus urticae (Acari: Tetranychidae) resistant to organophosphates. Pest Manag Sci. 2010;66(2):220–8.

  73. 73.

    Karatolos N, Williamson MS, Denholm I, Gorman K, ffrench-Constant R, Nauen R. Resistance to spiromesifen in Trialeurodes vaporariorum is associated with a single amino acid replacement in its target enzyme acetyl-coenzyme A carboxylase. Insect Mol Biol. 2012;21(3):327–34.

  74. 74.

    Smith KM, Hills GJ, Munger F, Gilmore JE. A suspected virus disease of the citrus red mite Panonychus citri (McG.). Nature. 1959;184(4679):70–0.

  75. 75.

    Bird FT. A virus disease of the European red mite Panonychus ulmi (Koch). Can J Microbiol. 1967;13(8):1131–1.

  76. 76.

    Reed DK, Desjardins PR. Morphology of a non-occluded virus isolated from citrus red mite, Panonychus citri. Experientia. 1982;38(4):468–9.

  77. 77.

    Putman WL. Occurrence and transmission of a virus disease of the european red mite, Panonychus ulmi. The Canadian Entomologist. 1970;102(03):305–21.

  78. 78.

    Khalighi M, Dermauw W, Wybouw N, Bajda S, Osakabe M, Tirry L, Van Leeuwen T. Molecular analysis of cyenopyrafen resistance in the two-spotted spider mite Tetranychus urticae. Pest Manag Sci. 2015:n/a-n/a.

  79. 79.

    Joshi NA, Fass JN: Sickle. A sliding-window, adaptive, quality-based trimming tool for FastQ files (Version 1.33). Available at https://github.com/najoshi/sickle. 2011.

  80. 80.

    Schulz MH, Zerbino DR, Vingron M, Birney E. Oases: robust de novo RNA-seq assembly across the dynamic range of expression levels. Bioinformatics. 2012;28(8):1086–92.

  81. 81.

    Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Meth. 2012;9(4):357–U354.

  82. 82.

    Shiryev SA, Papadopoulos JS, Schaeffer AA, Agarwala R. Improved BLAST searches using longer words for protein seeding. Bioinformatics. 2007;23(21):2949–51.

  83. 83.

    Conesa A, Gotz S, Garcia-Gomez JM, Terol J, Talon M, Robles M. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics. 2005;21(18):3674–6.

  84. 84.

    Hunter S, Jones P, Mitchell A, Apweiler R, Attwood TK, Bateman A, et al. InterPro in 2011: new developments in the family and domain prediction database. Nucleic Acids Res. 2012;40(10):4725–5.

  85. 85.

    Myhre S, Tveit H, Mollestad T, Laegreid A. Additional Gene Ontology structure for improved biological reasoning. Bioinformatics. 2006;22(16):2020–7.

  86. 86.

    Hall TA. BioEdit: a user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT. Nucleic Acids Symp Ser. 1999;41:95–8.

  87. 87.

    Edgar RC. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 2004;32(5):1792–7.

  88. 88.

    Good RT, Gramzow L, Battlay P, Sztal T, Batterham P, Robin C. The molecular evolution of cytochrome P450 genes within and between Drosophila species. Genome Biol Evol. 2014;6(5):1118–34.

  89. 89.

    Jobb G, von Haeseler A, Strimmer K. TREEFINDER: a powerful graphical analysis environment for molecular phylogenetics. BMC Evol Biol. 2004;4:18.

  90. 90.

    Tamura K, Stecher G, Peterson D, Filipski A, Kumar S. MEGA6: Molecular evolutionary genetics analysis version 6.0. Mol Biol Evol. 2013;30(12):2725–9.

  91. 91.

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

  92. 92.

    Anders S, McCarthy DJ, Chen Y, Okoniewski M, Smyth GK, Huber W, et al. Count-based differential expression analysis of RNA sequencing data using R and Bioconductor. Nat Protoc. 2013;8(9):1765–86.

  93. 93.

    Rozen S, Skaletsky H. Primer3 on the WWW for general users and for biologist programmers. Methods Mol Biol. 2000;132:365–86.

  94. 94.

    Pfaffl MW. A new mathematical model for relative quantification in realtime RT-PCR. Nucleic Acids Res. 2001;29(9):e45.

  95. 95.

    Pfaffl MW, Horgan GW, Dempfle L. Relative expression software tool (REST (c)) for group-wise comparison and statistical analysis of relative expression results in real-time PCR. Nucleic Acids Res. 2002;30(9):e36.

Download references

Acknowledgements

WD is a postdoctoral fellow of the Fund for Scientific Research Flanders. This work was partially supported by the Research Foundation Flanders (FWO) (grant G.0093.12 N to LT and TVL). RG was supported by National Institutes of Health Genetics Training Grant T32 GM07464. We thank Jan van Arkel for providing a photograph of an adult female of P. ulmi. We thank Johan Witters for morphological identification of P. ulmi species used in this study.

Author information

Correspondence to Wannes Dermauw or Thomas Van Leeuwen.

Additional information

Competing interests

The authors declare they have no competing interests.

Author’s contributions

SB, WD, RMC and TVL designed research. SB and WD performed experiments. WD, SB and RG analysed data. SB, WD and TVL wrote the manuscript, with input from RMC, RG, RN and LT. All authors read and approved the final manuscript.

Sabina Bajda and Wannes Dermauw contributed equally to this work.

Additional files

Additional file 1: Table S1.

Alignment summary for the CLC and Velvet/Oases assemblies. (XLSX 16 kb)

Additional file 2: File S1.

Complete transcriptome assembly of P. citri reassembled from raw reads data [33]. (Format: FASTA, examination of this file requires the BioEdit program available at http://www.mbio.ncsu.edu/bioedit/bioedit.html. (FAS 19163 kb)

Additional file 3: File S2.

Complete transcriptome assembly of P. citri reassembled from raw reads data [34]. (Format: FASTA, examination of this file requires the BioEdit program available at http://www.mbio.ncsu.edu/bioedit/bioedit.html). (FAS 21895 kb)

Additional file 4: Table S2.

BLAST2GO analysis of P. ulmi transcriptome. (XLSX 12041 kb)

Additional file 5: Figure S1.

Taxonomic distribution of top BLAST hits of the P. ulmi transcriptome. Analysis of the taxonomic distribution of the 11,673 top BLAST hits obtained by BLASTx and BLASTn against the non-redundant protein and nucleotide database (NCBI) respectively, and by BLASTx against the T. urticae proteome (http://bioinformatics.psb.ugent.be/orcae/overview/Tetur). Wheel diagram depicts the percentage distribution within different taxonomic groups. (PDF 25 kb)

Additional file 6: Figure S2.

Enzyme Classification (EC) analysis of the P. ulmi transcriptome. The wheel diagram depicts percentage distribution of EC numbers in general EC terms. (PDF 26 kb)

Additional file 7: Table S3.

Annotation details and accession numbers of Cytochrome P450s. (XLSX 154 kb)

Additional file 8: Figure S3.

Phylogenetic analysis of P. ulmi mitochondrial CYPs. Maximum likelihood phylogenetic tree of P. ulmi, P. citri, T. urticae and D. melanogaster mitochondrial CYP protein sequences. The tree was rooted with T. urticae CYP307A1, P. ulmi PcP450_42 and P. citri Pc0cP450_21 (CYP2 clan members). The scale bar represents 0.5 amino-acid substitutions per site. Numbers at the branch point of each node represent the bootstrap value resulting from 1000 pseudoreplicates (LR-ELW). Colour and shape codes are as follows: P. ulmi, green triangle, P. citri, yellow triangle, T. urticae, red square and D. melanogaster, black dot. (PDF 1.35 MB)

Additional file 9: Table S4.

Annotation details and accession numbers of the carboxylcholinesterases. (XLSX 171 kb)

Additional file 10: Table S5.

Annotation details and accession numbers of Glutathione S- transferases. (XLSX 57 kb)

Additional file 11: Table S6.

Annotation details and accession numbers of ATP-binding cassette family transporters. (XLSX 441 kb)

Additional file 12: Table S7.

Annotation details and accession numbers of HTGs. (XLSX 159 kb)

Additional file 13: File S3.

Alternative transcriptome assembly (Velvet/Oases) of P. ulmi used for obtaining target-site sequences as full length as possible and as an extra resource for mining of HTGs. (Format: FASTA, examination of this file requires the BioEdit program available at http://www.mbio.ncsu.edu/bioedit/bioedit.html. (FAS 20474 kb)

Additional file 14: File S4.

Target site nucleotide sequences of the P. ulmi PSR-TK and HS-strain. (Format: FASTA, examination of this file requires the BioEdit program available at http://www.mbio.ncsu.edu/bioedit/bioedit.html). (FAS 67 kb)

Additional file 15: Figure S4.

ClustalW alignment of ACCase protein sequences of the PSR-TK and HS P. ulmi strain and T. urticae. The carboxyltransferase domain is shaded in grey while the residue corresponding to the E645K substitution associated with spiromesifen resistance in T. vaporarorium is shaded in yellow. (DOC 261 kb)

Additional file 16: Table S8.

The number of mapped RNA reads per contig for both P. ulmi strains (HS and PSR-TK). Mapped reads were counted using htseq-count, included in the HTSeq python package [92]. (XLSX 2544 kb)

Additional file 17: Table S9.

Differentially expressed contigs in the P. ulmi PSR-TK strain relative to the HS strain (Benjamini-Hochberg adjusted p-value ≤ 0.05 and │FC│ > 2). (XLSX 30 kb)

Additional file 18: Table S10.

Differentially expressed contigs in the P. ulmi PSR-TK strain with a BLASTx (E-value treshold e-15) hit against amino acid translations of differentially expressed genes in SR-VP and SR-TK, two spirodiclofen resistant T. urticae strains [20]. (XLSX 48 kb)

Additional file 19: File S5.

Transcriptome assembly (CLC Genomics Workbench) of P. ulmi described in this study. (Format: FASTA, examination of this file requires the BioEdit program available at http://www.mbio.ncsu.edu/bioedit/bioedit.html). (FAS 28273 kb)

Additional file 20: Table S11.

qPCR primers used in this study. Sequences are given for the forward (F) and reverse primer (R). (DOCX 13 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. 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.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Keywords

  • Cyclic keto-enols
  • Tetranychidae
  • Target-site mutation
  • biRNA virus
  • Spirotetramat
  • Intradiol ring-cleavage dioxygenases
  • UDP-glycosyltransferases
  • Cytochrome P450 mono-oxygenases
  • Carboxyl/cholinesterases
  • Gluthathione-S-transferases