Transcriptomics-Based Approach Identifies Spinosad-Associated Targets in the Colorado Potato Beetle, Leptinotarsa decemlineata

Simple Summary The Colorado potato beetle Leptinotarsa decemlineata is a potato pest that can cause substantial damages to potato crops worldwide. Multiple approaches have been leveraged to control this pest including the use of a variety of insecticides. Resistance to different insecticides aimed at controlling this insect has been reported and much work has been conducted in recent years to elucidate the underlying molecular changes associated with insecticide resistance in L. decemlineata. However, information is sparse regarding the molecular impact associated with spinosad treatment in this insect pest. The current study thus explores transcriptional changes associated with spinosad response in L. decemlineata exposed to this compound using high-throughput sequencing. Results presented show multiple transcripts of interest that exhibit differential expression in spinosad-treated L. decemlineata and provide a preliminary footprint of transcripts affected by this insecticide in this potato pest. Select targets identified in this signature should be further explored in follow-up studies to better characterize their contribution, if any, in the process of spinosad resistance. Abstract The Colorado potato beetle Leptinotarsa decemlineata is an insect pest that threatens potato crops globally. The primary method to control its damage on potato plants is the use of insecticides, including imidacloprid, chlorantraniliprole and spinosad. However, insecticide resistance has been frequently observed in Colorado potato beetles. The molecular targets and the basis of resistance to imidacloprid and chlorantraniliprole have both been previously quantified. This work was undertaken with the overarching goal of better characterizing the molecular changes associated with spinosad exposure in this insect pest. Next-generation sequencing was conducted to identify transcripts that were differentially expressed between Colorado potato beetles exposed to spinosad versus control insects. Results showed several transcripts that exhibit different expression levels between the two conditions, including ones coding for venom carboxylesterase-6, chitinase 10, juvenile hormone esterase and multidrug resistance-associated protein 4. In addition, several microRNAs, such as miR-12-3p and miR-750-3p, were also modulated in the investigated conditions. Overall, this work reveals a molecular footprint underlying spinosad response in Colorado potato beetles and provides novel leads that could be targeted as part of RNAi-based approaches to control this insect pest.


Introduction
The Colorado potato beetle Leptinotarsa decemlineata (Say) (Coleoptera: Chrysomelidae) can significantly harm potato crops and is an insect pest of substantial economic importance for the agricultural industry worldwide [1]. Complete yield loss of potato crops was reported when defoliation was greater than 75%, underscoring the damage that this insect can cause to potato fields [2]. Several approaches, including crop rotation [3,4], biopesticides [5], and transgenic plants [6] have been leveraged to control this insect with various levels of success. The use of insecticides remains a primary means by which Colorado potato beetles are targeted even though greater resistance against select compounds has been observed in recent years. Decreased susceptibility against the neonicotinoids thiamethoxam and clothianidin was reported in populations of Colorado potato beetles collected in Canada [7]. Resistance to the neonicotinoid imidacloprid has also been reported in several studies [8,9]. Work on a population of Colorado potato beetles collected in Long Island notably displayed a 309-fold increase in imidacloprid resistance when compared with control insects [10]. Interestingly, imidacloprid-resistant insects sampled from the same location displayed low levels of resistance towards the compound spinosad.
Spinosad, an insecticide obtained following fermentation of the actinomycete Saccharopolyspora spinosa and comprising spinosyns A and D [11] has raised interest for its use in controlling Colorado potato beetles, even though the complete molecular signature linked with its use in this insect has yet to be determined. Spinosad is particularly of interest to organic growers and has been registered to control several insect pests in multiple countries on certified organic horticultural crops. However, resistance to spinosad has been documented in an increasing number of insects. Spinosad resistance has been found for the tomato borer Tuta absoluta (Lepidoptera: Gelechiidae) [12], the diamondback moth Plutella xylostella (Lepidoptera: Plutellidae) [13], the leafworm moth Spodoptera litura (Lepidoptera: Noctuidae) [14] and the fall armyworm Spodoptera frugiperda (Lepidoptera: Noctuidae) [15] to name a few. Pioneering work has shown that mutations in subunit α6 of the nicotinic acetylcholine receptor (nAChR) were associated with increased spinosad resistance in Drosophila melanogaster (Diptera: Drosophilidae) [16]. Work undertaken in P. xylostella exhibiting different degrees of response to spinosad showed substantial levels of γ-aminobutyric acid receptor (GABAR) transcripts in the spinosad-resistant strain when compared with a sensitive counterpart [17]. Spinosad treatment of S. frugiperda Sf9 cells was associated with induction of autophagy as well as programmed cell death [18,19]. High-throughput approaches have also been leveraged to better understand the molecular changes underlying spinosad resistance in insects. For example, whole transcriptome analysis conducted for strains of olive fly Bactrocera oleae (Diptera: Tephritidae) that were sensitive versus resistant to spinosad revealed differentially expressed genes associated with differential energy requirements in the resistant strain [20]. While these studies provide evidence of the molecular impact resulting from spinosad resistance in diverse pests, information is sparse regarding the effect of this compound on the Colorado potato beetle.
In this study, we undertook a transcriptomics-based approach to identify differentially expressed transcripts related to spinosad response in L. decemlineata. We quantified differentially expressed mRNAs and miRNAs in Colorado potato beetles treated with spinosad as compared to controls. We present a rationale for the modulation of these targets in the context of insecticide responses in Colorado potato beetles.

Insects
L. decemlineata adults were obtained from the Fredericton Research and Development Centre colony in Fredericton (NB, Canada) in July 2018. Insects were brought back to Université de Moncton (Moncton, NB, Canada) in plastic containers supplemented with potato (Solanum tuberosum var. Shepody) leaves. Upon arrival, insects were housed in a cage with potato plants until use in experiments. A group of Insects 2020, 11, 820 3 of 16 insects was acclimated in an incubator (Thermo Fisher Scientific, Waltham, MA, USA) set at 25 • C for 5 days under 16 L:8 D cycles and containing potato plants. A volume of 0.5 µL spinosad (Sigma #33706, dissolved in acetone, Sigma-Aldrich, St. Louis, MO, USA) at a concentration of 1 µg/µL was pipetted topically on the abdomen of 15 beetles. An equal volume of acetone was applied to 15 beetles used as controls. Insects were placed back into the incubator for 4 h. All insects displayed signs of activity after the incubation period and prior to storage. It is relevant to point out that more insects (15) were treated and stored for each condition than the number of insects required to generate the results presented in this specific study. Additional insects were stored for use in follow-up work relevant to spinosad response in L. decemlineata. Parallel exposure tests were also performed on groups of five L. decemlineata exposed to doses of 0.5, 1.0, 2.0, 5.0, 10.0, 20.0 and 40.0 µg of spinosad. Activity impairment and insect lethality was observed after 24 h in ≥40% of insects at the five highest doses, and the lowest dose (0.5 µg) was thus selected for this study (Supplementary Materials Table S1). Insect mobility impairment was assessed via moderate agitation of the dish and subsequent evaluation of the ability of the insect to right itself [21]. Insects were sampled by being placed rapidly in liquid nitrogen and stored at −80 • C until use.

RNA Isolation
Small and large RNA fractions from control and spinosad-treated L. decemlineata were generated, as previously reported [22], with the miRVana miRNA Isolation Kit (Thermo Fisher Scientific) following manufacturer's instructions and using two insects as the starting material for each replicate. Briefly, insects were homogenized in lysis/binding buffer using a mechanical homogenizer. Extraction was conducted using an acid-phenol/chloroform solution. Filter cartridges were used to obtain fractions enriched with small RNAs as well as fractions containing large RNAs depleted of small RNAs. A NanoVue Plus Spectrophotometer (Thermo Fisher Scientific) was used to measure RNA concentrations. Small and large RNA isolates were stored at −80 • C until analysis by next-generation sequencing for the spinosad-associated expression of miRNAs and mRNAs, respectively.

Large RNA Library Construction and Sequencing
Large RNA fractions (n = 3) generated with two insects per replicate were quantified and RNA integrity was analyzed on an Agilent TapeStation (Agilent, Santa Clara, CA, USA). Samples with RNA integrity numbers above 8 were selected for library preparation. Poly(A) enrichment was performed from 20 µg of total RNA to enrich for mRNA using a Dynabeads mRNA Direct Micro Purification Kit (Thermo Fisher Scientific). An amount of 100 ng of poly(A)-enriched RNA was fragmented with RNase III and purified using the magnetic bead purification module Ion Total RNA-Seq v2 kit (Thermo Fisher Scientific). Size distribution of the fragmented RNA was assessed on the Agilent TapeStation using High Sensitivity RNA ScreenTape. A total of 50 ng of fragmented poly(A)-enriched RNA was used to prepare whole transcriptome libraries using Ion Total RNA-Seq v2 kit (Thermo Fisher Scientific). Poly(A) transcripts were ligated onto Ion Torrent adapters followed by reverse transcription. Resulting cDNA transcripts were purified and barcoded libraries were generated using Ion Xpress RNA barcodes (Thermo Fisher Scientific). Conditions recommended by Ion Torrent were followed throughout the protocol. Yield and size distribution of libraries were analyzed on an Agilent TapeStation using D1000 ScreenTape. Libraries were then equally pooled by pooling 7 pM libraries from each sample, and amplified onto ion sphere particles (ISPs) supplied by the Ion PI Hi-Q OT2 kit (Thermo Fisher Scientific). ISPs enriched with template libraries were loaded onto the Ion PI chip v3 and next-generation sequencing (NGS) was performed on an Ion Proton System (Thermo Fisher Scientific).
Raw poly(A) RNA-Seq reads were processed with Torrent Suite Software v5.10.1, which runs the Torrent Mapping Alignment Program (TMAP) using the L. decemlineata assembly genome version Ldec 2.0 from NCBI [23]. The number of genes in the L. decemlineata genome predicted in this assembly genome totaled 24,671 gene transcripts and 93,782 predicted exons. Reads were then submitted to a two-step mapping process. The first step consisted of mapping using STAR v.2.5.2b [24], while the unmapped reads were subsequently mapped using Bowtie 2 v2.2.5 [25]. The two aligned files were merged using Picard tools (release 2.5) [26]. Raw counts were extracted with featureCounts v1.4.6-p4 [27]. Gene expression in fragments per kilobase of transcript per million mapped reads (FPKM) values was subsequently calculated with Cufflinks v2.2.1 [28], an approach supported for single-end reads assessment on this sequencing platform and exemplified elsewhere [29,30]. To test for statistical differences in differentially expressed genes between the control and treatment conditions, unpaired Student's t-tests were performed in the R statistical environment v3.4.1 [31] after removing genes exhibiting low expression (FPKM < 1 in two samples for both conditions).

Small RNA Library Construction and Sequencing
Small RNA fractions quality (n = 4), generated with two insects per replicate, was assessed on an Agilent TapeStation using a standard RNA ScreenTape assay. Small RNA libraries were prepared from 250 ng of enriched small RNA using the RNA-Seq v2 kit (Thermo Fisher Scientific). Small RNA transcripts were ligated to Ion Torrent adapters followed by reverse transcription (Ion Total RNA-seq v2 kit, Thermo Fisher Scientific). Resulting cDNA products were next subjected to double size selection. This step removed the majority of larger cDNA fragments and enriched the cDNA fragments from small RNA transcripts. Barcoded libraries were generated from the enriched cDNA fragments using Ion Xpress RNA barcodes (Thermo Fisher Scientific). Ion Torrent recommended conditions were followed throughout the protocol. Quality of the libraries was analyzed using a TaqMan library quantification assay (Ion Library TaqMan Quantitation Kit, Thermo Fisher Scientific). Barcoded libraries were equally pooled and amplified onto ISPs supplied with the Ion PI Hi-Q OT2 kit (Thermo Fisher Scientific). ISPs enriched with template libraries were loaded onto Ion PI chip v3 and sequenced as above.
Small RNA-Seq raw reads were processed with Torrent Suite Software v5.10.1, which runs the Torrent Mapping Alignment Program (TMAP) using the L. decemlineata genome mentioned above. Removal of the adapters and selection of reads 15 to 30 bp in length were performed using Cutadapt v1.8.1 [32] and Picard tools (release 2.5) [26]. Reads were mapped to miRNA sequences of L. decemlineata obtained from the Ldec 2.0 genome [23] using Bowtie 2 v2.2.5 [25], and read counts were obtained with featureCounts v1.4.6-p4 [27]. R statistical environment v3.4.1 and Bioconductor package edgeR v3.18.1 [33] were used to identify differentially expressed miRNAs. Transcripts with fewer than 10 reads in more than half the samples under both conditions were removed. Analysis was performed using the trimmed mean of M-values (TMM) normalization method [34]. Transcriptomics data sets relevant to the current study are accessible through the Sequence Read Archive (SRA) of NCBI under the bioproject identification PRJNA635016.

Synthesis of cDNA
First strand synthesis for mRNA amplification was performed by mixing 1 µg of total RNA-which was obtained from large RNA isolates that were previously RNA isolates (3) used for NGS work or new isolates (2) not used as part of NGS analysis-with 1 µL oligo dT and 1 µL 10 mM dNTPs. Diethyl pyrocarbonate (DEPC)-treated water was next added to a volume of 12 µL. The solution was incubated at 65 • C for 5 min. These reagents were subsequently added to the mixture: 4 µL 5× First Strand Buffer, 2 µL 0.1 M DTT and 1.5 µL DEPC-treated water. This solution was incubated at 37 • C for 2 min, then 0.5 µL of M-MLV Reverse Transcriptase (RT) was added. Final incubation steps at 37 • C for 50 min and 70 • C for 15 min were ultimately performed.

qRT-PCR Amplification of Spinosad-Associated Transcripts Identified via NGS
Primers for amplification and quantification of transcripts identified via NGS were conceived based on L. decemlineata sequences for each transcript. Primers are presented in Table 1. Initial reactions contained 5 µL of diluted cDNA template (10 −1 ), 1 µL 25 µM forward primer, 1 µL 25 µM reverse primer, 5.5 µL DEPC-treated water and 12.5 µL 2X Taq FroggaMix (FroggaBio, Concord, ON, Canada). Amplification protocol consisted of a denaturing step at 95 • C for 5 min, followed by 40 Insects 2020, 11, 820 5 of 16 cycles at 95 • C for 15 s, at a temperature gradient between 50 and 65 • C for 30 s and at 72 • C for 45 s. Amplified products were resolved on a 2% agarose gel and sequenced at the Université Laval sequencing platform (Quebec City, QC, Canada). Targets of interest were also amplified via qRT-PCR at different annealing temperatures and on serial cDNA dilutions to measure primer set efficiencies. Reactions conducted to quantify transcript levels consisted of 2.5 µL of cDNA template (10 −1 ), 0.5 µL DEPC-treated water, 1 µL 5 µM forward primer, 1 µL 5 µM reverse primer and 5 µL of iTaq Universal SYBR Green Supermix (Bio-Rad, Hercules, CA, USA). The amplification protocol comprised an initial step at 95 • C for 3 min, followed by 40 cycles at 95 • C for 15 s and at the identified optimal annealing temperature for 30 s. Transcript levels of α-tubulin were used as reference and were amplified as above using 55 • C as the annealing temperature. Significant differences between transcript target expression measured in the control versus spinosad-exposed insects were assessed with an unpaired Student's t-test. Relative normalized transcript expression and a statistical analysis of the data were conducted with CFX Maestro software v1.1 (Bio-Rad). Table 1. Primers for the qRT-PCR-based quantification of spinosad-associated transcripts, as measured by next-generation sequencing (NGS). Abbreviations: Forward (Fwd) and reverse (Rev) primers.

Primer
Sequence Efficiency Temperature

Functional Assessment of Spinosad-Modulated mRNAs
Functional annotation of the most differentially expressed mRNA transcripts in spinosad-exposed L. decemlineata, as determined by next-generation sequencing was performed using Blast2GO Basic v5.2.5 [35]. Up-and down-regulated L. decemlineata transcript gene sequences were mapped using the Blast Diptera database (taxa: 7147) to identify homologous dipteran transcripts to use for functional classification with Blast2GO. Dipteran genes obtained with e-values below the cut-off of 1.0 E-3 were used for subsequent functional annotation.

mRNA Expression in Spinosad-Exposed L. decemlineata by NGS
Next-generation sequencing revealed 182,009,095 reads among all samples. Filtering was performed to discard low-quality reads and generated an average of 27,398,426 and 21,001,758 mapped reads for control and spinosad-treated insects, respectively. Reads were mapped to an average of 13,281 genes among all samples. Transcripts with absolute log2 fold-changes greater than 1.0 and p values below 0.05 were further identified. Up-regulation in spinosad-treated L. decemlineata was measured for 21 transcripts, and down-regulation was monitored for nine transcripts. Results are presented in Figure 1 and Table 2. Transcripts showing the greatest up-regulation in insects treated with spinosad as determined by next-generation sequencing were peritrophin-1-like (LOC111514273), polygalacturonase/glycoside hydrolase family 28 (LOC111505285), cathepsin B-like (LOC111505933), lysosomal alpha-mannosidase-like (LOC111509352) and venom carboxylesterase-6-like (LOC111507304).
Transcripts exhibiting the most reduced levels in spinosad-exposed insects when compared with control insects were cathepsin L proteinase (LOC111513727), voltage-dependent calcium channel type A subunit alpha-1-like (LOC111507833), putative nuclease HARBI1 (LOC111517073), senecionine N-oxygenase-like (LOC111507616) and lipase 3-like (LOC111513879). average of 13,281 genes among all samples. Transcripts with absolute log2 fold-changes greater than 1.0 and p values below 0.05 were further identified. Up-regulation in spinosad-treated L. decemlineata was measured for 21 transcripts, and down-regulation was monitored for nine transcripts. Results are presented in Figure 1 and Table 2. Transcripts showing the greatest upregulation in insects treated with spinosad as determined by next-generation sequencing were peritrophin-1-like (LOC111514273), polygalacturonase/glycoside hydrolase family 28 (LOC111505285), cathepsin B-like (LOC111505933), lysosomal alpha-mannosidase-like (LOC111509352) and venom carboxylesterase-6-like (LOC111507304). Transcripts exhibiting the most reduced levels in spinosad-exposed insects when compared with control insects were cathepsin L proteinase (LOC111513727), voltage-dependent calcium channel type A subunit alpha-1-like (LOC111507833), putative nuclease HARBI1 (LOC111517073), senecionine N-oxygenase-like (LOC111507616) and lipase 3-like (LOC111513879).   Transcripts coding for targets with proposed relevance to insecticide responses in L. decemlineata, including cytochrome P450s [36] and cuticular proteins [37], were also identified as above using less stringent cut-offs. Transcripts displaying absolute log2 fold-changes above 0.25 and p values below 0.1 were obtained from NGS data. Results are shown in Table 3. Transcripts showing the greatest up-regulation in spinosad-exposed insects versus control insects consisted of the ones coding for cytochrome P450 9e2 (LOC111506250) and cytochrome P450 6k1 (LOC111510879), while transcripts exhibiting the strongest down-regulation under the same conditions consisted of a mitochondrial cytochrome P450 301a1 (LOC111516463) and cytochrome P450 4c1 (LOC111509981). Table 3. Differential expression of select transcripts with putative relevance to insecticide responses in L. decemlineata treated with spinosad as measured by high-throughput sequencing.

qRT-PCR-Based Expression of Transcripts Modulated by Spinosad in L. decemlineata as Determined by NGS
Expression of six targets with varying levels following NGS was measured by qRT-PCR in control and insecticide-exposed insects. Results are shown in Figure 2. Transcript levels of cathepsin B (LOC111505933), lysosomal alpha-mannosidase (LOC111509352) and multidrug resistance-associated protein 4 (LOC111503007) were elevated in spinosad-treated insects to relative normalized expression values that were 1.62-, 1.95-and 1.76-fold higher than the ones measured in control insects, respectively. Transcript levels of cathepsin L proteinase (LOC111513727), senecionine N-oxygenase (LOC111507616) and voltage-dependent calcium channel type A subunit alpha-1 (LOC111507833) displayed reduced expression in spinosad-exposed insects to relative normalized expression values that were 0.09-, 0.32and 0.52-fold higher than the ones observed in untreated insects, respectively. Changes observed by qRT-PCR displayed similar trends in expression to the ones measured by NGS for all transcripts, although most reported changes were not statistically significant.

qRT-PCR-Based Expression of Transcripts Modulated by Spinosad in L. decemlineata as Determined by NGS
Expression of six targets with varying levels following NGS was measured by qRT-PCR in control and insecticide-exposed insects. Results are shown in Figure 2. Transcript levels of cathepsin B (LOC111505933), lysosomal alpha-mannosidase (LOC111509352) and multidrug resistanceassociated protein 4 (LOC111503007) were elevated in spinosad-treated insects to relative normalized expression values that were 1.62-, 1.95-and 1.76-fold higher than the ones measured in control insects, respectively. Transcript levels of cathepsin L proteinase (LOC111513727), senecionine N-oxygenase (LOC111507616) and voltage-dependent calcium channel type A subunit alpha-1 (LOC111507833) displayed reduced expression in spinosad-exposed insects to relative normalized expression values that were 0.09-, 0.32-and 0.52-fold higher than the ones observed in untreated insects, respectively. Changes observed by qRT-PCR displayed similar trends in expression to the ones measured by NGS for all transcripts, although most reported changes were not statistically significant.

miRNA Expression in Control and Spinosad-Exposed Insects via NGS
High-throughput sequencing revealed miRNAs with elevated levels in L. decemlineata, including miRNA miR-8-3p (82,033 reads in control and 48,570 reads in spinosad-treated insects), miR-14-5p (57,408 reads in control and 69,918 reads in spinosad-treated insects) and miR-317-5p (49,768 reads in control and 34,526 reads in spinosad-treated insects). A list of the most strongly expressed miRNAs in L. decemlineata is shown in Table 4. Furthermore, miRNA transcripts with absolute log2 fold-changes above 0.5 and p values below 0.05 were identified. Ratios of miRNA expression in spinosad-treated versus control L. decemlineata revealed differential expression of 14 miRNAs. A total of nine up-regulated miRNAs, including miR-750-3p, miR-2796-5p and miR-3791-5p, as well as five down-regulated, such as miR-8-3p, miR-9a-5p and miR-12-3p, were identified. MiRNAs modulated in spinosad-exposed insects are presented in Table 5.

Functional Classification of Differentially Expressed mRNA Transcripts Following Spinosad Treatment
Functional annotation of spinosad-modulated transcripts as determined by next-generation sequencing was performed. Biological processes were obtained for select transcripts that displayed substantial increased or decreased expression in L. decemlineata exposed to spinosad. Identified transcripts were associated, using Blast2GO, with GO terms that included primary metabolic process (GO:0044238), proteolysis (GO:0006508), nucleoside transmembrane transport (GO:1901642), response to hypoxia (GO:0001666), oxidation-reduction process (GO:0055114), lipid metabolic process (GO:0006629) and ion transport (GO:0006811), to name a few. A complete list of the processes identified is presented in Table 6. Table 6. Processes associated with select spinosad-modulated mRNA transcripts in L. decemlineata. L. decemlineata accession numbers and corresponding targets obtained from the dipteran database are shown.

Discussion
Successful management of insecticide resistance in insect pests that threaten agricultural production requires a better understanding of the transcriptional changes that occur in response to exposure to a given chemical. Multiple reports have leveraged transcriptomics-based approaches in order to highlight the molecular changes observed in L. decemlineata exposed to diverse compounds including chlorothalonil and imidacloprid [38]. However, similar studies are lacking with respect to spinosad response in L. decemlineata, despite its increased usage and need by organic producers. This work thus aimed to better characterize the modulated transcripts underlying spinosad response in this insect pest.
Multiple transcripts were identified in this study as differentially expressed in spinosad-exposed Colorado potato beetles when compared with control insects using high-throughput sequencing. Varying levels of select transcripts with functions linked to insecticide resistance were observed. For example, elevated transcript levels of venom carboxylesterase-6 (LOC111507304), a carboxylesterase (CBE), were measured in spinosad-treated L. decemlineata. Multiple reports have showcased the involvement of CBEs in the resistance of insects towards different chemicals including organophosphates, carbamates and pyrethroids [39]. Assessment of CBE activity in the diamondback moth P. xylostella demonstrated a positive correlation with several insecticides including spinosad [40]. Members of a CBE and cholinesterase superfamily have also been previously investigated in L. decemlineata, and it was concluded that their expression could be influenced following insect treatment with fipronil or cyhalotrin [41]. Interestingly, mutations in a carboxylesterase gene of the house fly Musca domestica (Diptera: Muscidae) are potentially correlated with spinosad resistance [42]. Besides CBE, up-regulation of a transcript coding for multidrug resistance-associated protein 4 (LOC111503007) was also observed following spinosad treatment in the current work. Modulation of multidrug resistance-associated proteins with the primary function of transporting a broad array of compounds (e.g., xenobiotics) across cellular membranes [43] has been observed in several insects exposed to insecticides. Exposition of the bird cherry-oat aphid Rhopalosiphum padi (Hemiptera: Aphididae) to imidacloprid or chlorpyrifos was associated with elevated transcript levels of ABCC1, coding for the multidrug resistance-associated protein 1 (MRP1) [44]. Furthermore, the same study reported higher levels of this target in insect strains that were resistant to these chemicals. It is interesting to note that several genes coding for ATP-binding cassette transporters displayed increased transcript levels in imidacloprid-resistant L. decemlineata [45]. Furthermore, significant elevation of ATP-binding cassette subfamily G transporter transcript levels was observed in a commercially managed population of L. decemlineata that was imidacloprid-resistant when compared with an organically managed susceptible counterpart [46], which also supports the relevance of such transporters in insecticide resistance with respect to this insect pest. Functional annotation performed in the current work highlighted processes associated with multidrug resistance-associated protein 4, including response to oxidative stress. This aligns with previous work showing that spinosad can induce oxidative stress in various models including Sf9 cells [47] and the Nile tilapia Oreochromis niloticus [48]. Future work is nevertheless warranted to build on the functional insights obtained here regarding the identified transcripts with potential relevance to spinosad. Subsequent in-depth functional enrichment analysis and a closer investigation of possible gene orthologues existing in other species with better defined functions would suit this objective. Overall, the data presented here reveal that spinosad treatment in Colorado potato beetles was linked with differential expression of multiple transcripts associated with insecticide responses.
Spinosad exposure also led to differential expression in L. decemlineata of transcripts underlying insect development. Chitinase 10 (LOC111512533) was up-regulated in insects treated with this compound. Chitinases can degrade the polysaccharide chitin, a key exoskeletal component in insects, and their underlying activity is notably important in insect development and molting [49]. Previous studies have explored the potential contribution of chitinases towards insecticide responses. For example, treatment of mosquito Culex pipiens (Diptera: Culicidae) larvae with the chemical lambda-cyhalotrin has been associated with elevated chitinase activity [50], while treatment of D. melanogaster with the tetranortriterpenoid insecticide azadirachtin has been associated with reduction in chitinase activity [51]. It is interesting to note that in our study, peritrophin-1 (LOC111514273) was strongly up-regulated in L. decemlineata treated with spinosad. Peritrophins are chitin-binding proteins and key components of the peritrophic membrane in the insect gut [52]. Reduction of chitin and peritrophins using an approach based on double-stranded RNA (dsRNA) has resulted in increased sensitivity to imidacloprid in the eastern subterranean termite Reticulitermes flavipes (Isoptera: Rhinotermitidae) [53]. Besides chitinase 10 and peritrophin-1, elevated transcript levels of juvenile hormone esterase (LOC111512328) were also detected in this study. Exposure of P. xylostella to pyriproxyfen has been correlated with elevated transcript levels and the activity of juvenile hormone esterase [54]. In addition, treatment of diamondback moth larvae with 3-octylthio-1,1,1-trifluoropropan-2-one (OTFP), a trifluoromethyl ketone that displayed inhibitory activity against juvenile hormone esterase, was associated with differential responses to various chemicals such as diafenthiuron and indoxacarb [55]. Collectively, these changes highlight key modulation in Colorado potato beetles of select targets involved in insect development following spinosad exposure.
Data presented here also showed differential expression of several miRNAs in Colorado potato beetles exposed to spinosad when compared with control insects. Since the pioneering report that identified 460 putative miRNAs in L. decemlineata [56], several studies have reported miRNA modulation in Colorado potato beetles exposed to different insecticides. A next-generation sequencing-based approach revealed the modulation of 33 miRNAs, such as miR-100 and miR-989, in L. decemlineata treated with imidacloprid [57]. The former was notably shown to target the cytochrome P450 CYP6CY3, which had been linked with neonicotinoid resistance [58,59]. Numerous specific differentially expressed cytochrome P450s, part of the CYP450 gene superfamily from which multiple members have been investigated for their underlying role in insecticide resistance in L. decemlineata [60,61], were identified via high-throughput sequencing in the current work, supporting further investigation of the miRNA-mediated regulation of cytochrome P450s following spinosad exposure. A different signature of miRNAs displaying up-regulation of miR-1-3p and miR-305-5p was observed in L. decemlineata exposed to chlorantraniliprole [22]. Recent work also reported the potential involvement of miR-9a-5p and miR-965-5p in hexanoic acid response in L. decemlineata [62]. The current study highlighted a signature of different miRNAs that were responsive to spinosad, including miR-9a-5p, miR-12-3p and miR-750-3p, which further highlights the variations in miRNA footprints elicited by specific insecticides in L. decemlineata. Subsequent work aimed at better clarifying the targets regulated by these miRNAs is warranted to assess the potential relevance of these leads in spinosad response in this insect pest. It is also worth pointing out that follow-up investigations using longer spinosad exposures would contribute to a more in-depth characterization of molecular changes, both shortand mid-term, associated with spinosad response in L. decemlineata. Treatment protocols using higher doses of spinosad, combined with increased stringency of the thresholds used for NGS data analysis, could also provide a more thorough and rigorous overview of the transcriptional variations linked with response to this chemical and build on the data reported in this study given the limited number of transcripts presented.

Conclusions
Given the accumulating reports of spinosad resistance in select populations of Colorado potato beetles [63] and the long-term objective of better characterizing the changes in gene expression potentially affected by spinosad, the current study investigated the molecular changes linked to spinosad response in Colorado potato beetles using a high-throughput sequencing approach. A signature of mRNAs coding for diverse targets with relevance to spinosad response were identified, including select chitinases, multidrug resistance-associated protein and juvenile hormone esterase. Expression levels of select transcripts identified using the NGS approach were subsequently quantified via qRT-PCR in an expanded number of replicates. Comparable trends in expression status, albeit not all statistically significant, were observed following NGS-versus qRT-PCR-based quantification. In addition, spinosad-associated differential expression of miRNAs, including miR-12-3p and miR-750-3p, was also observed for the first time in L. decemlineata exposed to spinosad. Overall, this work showed signatures of transcripts that responded to spinosad in L. decemlineata and further revealed potential targets to consider when designing novel approaches aimed at limiting the damages this insect can cause to potato fields.
Supplementary Materials: The following are available online at http://www.mdpi.com/2075-4450/11/11/820/s1, Table S1: Impact of spinosad on L. decemlineata. Mortality/impairment rates of insects exposed to varying doses of spinosad are shown. is acknowledged. Salary support for PD was covered in part via the New Brunswick Innovation Foundation (NBIF) Research Assistantships Initiative (RAI2018-062). Salary support for PB was covered in part via a Canadian Horticultural Council-Agri-Science Cluster for Horticulture 3 grant.