CA Tabuloc, CR Carlson, F Ganjisaffar, CC Truong, C Chen, KM Lewald, S Hidalgo, N Nicola, C Jones, AA Sial, FG Zalom, JC Chiu

Data description

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 pyrethroid insecticide (specifically zeta-cypermethrin) or spinosad 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 (ZT16). 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.

Zeta-cypermethrin-susceptible lines = S7, S8
Zeta-cypermethrin-resistant lines = S3, S4
Spinosad-susceptible lines = C2, C5
Spinosad-resistant lines = C3, C4

Differential gene expression (DEG analysis)

Differential gene expression analysis was performed using sequencing reads derived from Illumina 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.

Information for data tables

-Suppl. Table 1: Statistics for zeta-cypermethrin bioassays

Bioassays were performed to identify isofemale lines resistant to zeta-cypermethrin (Mustang® Maxx). Eight isofemale lines were tested (indicated as S#). Two isofemale lines from an untreated orchard (Wolfskill, W#) served as the susceptible control. Each biological replicate consisted of 5 males and 5 females (n=8). Significance was determined by one-way ANOVA followed by Tukey’s multiple comparison test.

-Suppl. Table 2: Statistics for spinosad bioassays

Bioassays were performed to identify isofemale lines resistant to spinosad (Entrust®). Eight isofemale lines were tested (indicated as C#). Two isofemale lines from an untreated orchard (Wolfskill, W#) served as the susceptible control. Each biological replicate consisted of 5 males and 5 females (n=8). Significance was determined by one-way ANOVA followed by Tukey’s multiple comparison test.

-Suppl. Table 4: DEGs in zeta-cypermethrin-resistant Drosophila suzukii isofemale line (line S3) relative to susceptible lines (lines S7 and S8) developed from the same field-collected population

List of differentially expressed genes (DEGs) in zeta-cypermethrin resistant line S3 as compared to susceptible lines S7 and S8. “SWD gene ID” and “SWD Gene Description” are based off of NCBI Drosophila suzukii Annotation Release 102 based on the LBDM_Dsuz_2.1.pri assembly (accession no. GCF_013340165.1) (Paris et al., 2020). “Fbgn” is the FlyBase Gene ID that is most similar (by sequence similarity) to the SWD gene ID. “Dmel gene ID” is the Drosophila melanogaster gene symbol associated with each Dmel gene ID. “Log2FoldChange” is the difference in expression between the resistant line (S3) to both susceptible lines (S7 and S8–pooled together). Positive values indicate that the gene is upregulated in the resistant line while negative values indicate that the gene is downregulated in the resistant line. “padj” is the Benjamini-Hochberg adjusted p-value, and -log10(padj) is the log transformation of the adjusted p-values, where higher values indicate more significant fold change differences. Data is sorted into 3 tabs: genes that are upregulated, downregulated, and not significant. Genes with fold change differences between resistant vs susceptible populations with a Benjamini-Hochberg adjusted p-value < 0.05 were considered differentially expressed.

-Suppl. Table 5: DEGs in zeta-cypermethrin-resistant Drosophila suzukii isofemale line (line S4) relative to susceptible lines (lines S7 and S8) developed from the same field-collected population

List of differentially expressed genes (DEGs) in zeta-cypermethrin resistant line S4 as compared to susceptible lines S7 and S8. “SWD gene ID” and “SWD Gene Description” are based off of NCBI Drosophila suzukii Annotation Release 102 based on the LBDM_Dsuz_2.1.pri assembly (accession no. GCF_013340165.1) (Paris et al., 2020). “Fbgn” is the FlyBase Gene ID that is most similar (by sequence similarity) to the SWD gene ID. “Dmel gene ID” is the Drosophila melanogaster gene symbol associated with each Dmel gene ID. “Log2FoldChange” is the difference in expression between the resistant line (S3) to both susceptible lines (S7 and S8–pooled together). Positive values indicate that the gene is upregulated in the resistant line while negative values indicate that the gene is downregulated in the resistant line. “padj” is the Benjamini-Hochberg adjusted p-value, and -log10(padj) is the log transformation of the adjusted p-values, where higher values indicate more significant fold change differences. Data is sorted into 3 tabs: genes that are upregulated, downregulated, and not significant. Genes with fold change differences between resistant vs susceptible populations with a Benjamini-Hochberg adjusted p-value < 0.05 were considered differentially expressed.

-Suppl. Table 6: Statistics comparing the expression of metabolic genes in zeta-cypermethrin-resistant vs. susceptible Drosophila suzukii

Relative expression (FPKM) of cytochrome P450 genes (Cyp), heat shock proteins (Hsp), the carboxylesterase Cricklet (Clt), and glutathione-s-transferase E3 (GstE3) in the susceptible (S7 and S8) vs resistant (S3 and S4) groups. There are a total of 3 biological replicates of 8-10 females per line. Significance was determined by 2-way ANOVA.

-Suppl. Table 7: Enrichment of DEGs in zeta-cypermethrin-resistant Drosophila suzukii (line S3)

Enrichment pathways within the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) Biological Processes (bio proc) categories for genes up- or down-regulated in line S3 (Suppl. Table 4). The table was produced by ShinyGo (Ge et al., 2020). “Enrichment FDR” (where FDR is the false discovery rate) is adjusted from the hypergeometric test. “nGenes” is number of differentially expressed genes that belong to the indicated pathway while “Pathway Genes” indicate the total number of genes that make up that pathway. FDR tells us how likely the enrichment is by chance. “Fold Enrichment” indicates how drastically genes of a certain pathway is overrepresented. “Pathway” is the name of the enriched pathway, and “URL” provides information on that pathway. Finally, “Genes” provides a list of differentially expressed genes that are in each pathway.

-Suppl. Table 8: Enrichment of DEGs in zeta-cypermethrin-resistant Drosophila suzukii (line S4)

Enrichment pathways within the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) Biological Processes (bio proc) categories for genes up- or down-regulated in line S4 (Suppl. Table 5). The table was produced by ShinyGo (Ge et al., 2020). “Enrichment FDR” (where FDR is the false discovery rate) is adjusted from the hypergeometric test. “nGenes” is number of differentially expressed genes that belong to the indicated pathway while “Pathway Genes” indicate the total number of genes that make up that pathway. FDR tells us how likely the enrichment is by chance. “Fold Enrichment” indicates how drastically genes of a certain pathway is overrepresented. “Pathway” is the name

-Suppl. Table 9: Genes in all WGCNA modules for zeta-cypermethrin-resistant Drosophila suzukii (line S3)

To identify which genes in line S3 are most correlated with zeta-cypermethrin resistance, WGCNA was performed. Modules with a p-value < 0.05 were considered significant. The data table was produced by the WGCNA package (Version 1.72.1) (Langfelder & Horvath, 2008) on R. “Transcript” is the NCBI Reference Sequence ID, and “LOC” is the SWD gene ID based off of NCBI Drosophila suzukii Annotation Release 102 based on the LBDM_Dsuz_2.1.pri assembly (accession no. GCF_013340165.1) (Paris et al., 2020). “Fbgn” is the FlyBase Gene ID that is most similar (by sequence similarity) to the SWD gene ID. “GeneSymbol” is the Drosophila melanogaster gene symbol associated with each Fbgn. “Module Color” is the clusters of highly interconnected genes. Modules correspond to clusters of genes with high absolute correlations. Gene significance (GS) is the correlation between gene expression and insecticide susceptibility. “p.GS.Susceptibility” is the p-value associate with the GS. Module membership (MM) is the association between gene expression and each module eigengene. The eigengene can be thought of as a weighted average expression profile. Finally, “p.MM” is the p-value associated with the Module Membership.

-Suppl. Table 10: Enrichment of the turquoise module in zeta-cypermethrin-resistant Drosophila suzukii (line S3)

Enrichment pathways within the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) Biological Processes (bio proc) categories for genes within the turquoise cluster. Turquoise was the cluster most correlated with zeta-cypermethrin resistance in line S3. The table was produced by ShinyGo (Ge et al., 2020). “Enrichment FDR” (where FDR is the false discovery rate) is adjusted from the hypergeometric test. “nGenes” is number of differentially expressed genes that belong to the indicated pathway while “Pathway Genes” indicate the total number of genes that make up that pathway. FDR tells us how likely the enrichment is by chance. “Fold Enrichment” indicates how drastically genes of a certain pathway is overrepresented. “Pathway” is the name of the enriched pathway, and “URL” provides information on that pathway. Finally, “Genes” provides a list of differentially expressed genes that are in each pathway.

-Suppl. Table 11: Genes in all WGCNA modules for zeta-cypermethrin-resistant Drosophila suzukii (line S4)

To identify which genes in line S4 are most correlated with zeta-cypermethrin resistance, WGCNA was performed. Modules with a p-value < 0.05 were considered significant. Data table was produced by the WGCNA package (Version 1.72.1) (Langfelder & Horvath, 2008) on R. “Transcript” is the NCBI Reference Sequence ID, and “LOC” is the SWD gene ID based off of NCBI Drosophila suzukii Annotation Release 102 based on the LBDM_Dsuz_2.1.pri assembly (accession no. GCF_013340165.1) (Paris et al., 2020). “Fbgn” is the FlyBase Gene ID that is most similar (by sequence similarity) to the SWD gene ID. “GeneSymbol” is the Drosophila melanogaster gene symbol associated with each Fbgn. “Module Color” is the clusters of highly interconnected genes. Modules correspond to clusters of genes with high absolute correlations. Gene significance (GS) is the correlation between gene expression and insecticide susceptibility. “p.GS.Susceptibility” is the p-value associate with the GS. Module membership (MM) is the association between gene expression and each module eigengene. The eigengene can be thought of as a weighted average expression profile. Finally, “p.MM” is the p-value associated with the Module Membership.

-Suppl. Table 12: Enrichment of the turquoise module in zeta-cypermethrin-resistant Drosophila suzukii (line S4)

Enrichment pathways within the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) Biological Processes (bio proc) categories for genes within the turquoise cluster. Turquoise was the cluster most correlated with zeta-cypermethrin resistance in line S4. The table was produced by ShinyGo (Ge et al., 2020). “Enrichment FDR” (where FDR is the false discovery rate) is adjusted from the hypergeometric test. “nGenes” is number of differentially expressed genes that belong to the indicated pathway while “Pathway Genes” indicate the total number of genes that make up that pathway. FDR tells us how likely the enrichment is by chance. “Fold Enrichment” indicates how drastically genes of a certain pathway is overrepresented. “Pathway” is the name of the enriched pathway, and “URL” provides information on that pathway. Finally, “Genes” provides a list of differentially expressed genes that are in each pathway.

-Suppl. Table 13: DEGs in spinosad-resistant Drosophila suzukii isofemale line (line C3) relative to susceptible lines (lines C2 and C5) developed from the same field-collected population

List of differentially expressed genes (DEGs) in zeta-cypermethrin resistant line C3 as compared to susceptible lines C2 and C5. “SWD gene ID” and “SWD Gene Description” are based off of NCBI Drosophila suzukii Annotation Release 102 based on the LBDM_Dsuz_2.1.pri assembly (accession no. GCF_013340165.1) (Paris et al., 2020). “Fbgn” is the FlyBase Gene ID that is most similar (by sequence similarity) to the SWD gene ID. “Dmel gene ID” is the Drosophila melanogaster gene symbol associated with each Dmel gene ID. “Log2FoldChange” is the difference in expression between the resistant line (C3) to both susceptible lines (C2 and C5–pooled together). Positive values indicate that the gene is upregulated in the resistant line while negative values indicate that the gene is downregulated in the resistant line. “padj” is the Benjamini-Hochberg adjusted p-value, and -log10(padj) is the log transformation of the adjusted p-values, where higher values indicate more significant fold change differences. Data is sorted into 3 tabs: genes that are upregulated, downregulated, and not significant. Genes with fold change differences between resistant vs susceptible populations with a Benjamini-Hochberg adjusted p-value < 0.05 were considered differentially expressed.

-Suppl. Table 14: DEGs in spinosad-resistant Drosophila suzukii isofemale line (line C4) relative to susceptible lines (lines C2 and C5) developed from the same field-collected population

List of differentially expressed genes (DEGs) in zeta-cypermethrin resistant line C4 as compared to susceptible lines C2 and C5. “SWD gene ID” and “SWD Gene Description” are based off of NCBI Drosophila suzukii Annotation Release 102 based on the LBDM_Dsuz_2.1.pri assembly (accession no. GCF_013340165.1) (Paris et al., 2020). “Fbgn” is the FlyBase Gene ID that is most similar (by sequence similarity) to the SWD gene ID. “Dmel gene ID” is the Drosophila melanogaster gene symbol associated with each Dmel gene ID. “Log2FoldChange” is the difference in expression between the resistant line (C4) to both susceptible lines (C2 and C5–pooled together). Positive values indicate that the gene is upregulated in the resistant line while negative values indicate that the gene is downregulated in the resistant line. “padj” is the Benjamini-Hochberg adjusted p-value, and -log10(padj) is the log transformation of the adjusted p-values, where higher values indicate more significant fold change differences. Data is sorted into 3 tabs: genes that are upregulated, downregulated, and not significant. Genes with fold change differences between resistant vs susceptible populations with a Benjamini-Hochberg adjusted p-value < 0.05 were considered differentially expressed.

-Suppl. Table 15: Statistics comparing the expression of metabolic and cuticular genes in spinosad-resistant vs. susceptible Drosophila suzukii

Relative expression (FPKM) of metabolic and cuticular genes (twdl: tweedle; cpr: cuticular protein) in the susceptible (C2 and C5) and resistant (C3 and C4) groups.There are a total of 3 biological replicates of 8-10 females per line. Significance was determined by 2-way ANOVA.

-Suppl. Table 16: Enrichment of DEGs in spinosad-resistant Drosophila suzukii (line C3)

Enrichment pathways within the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) Biological Processes (bio proc) categories for genes up- or down-regulated in line C3 (Suppl. Table 13). The table was produced by ShinyGo (Ge et al., 2020). “Enrichment FDR” (where FDR is the false discovery rate) is adjusted from the hypergeometric test. “nGenes” is number of differentially expressed genes that belong to the indicated pathway while “Pathway Genes” indicate the total number of genes that make up that pathway. FDR tells us how likely the enrichment is by chance. “Fold Enrichment” indicates how drastically genes of a certain pathway is overrepresented. “Pathway” is the name of the enriched pathway, and “URL” provides information on that pathway. Finally, “Genes” provides a list of differentially expressed genes that are in each pathway.

-Suppl. Table 17: Enrichment of DEGs in spinosad-resistant Drosophila suzukii (line C4)

Enrichment pathways within the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) Biological Processes (bio proc) categories for genes up- or down-regulated in line C4 (Suppl. Table 14). The table was produced by ShinyGo (Ge et al., 2020). “Enrichment FDR” (where FDR is the false discovery rate) is adjusted from the hypergeometric test. “nGenes” is number of differentially expressed genes that belong to the indicated pathway while “Pathway Genes” indicate the total number of genes that make up that pathway. FDR tells us how likely the enrichment is by chance. “Fold Enrichment” indicates how drastically genes of a certain pathway is overrepresented. “Pathway” is the name of the enriched pathway, and “URL” provides information on that pathway. Finally, “Genes” provides a list of differentially expressed genes that are in each pathway.

-Suppl. Table 18: Genes in all WGCNA modules for spinosasd-resistant Drosophila suzukii (line C3)

To identify which genes in line C3 are most correlated with spinosad resistance, WGCNA was performed. Modules with a p-value < 0.05 were considered significant. Data table was produced by the WGCNA package (Version 1.72.1) (Langfelder & Horvath, 2008) on R. “Transcript” is the NCBI Reference Sequence ID, and “LOC” is the SWD gene ID based off of NCBI Drosophila suzukii Annotation Release 102 based on the LBDM_Dsuz_2.1.pri assembly (accession no. GCF_013340165.1) (Paris et al., 2020). “Fbgn” is the FlyBase Gene ID that is most similar (by sequence similarity) to the SWD gene ID. “GeneSymbol” is the Drosophila melanogaster gene symbol associated with each Fbgn. “Module Color” is the clusters of highly interconnected genes. Modules correspond to clusters of genes with high absolute correlations. Gene significance (GS) is the correlation between gene expression and insecticide susceptibility. “p.GS.Susceptibility” is the p-value associate with the GS. Module membership (MM) is the association between gene expression and each module eigengene. The eigengene can be thought of as a weighted average expression profile. Finally, “p.MM” is the p-value associated with the Module Membership.

-Suppl. Table 19: Genes in all WGCNA modules for spinosasd-resistant Drosophila suzukii (line C4)

To identify which genes in line C4 are most correlated with spinosad resistance, WGCNA was performed. Modules with a p-value < 0.05 were considered significant. Data table was produced by the WGCNA package (Version 1.72.1) (Langfelder & Horvath, 2008) on R. “Transcript” is the NCBI Reference Sequence ID, and “LOC” is the SWD gene ID based off of NCBI Drosophila suzukii Annotation Release 102 based on the LBDM_Dsuz_2.1.pri assembly (accession no. GCF_013340165.1) (Paris et al., 2020). “Fbgn” is the FlyBase Gene ID that is most similar (by sequence similarity) to the SWD gene ID. “GeneSymbol” is the Drosophila melanogaster gene symbol associated with each Fbgn. “Module Color” is the clusters of highly interconnected genes. Modules correspond to clusters of genes with high absolute correlations. Gene significance (GS) is the correlation between gene expression and insecticide susceptibility. “p.GS.Susceptibility” is the p-value associate with the GS. Module membership (MM) is the association between gene expression and each module eigengene. The eigengene can be thought of as a weighted average expression profile. Finally, “p.MM” is the p-value associated with the Module Membership.

-Suppl. Table 20: Enrichment of the green module in spinosasd-resistant Drosophila suzukii (line C4)

Enrichment pathways within the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) Biological Processes (bio proc) categories for genes within the green cluster. Green was the cluster most correlated with spinosad resistance in line C4. The table was produced by ShinyGo (Ge et al., 2020). “Enrichment FDR” (where FDR is the false discovery rate) is adjusted from the hypergeometric test. “nGenes” is number of differentially expressed genes that belong to the indicated pathway while “Pathway Genes” indicate the total number of genes that make up that pathway. FDR tells us how likely the enrichment is by chance. “Fold Enrichment” indicates how drastically genes of a certain pathway is overrepresented. “Pathway” is the name of the enriched pathway, and “URL” provides information on that pathway. Finally, “Genes” provides a list of differentially expressed genes that are in each pathway.

-Suppl. Table 21: Statistics for quantitative PCR (qPCR)-based assays with 2022 D. suzukii populations

Results of a one sample T and Wilcoxon Test comparing the expression of cytochrome P450 (cyp) genes, tweedle (twdl) genes, and ecdysone receptor (ecR) in 2022 field-collected populations compared to the average gene expression of the susceptible (S7 and C2) and resistant (S3 and C3) isofemale lines (established from 2019 collections) (n=5 biological replicates of 8-10 females). One-way ANOVA followed by Holm-Sidak’s multiple comparisons test was also performed to assess expression differences between the 2022 F1 populations and the 2019 resistant controls.

-Suppl. Table 22: Statistics for bioassays with D. melanogaster mutants

To confirm knockdown of genes associated with insecticide resistance, the relative expression of cyp4d14, cyp4d8, and cpr6D in wild-type (w1118) and mutant fly lines was compared to a hypothetical mean of 1 using a one-sample t and Wilcoxon Test (n=4). Next, mortality of w1118 and mutant flies when exposed to sublethal concentrations of zeta-cypermethrin and spinosad (n=5 biological replicates of 5 males and 5 females) were compared using two-way ANOVA followed by Holm-Sidak’s multiple comparisons test.

Note: Missing values are denoted by NA.

Sharing/Access

Raw Data can be found under BioProject PRJNA983428