Identification of Differentially Expressed miRNAs in Colorado Potato Beetles (Leptinotarsa decemlineata (Say)) Exposed to Imidacloprid

The Colorado potato beetle (Leptinotarsa decemlineata (Say)) is a significant pest of potato plants that has been controlled for more than two decades by neonicotinoid imidacloprid. L. decemlineata can develop resistance to this agent even though the molecular mechanisms underlying this resistance are not well characterized. MicroRNAs (miRNAs) are short ribonucleic acids that have been linked to response to various insecticides in several insect models. Unfortunately, the information is lacking regarding differentially expressed miRNAs following imidacloprid treatment in L. decemlineata. In this study, next-generation sequencing and quantitative real-time polymerase chain reaction (qRT-PCR) were used to identify modulated miRNAs in imidacloprid-treated versus untreated L. decemlineata. This approach identified 33 differentially expressed miRNAs between the two experimental conditions. Of interest, miR-282 and miR-989, miRNAs previously shown to be modulated by imidacloprid in other insects, and miR-100, a miRNA associated with regulation of cytochrome P450 expression, were significantly modulated in imidacloprid-treated beetles. Overall, this work presents the first report of a miRNA signature associated with imidacloprid exposure in L. decemlineata using a high-throughput approach. It also reveals interesting miRNA candidates that potentially underly imidacloprid response in this insect pest.


Introduction
The Colorado potato beetle (CPB) (Leptinotarsa decemlineata (Say)) is a significant insect pest harming potato crops worldwide [1]. CPBs are often considered to be the primary insect associated with potato plant defoliation [2], leading to up to 75% of foliage consumption and impacting revenues for growers [3,4]. Pest control strategies targeting CPBs typically involve pesticides, even though resistance against a variety of such compounds, including spinosad, thiamethoxam, and deltamethrin, has been observed [5][6][7]. Neonicotinoids, which target the nicotinic acetylcholine receptors (nAChRs), are a class of insecticides that are used extensively against CPBs. The neonicotinoid imidacloprid (1-(6-chloro-3-pyridylmethyl)-N-nitroimidazolidin-2-ylideneamine) was registered for such a purpose over 20 years ago, and was initially successful in managing CPB populations [6]. Unfortunately, Int. J. Mol. Sci. 2017, 18, 2728 2 of 12 early studies reported adult CPBs exhibiting as much as a 155-fold increase in imidacloprid resistance in select populations [8]. Efforts have been deployed in recent years to understand the molecular players associated with imidacloprid resistance in CPBs. Multiple studies have notably highlighted the differential expression of cytochrome P450s in imidacloprid-treated CPBs [9,10], positioning these enzymes as targets that could sensitize CBPs to this agent. Nevertheless, the complete molecular picture associated with imidacloprid resistance in CPBs remains incompletely characterized.
MicroRNAs (miRNAs) are short conserved non-coding RNAs capable of regulating the expression of multiple target mRNA transcripts. These molecules are involved in the response to various stresses including starvation, anoxia, and freezing [11,12]. Previous work has also revealed miRNA modulation in insects following exposure to different chemicals including pyrethroids, chlorantraniliprole and fenpropathrin [13][14][15]. In addition, recent studies have highlighted the likely regulation of pyrethroid resistance in the mosquito Culex pipiens pallens (L.) by miRNAs [16,17]. These results thus support the potential involvement of miRNAs in insecticide resistance. Unfortunately, data is currently lacking regarding a signature of differentially expressed miRNAs in response to imidacloprid treatment in CPBs.
The current study was undertaken to characterize the miRNA footprint observed in CPBs following imidacloprid exposure, and to identify novel miRNA-associated molecular leads that could be utilized in the development of alternative strategies to control CPBs using an approach relying on next-generation sequencing and qRT-PCR. The miRNA target prediction was also performed using predictive tools to better assess the likely consequences of these imidacloprid-associated miRNAs. Transcript targets linked to transcriptional regulation and glucose metabolism were observed and are further discussed.

Small RNA Sequence Analysis
High-throughput sequencing resulted in 59,406,387 reads among all samples. Low quality reads were discarded. Reads ranging from 16 to 60 nucleotides were conserved for analysis. Small RNA libraries displayed comparable read distribution profiles in the 16 and 24 nucleotides region, with maximum peaks observed at 21 nucleotides. This approach resulted in 22,534,030 and 13,763,160 total reads for control and imidacloprid-treated insects, respectively (Table 1). Within these reads, 8,825,610 (control) and 6,211,853 (imidacloprid) unique reads were obtained. MiRBase was utilized to identify the small RNAs of interest in L. decemlineata samples using red flour beetle Tribolium castaneum reference databases [18]. The miRNAs were annotated using sRNAbench in sRNAtoolbox [19]. Small RNAs unique reads were mapped to miRNAs (26,447 for control and 24,221 for imidacloprid-treated). SnRNAs (4735 for control and 3616 for imidacloprid-treated), snRNAs (4705 for control and 2937 for imidacloprid-treated) and tRNAs (40,144 for control and 22,115 for imidacloprid-treated) were also observed. The most frequent miRNAs detected with this approach were miR-14-3p (95,802 reads in control and 81,633 reads in imidacloprid-exposed insects) and miR-8-3p (68,314 reads in control and 79,134 reads in imidacloprid-exposed insects) ( Table 2). Based on next-generation sequencing results, log2 ratios of miRNA expression in control and imidacloprid-exposed L. decemlineata demonstrated 33 differentially expressed miRNAs. A total of 14 upregulated and 19 downregulated miRNAs displaying absolute log2 fold-changes greater than 0.3 (p < 0.05, n = 3) were identified (Table 3).

Discussion
The molecular changes underlying neonicotinoid resistance in insects have been investigated with great interest in recent years. Work conducted on various insect models treated with imidacloprid has revealed targets potentially involved in the response to this agent. These targets include select cytochrome P450s, ATP-binding cassette (ABC) transporters, and nAChRs [20][21][22]. Numerous examples exist that highlight the regulation of these imidacloprid-relevant targets by miRNAs [23,24]. While recent studies have started to associate miRNAs with insecticide resistance in various insects [14,25,26], no work so far has attempted to characterize modulated miRNAs in imidacloprid-treated CPBs. Using a high-throughput sequencing approach, the current study reveals a signature of differentially expressed miRNAs in CPBs exposed to this neonicotinoid.
The information is sparse regarding the molecular changes associated with imidacloprid response in insects. Pioneering work performed on beehives submitted to low levels of imidacloprid revealed a set of differentially expressed miRNAs in larvae gathered from imidacloprid-exposed versus unexposed hives [27]. The present study also identified similar changes in expression levels of select miRNAs, including miR-282 upregulation and miR-989 downregulation, in imidacloprid-treated CPBs. MiR-989 has been linked with key physiological processes, as Drosophila melanogaster (Meigen) miR-989 mutants notably exhibit impaired border cell migration [28], and Bombyx mori (L.) strongly upregulates miR-989 during the pupal phases [29]. It is interesting to note that levels of miR-989-3p, like miR-305-5p and miR-927a-5p, remained unchanged in insects treated with imidacloprid for different lengths of time, strengthening its potential involvement in response to this insecticide. MiR-282, on the other hand, can influence D. melanogaster viability and longevity by targeting the adenylate cyclase rutabaga [30]. While modulation of miR-282 and miR-989 in imidacloprid-treated CPBs presented here supports their probable role in response to this pesticide, additional work is required to fully delineate their functions in insects exposed to imidacloprid.
Several reports have highlighted the upregulation of cytochrome P450s in CPBs treated with imidacloprid. A high-throughput sequencing approach performed in imidacloprid-sensitive and imidacloprid-resistant CPBs revealed that transcripts for 41 cytochrome P450s were more than 2-fold upregulated in pesticide-exposed insects. In addition, cytochrome P450 activity was increased in the resistant beetles [10]. Gene expression analysis performed on CPB populations susceptible or resistant to imidacloprid reported differential cytochrome P450 levels. This further strengthens an underlying role for these enzymes in imidacloprid resistance [9]. Several studies have reported miRNA-mediated regulation of cytochrome P450 expression. Recent work conducted in the tobacco aphid Myzus persicae nicotianae (Blackman) showed that expression of CYP6CY3, a cytochrome P450 previously associated with neonicotinoid resistance [31], was regulated by miR-100 [32]. The present study highlighted downregulation of miR-100-5p levels in imidacloprid-treated CPBs, supporting a potential role for the miR-100-cytochrome P450 axis in the imidacloprid response.
Functional annotation was undertaken using the predicted targets of differentially expressed miRNAs following imidacloprid exposure in order to better characterize the potential impact of this compound in CPBs (Table 5). This approach notably highlighted gene transcription as a likely process affected by the modulated miRNAs. Some miRNAs, such as miR-92a-3p and miR-927a-5p, were predicted to impact key players involved in transcriptional regulation, including the cAMP response element binding (Creb) transcription factor and the ecdysone receptor (EcR). Previous work has reported that Creb levels were differentially expressed in honey bees (Apis mellifera carnica (Pollmann)) exposed to various neonicotinoids [33], supporting a potential role for Creb in the response to imidacloprid. Several studies have investigated the effect of imidacloprid treatment on gene expression in multiple insect models. A microarray-based approach performed in the mosquito Aedes aegypti (L.) revealed 344 and 108 differentially expressed genes in imidacloprid-resistant larvae and adults, respectively, when compared to their imidacloprid-sensitive counterparts [34]. In addition, a high-throughput sequencing approach in D. melanogaster that were sensitive or resistant to the same agent highlighted 357 transcripts that were either upregulated or downregulated in imidacloprid-resistant flies [35]. Multiple transcript targets associated with imidacloprid-modulated miRNAs were also linked to glucose metabolism. Previous reports have linked imidacloprid treatment with modulation of various targets involved in glucose homeostasis. Honey bee larvae exposed to imidacloprid exhibited increased transcript levels of multiple genes involved in the glycolytic and gluconeogenic pathways. Transcript levels of phosphoenolpyruvate carboxykinase (PEPCK), a rate-limiting enzyme for the latter cascade, were notably upregulated in imidacloprid-exposed larvae [27]. In addition, a study performed on Madagascar hissing cockroaches (Gromphadorhina portentosa (Schaum)) showed variations in glucose uptake and carbohydrate metabolism in cockroaches exposed to imidacloprid [36]. These leads generate potential processes via which the identified miRNAs can influence imidacloprid response. However, further characterization of the link that exists between these processes and imidacloprid-associated miRNAs is envisioned. Ultimately, such transcripts could be targeted by diverse means including RNA interference (RNAi)-based approaches that form part of a strategy aimed at influencing imidacloprid response in CPBs. Such approaches are currently under investigation in various pests including CPBs [37].
In conclusion, the present study revealed a set of modulated miRNAs in CPBs following exposure to the neonicotinoid imidacloprid, and is aligned with recent work aimed at better understanding the functions associated with miRNAs in these insects [38]. These results present the first example of a miRNA signature in L. decemlineata treated with this agent. Assessment of miRNA expression following different imidacloprid treatment durations also revealed varying levels of select miRNAs in short-versus long-term insecticide exposures. This warrants further investigation of early and late molecular changes associated with response to this agent. Predicted miRNA transcript targets further suggest that differential miRNA expression could impact molecular players involved in key processes. It is important to note that while the miRNA changes uncovered in this study generated valuable information regarding miRNAs underlying insect response to imidacloprid, future work is needed to investigate miRNA signatures observed in response to additional experimental conditions. These conditions include repeated topical imidacloprid exposures, usage of a diet comprised of imidacloprid-treated potato leaves, or field CPBs naturally exposed to this agent and exhibiting various degrees of tolerance against imidacloprid. Further investigations will be undertaken to examine the expression status of enzymes such as cytochrome P450s in CPB populations that do or do not respond to imidacloprid. Overall, this work provides novel information on a potential role for miRNAs in imidacloprid response and resistance, as well as contributing to the growing knowledge of the molecular levers underlying the response to this agent in Colorado potato beetles.

Insect Collection and Treatment
Adult Colorado potato beetles that had overwintered were collected in nursery potato fields at the Fredericton Research and Development Centre in Fredericton (NB, Canada) in June 2016. Potato fields where collection occurred were not treated with any types of insecticides. This nevertheless cannot rule out the potential exposure of beetles to insecticides in years prior to the study, nor exposure attributable to insect migration. Beetles were placed in plastic containers (50 per) containing potato leaves, and closed using a lid with insect screening for ventilation. Insects were subsequently transported to the Université de Moncton (Moncton, NB, Canada). A group of 30 insects was placed in an incubator (Thermo Fisher Scientific, Waltham, MA, USA) set at 25 • C for 5 days under 16L:8D cycles. Insects were provided with potato plants. Following acclimation, 5 µL (0.5 µg) of analytical grade imidacloprid (100 µg/mL in acetonitrile, Sigma-Aldrich, St. Louis, MO, USA) was applied topically on the abdomen of 15 beetles, as previously described [10]. A volume of 5 µL acetonitrile was applied on 15 beetles in parallel, which were used as controls. Imidacloprid-treated and control insects were returned to the incubator for 8 h. A similar experimental protocol was undertaken where insects were treated with imidacloprid for 24 h. Insects were active after the incubation period and prior to storage. A parallel bioassay performed on insects treated with increasing imidacloprid doses (0.05 µg to 5.0 µg) and not destined for high-throughput sequencing highlighted a marked impact on insect activity, albeit not viability, at maximum doses of 2.5 µg and 5.0 µg of imidacloprid (LD 50 >5.0 µg/insect). All insects destined for sequencing were ultimately placed rapidly in liquid nitrogen and stored at −80 • C until RNA isolation.

Small RNA Isolation
Small RNA fractions were isolated from control and imidacloprid-treated L. decemlineata using the miRVana miRNA Isolation Kit (Thermo Fisher Scientific, Waltham, MA, USA) as per the manufacturer's instructions and as described previously [39]. Each isolate was prepared using two L. decemlineata. RNA isolates were generated in triplicates. The quality and integrity of small RNA fractions were confirmed using the High Sensitivity RNA ScreenTape on an Agilent TapeStation 2200 instrument (Agilent Technologies (Santa Clara, CA, USA)).

Small RNA Library Construction and Sequencing
Small RNA libraries were prepared and sequenced following Ion Torrent protocols and as reported previously [40]. Samples were loaded on an Ion PI Chip v2 and sequenced with an Ion Proton Sequencer (Thermo Fisher Scientific, Waltham, MA, USA). Three control and three imidacloprid-treated libraries were constructed. FASTQ files were generated for each sample and adaptor sequences were trimmed.
Reads with less than 16 and more than 60 nucleotides were discarded. Reads with Q scores of less than 20 were not included. Small RNA annotation was performed using the sRNAbench tool from sRNAtoolbox. Red flour beetle Tribolium castaneum was used as a reference [19]. The number of permitted mismatches when mapping to mature miRNAs was set to 3. Default values were used for the remaining parameters. Annotated sequences having less than 10 normalized read counts in either control or insecticide-treated samples were discarded, and log2 fold-changes were generated.

PCR and qRT-PCR Amplification of miRNAs
Initial PCR reactions were performed to confirm correct product amplification. Reactions consisted of 5 µL of cDNA (10 −1 ), 5.5 µL of DEPC-treated water, 1 µL of 25 µM miRNA-specific forward primer, 1 µL of 25 µM universal reverse primer and 12.5 µL of 2× Taq FroggaMix (FroggaBio, Toronto, ON, Canada) [42]. Forward and universal primers used are presented in Table 6. An initial denaturing step at 95 • C for 5 min was performed, followed by 35 cycles at 95 • C for 15 s and at a gradient of temperatures for 1 min. Products were separated on a 2% agarose gel and visualized using a ChemiDoc MP system (Bio-Rad, Hercules, CA, USA). Products were then sequenced at the Université Laval (Quebec City, QC, Canada) and identities were confirmed using Basic Local Alignment Search Tool (BLAST).
The qRT-PCR reactions were performed by mixing 2.5 µL of diluted cDNA template with 0.5 µL of DEPC-treated water, 1 µL of 5 µM miRNA-specific forward primer, 1 µL of 5 µM universal reverse primer and 5 µL of 2× iTaq Universal SYBR Green Supermix (Bio-Rad, Hercules, CA, USA) prepared in triplicate in 96-well plates. Amplification protocol consisted of an initial denaturing step at 95 • C for 3 min, followed by 40 cycles at 95 • C for 15 s and optimal annealing temperature for each miRNAs (Table 6) for 30 s. Reactions were conducted on a CFX Connect Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA). Efficiencies were calculated for each primer pair with the same protocol as above in serial cDNA dilutions. NormFinder software (Aarhus University Hospital, Aarhus, Denmark) [43] revealed miR-1-3p as the most stably expressed miRNA amongst a group of transcripts in insects exposed to imidacloprid for 8 h as determined by qRT-PCR, and was used as reference transcript.

miRNA Transcript Targets Prediction
TargetScanFly 6.2 [44] and the fruit fly miRanda algorithm [45] target prediction tools were used to identify likely mRNA candidates regulated by the imidacloprid-associated miRNAs identified by high-throughput sequencing. The top five transcript targets linked with imidacloprid-modulated miRNAs that were available in both prediction tools were obtained. A list of miRNA targets was subsequently generated and functionally annotated with the Database for Annotation, Visualization and Integrated Discovery (DAVID) bioinformatics resources v6.8 [46]. Biological processes associated with these targets were generated.

Quantification and Statistics
The Cq values following qRT-PCR runs were obtained using the Bio-Rad CFX Manager software (Bio-Rad, Hercules, CA, USA). miRNA levels were normalized using miR-1-3p levels amplified from the same sample. The miRNA levels were measured using the 2 −∆∆Cq method [47]. Ratios of normalized miRNA levels in control CPBs to average transcript expression in imidacloprid-treated CPBs were generated. Statistical differences between control and treated conditions were assessed with the Student's t-test. For high-throughput sequencing, annotated sequences were normalized using the trimmed mean of M-values (TMM) function in edgeR Bioconductor [48], and log2 ratios depicting miRNA expression in insecticide-treated versus control CPBs were calculated. Bioconductor "limma" package, with linear model fit and empirical Bayes statistics, was used to calculate differential expression [49]. Results for miRNAs with average normalized expression (ANE) greater than 10, log2 fold-change higher than 0.3 and p < 0.05 were investigated.