3.1 Selection, resistance and cross resistance to emamectin benzoate
Laboratory selections were conducted with emamectin benzoate on strain GR-IER-16-6 for eight sequential generations. Doses were adjusted to result in 60%-80% mortality in all cases, maintaining similar levels of selection pressure throughout. The probit analysis results are shown on Table 2. The slope ranged between 0.73 and 1.48 indicating low level variability in the response of the tested strains. In all cases low chi-square values were observed indicating good fit of the probit model to the actual bioassay data. The estimated LC50 value for emamectin benzoate in the susceptible reference strain ES-Sus was at 0.09 mg L-1. The selected strain GR-IER-16-6-S8 exhibited LC50 5.46 mgL-1, indicating a 60-fold resistance ratio (RR) and a substantial 7-fold increase compared to the parental strain GR-IER-16-6 (LC50 0.81 mgL-1, RR 9-fold, Table 2).
The cross-resistance pattern against the other registered avermectin insecticide, abamectin, was investigated (Table 2). The susceptible reference strain ES-Sus exhibited an LC50 at 0.05 mgL-1, the field collected strain GR-IER-16-6 0.24mgL-1 while the selected strain GR-IER-16-6-S8 0.66 mgL-1. The resistance ratio to abamectin increased after selection with emamectin benzoate from 4-fold to 13-fold, suggesting a moderate but clear cross resistance pattern between the two chemical compounds.
3.2 Effect of synergists
The probit analysis results from the synergistic bioassays are presented in Table 3. Pre-exposure to PBO resulted in significant reduction of LC50 levels to emamectin benzoate. More specifically, the LC50 was reduced from 5.46 mgL-1 to 0.91 mgL-1 after exposure to PBO, resulting in a 6-fold synergism ratio indicating that P450s monooxygenases may play a key role in emamectin benzoate toxicity and/or resistance. Pre-exposure to DEM resulted in partially synergized emamectin benzoate resistance (3.8-fold synergism ratio). Analogous synergistic effects were not observed in the susceptible reference strain ES-Sus. When the esterase inhibitor DEF was used, a low synergistic ratio of 1.6-fold was observed, indicating no significant effects and minimal involvement of esterase EST in the phenotype (Table 3).
3.3 Investigation of target site resistance mechanisms
The glutamate-gated chloride channel (GluCl), the gamma amino butyric acid (GABA) gated chloride channel and the histamine-gated chloride channel (HisCl) interact with avermectins (Clark et al. 1995; Ludmerer et al. 2002; Wolstenholme and Rogers 2005; Prichard et al. 2012) although GluCl is considered the main target-site in arthropods (Sparks and Nauen 2015). Several point mutations in GluCl (G314D/G326E, I32IT) (Kwon et al. 2010; Dermauw et al. 2012; Wang et al. 2016a; Xue et al. 2020), and GABA (I281T, A301S/G/N, T305L, V332I, T350M and R357Q) reviewed by Feyereisen et al. (2015) have been previously tightly associated in arthropods with abamectin and cyclodien/phenylpyrazoles resistance respectively. Gene fragments of GluCl and GABA (transmembrane domains TM1-TM3) encompassing known SNPs were amplified and sequenced in order to investigate whether the presence of point mutations could explain the resistance phenotype. From the aforementioned mutations in all three strains the A301S mutation was identified in segregation, indicating that this mutation is unlikely to be associated with emamectin benzoate resistance in the studied strains. In all the other previously published positions all strains harbored the wild type allele (Figure 1). Additionally, the 36-bp deletion in GluCl, previously associated with abamectin resistance in Plutella xylostella, was detected in all three strains, while a non-synonymous SNP at position 310 (E to V, T. absoluta numbering) was found in segregation and low frequency in GABA of the 16-6-P and the ES-Sus strains but not in the 16-6-S strain. Hence, their involvement in emamectin benzoate resistance is highly unlikely.
3.4 RNA sequencing
Resistance of T. absoluta to emamectin benzoate was studied at the molecular level using RNA sequencing (RNAseq) on three populations of different resistance status; (a) the laboratory susceptible strain (ES-Sus), (b) the parental strain derived from the population GR-IER-16-6, which was maintained in the laboratory without emamectin benzoate selection pressure (16-6-P), and (c) the resistant strain selected with emamectin benzoate for eight consecutive generations (16-6-S) (Table 2). The Principal Components Analysis (PCA) revealed that 21.44% of the total variation could be explained by principal component 1 (PC1) while 16.91% could be explained by PC2 (Supplementary File 2). The replicates from these three strains are clearly separated from each other, further confirming their quality and also the meaningfulness of downstream comparisons. It should be noted that replicates from populations 16-6-P and 16-6-S are only separated on PC2. In contrast, ES-Sus is obviously different from either 16-6-P or 16-6-S on both PC1 and PC2. Such finding is expected due to the fact that 16-6-S and 16-6-P share the same genetic background (originating from the same populations GR-IER-16-6, whereas the strain ES-Sus originates from Spain.
All RNAseq reads were pooled and assembled de novo with Trinity (Grabherr et al. 2011) into 549,601 transcripts >200 bp, originating from 233,453 unigenes. Running BUSCO (Waterhouse et al. 2018) on the transcriptome assembly showed that it is fairly complete, since it contains the complete sequence of 79.5% Insecta BUSCOs (Supplementary File 3). TransDecoder predicted proteins for 92,134 of the transcripts, with the majority of them having a significant hit against the Uniref50 database (Table 4). A set of unigenes was obtained after removing the following three types of proteins; (a) non-lepidopteran, contaminating sequences, (b) multiple isoforms of a unigene, and (c) proteins that were entirely contained within other, larger proteins. This final gene set included 32,502 proteins and, according to BUSCO, contained 75% of the Insecta BUSCO, only 4.3% of which was duplicated. An additional 16.5% of the BUSCOs were found as fragmented and 8.1% were missing (Supplementary File 4). As a result of this strict filtering scheme, the number of missing genes increased in the gene set. However, this was an acceptable trade-off in order to obtain a very clean unigene set.
A total of 143 P450s were identified in the aforementioned unigene set of T. absoluta. This number is higher than the 114 P450s previously identified in the genome of another lepidopteran insect, the cotton bollworm Helicoverpa armigera (Pearce et al. 2017). Such an increase is explained by the inclusion of P450 fragments in T. absoluta. We opted for including fragmented P450s because the present study has an exploratory nature and we aim at giving an exhaustive list of all candidate P450s. In this data set, representatives of all four major insect P450 clades were found. The majority of T. absoluta P450s (62 out of 143) belonged to Clan 4, 57 to Clan 3, 13 to the mitochondrial Clan, and 11 to Clan 2 (Table 5). The phylogenetic analysis allowed the classification of the T. absoluta P450s into clans and even into specific families, using the P450s from H. armigera as a reference (Figure 2). Moreover, it enabled the identification of P450s that are either duplicated or missing in T. absoluta, compared to H. armigera. More specifically, it appears that there is at least one H. armigera P450 (CYP4CG14) that is duplicated in T. absoluta. Moreover, it is worth mentioning that no CYP18B1 ortholog was found in T. absoluta (Figure 2), most probably due to low expression levels.
3.5 Differential expression summary
Differential expression analysis on the entire set of 233,453 unigenes showed that 4,199 unigenes were significantly over-expressed (log2|FC| >2, FDR <0.001) in the 16-6-S strain (emamectin benzoate-selected) against either 16-6-P (parental), or ES-Sus (reference susceptible) (Supplementary File 5). Similarly, another 4,887 were significantly under-expressed in either one of the same comparisons (Supplementary File 6). This list includes two cytochrome P450s (DN75966_c4_g1, DN60698_c0_g2), two GSTs (DN60153_c3_g2, DN65736_c3_g1) and three ABC transporters (DN48209_c0_g2, DN52519_c0_g1, DN80907_c5_g2) (Table 6). It should be noted that the majority of these unigenes do not contain a protein-coding gene, according to the TransDecoder results (see above). More specifically, only 350 of the over-expressed unigenes contain a protein-coding gene, of which 290 have a significant similarity to another protein in the Uniref50 database. Similarly, there are 937 under-expressed unigenes that contain a predicted peptide, of which 859 have a significant similarity in Uniref50.
Among the 4,199 over-expressed unigenes there were 15 that (a) are commonly over-expressed in 16-6-S against either 16-6-P, or ES-Sus, (b) are full-length, and (c) have a significant hit in Uniref50 (Supplementary File 5). The majority of these genes are similar to uncharacterized proteins (n =11), whereas the remaining four include a trypsin and a ribosomal protein. Using the same criteria we found 278 unigenes in the 4,887 under-expressed unigenes (Supplementary File 6).
Moreover, differential expression analysis revealed that two P450s were significantly over-expressed (log2|FC| >2, FDR <0.001) in the 16-6-S resistant strain (Table 6). One of these genes (DN75966_c4_g1) is similar to CYP9G5 and therefore belongs to Clan 3, whereas the other (DN60698_c0_g2) is similar to CYP4AU1 and therefore belongs to Clan 4 (Figure 2). Despite an almost 32-fold over-expression for the DN75966_c4_g1 in the selected (16-6-S) compared to the parental one (16-6-P) strain, it should be noted that it is fragmented because the ORF encodes is only 113 amino acids (Table 6), In contrast, DN60698_c0_g2 appears to be full-length.
3.6 RNAseq qPCR P450 validation
Among the seven detoxification genes found over-expressed in the 16-6-S resistant strain we selected two P450 genes for qPCR validation; DN75966_c4_g1 and DN60698_c0_g2. The former was significantly over-expressed against the parental 16-6-P strain, whereas the latter was over-expressed compared to the Es-Sus strain. Quantitive PCR confirmed the levels of expression of DN60698_c0_g2 (log2FC=2.29, p-value <0.001) and supported the significant up-regulation calculated from the transcriptomic analysis (log2FC=2.212, p-value <0.001; Table 6). In the contrary, DN75966_c4_g1 up-regulation in 16-6-S vs 16-6-P was not confirmed via qPCR.