Transcriptome analysis of Drosophila suzukii reveals molecular mechanisms conferring pyrethroid and spinosad resistance

Drosophila suzukii possess a serrated ovipositor that enables them to lay eggs in soft-skinned, ripening fruits, making this insect a serious threat to berry production. Since its 2008 introduction into North America, growers have used insecticides as the primary approach for D. suzukii management, resulting in detections of insecticide resistance in this pest. This study sought to identify the molecular mechanisms conferring insecticide resistance in these resistant populations. We sequenced the transcriptomes of two pyrethroid- and two spinosad-resistant isofemale lines. In both pyrethroid-resistant lines and one spinosad-resistant line, we identified overexpression of metabolic genes that are implicated in resistance in other insect pests. In the other spinosad-resistant line, we observed an overexpression of cuticular genes that have been linked to resistance. Our findings enabled the development of molecular diagnostics that we used to confirm persistence of insecticide resistance in California. To validate these findings, we leveraged D. melanogaster mutants deficient in either metabolic or cuticular genes that were upregulated in resistant D. suzukii to demonstrate that these genes are involved in promoting resistance. This study is the first to characterize the molecular mechanisms of insecticide resistance in D. suzukii and provides insights into how current management practices can be optimized.

Methods
Data generated from Illumina short-read sequencing (PE150)

RNA was extracted from isofemale lines developed from unsprayed Drosophila suzukii collected from populations resistant to either spinosad insecticide or resistant insecticide on either strawberries (S) or caneberries (C).  Female D. suzukii flies were entrained at 25C in 12-hour light:12-hour dark cycles for two full days. On the third day, flies were collected on dry ice sixteen hours after lights-on. This time point was selected because D. suzukii was previously observed to exhibit a low level of cytochrome P450 expression at this time (Hamby et al., 2013). This means any overexpression may be more easily observed. Fly bodies were separated from heads using frozen metal sieves (Newark Wire Cloth Company, Clifton, New Jersey). Eight to ten female bodies were used per sample

Pyrethroid-susceptible lines = S7, S8

Pyrethroid-resistant lines = S3, S4

Spinosad-susceptible lines = C2, C5

Spinosad-resistant lines = C3, C4

Differential gene expression analysis was performed using sequencing reads derived from llumina short-read sequencing. First, rRNA reads were removed using SortMeRNA v2.1 (Kopylova et al., 2012). Adapters (ILLUMINACLIP parameters 2:30:10) and low-quality ends (LEADING: 10, TRAILING:10, MINLEN:36) were trimmed using Trimmomatic v0.35 (Bolger et al., 2014). Cleaned reads were aligned to the NCBI Drosophila suzukii Annotation Release 102 based on the LBDM_Dsuz_2.1.pri assembly (accession no. GCF_013340165.1) (Paris et al., 2020) using STAR v2.7.9a (Dobin et al., 2013). Count data from STAR (--quantMode GeneCounts) served as input in the DESeq2 package (Love et al., 2014) in R to perform differential expression analysis on each resistant line vs all susceptible samples. Each resistant line was compared to susceptible samples separately as each line might exhibit resistance due to different mechanisms. Genes with fold change differences between resistant vs susceptible populations with a Benjamini-Hochberg adjusted p-value < 0.05 were considered differentially expressed. Expression levels of genes were also measured as fragments per kilobase of exon per million mapped (FPKM) values calculated with Stringtie v2.0.4 (Pertea et al., 2015). The consistency between biological replicates was calculated with Pearson’s correlation coefficient, which was determined with the ‘stats’ package in R version 4.2.1. Expression differences of key genes between the resistant and susceptible populations were calculated with two-way ANOVA followed by two-stage linear set-up procedure of Benjamini, Krieger, and Yekutieli on GraphPad Prism. 

Weighted Gene Co-expression Network Analysis

Gene expression (in FPKM) served as input for Weighted Gene Co-expression Network Analysis (WGCNA). Genes with an expression value of zero for more than six samples were excluded from analysis. To explore the modules most correlated with insecticide resistance, a correlation analysis using resistance status was performed with the WGCNA package (Version 1.72.1) (Langfelder & Horvath, 2008) on R. Modules with a p-value < 0.05 were considered significant. Functional enrichment analysis (described below) was performed on the module with the highest correlation with resistance. 

Functional enrichment analysis

Genes were functionally annotated using BLAST against the NCBI Drosophila melanogaster Annotation r6.32 based on the Release 6 plus ISO1 mitochondrial genome assembly (accession no. GCA_000001215.4) (dos Santos et al., 2015). Gene Ontology (GO) enrichment of genes were performed using ShinyGO 0.76.3 (Ge et al., 2020). GO terms and pathways were considered enriched if the false discovery rate (FDR) < 0.05. 

Raw Data can be found under BioProject PRJNA983428