Differentially expressed microRNAs associated with changes of transcript levels in detoxification pathways and DDT-resistance in the Drosophila melanogaster strain 91-R

Dichloro-diphenyl-trichloroethane (DDT) resistance among arthropod species is a model for understanding the molecular adaptations in response to insecticide exposures. Previous studies reported that DDT resistance may involve one or multiple detoxification genes, such as cytochrome P450 monooxygenases (P450s), glutathione S-transferases (GSTs), esterases, and ATP binding cassette (ABC) transporters, or changes in the voltage-sensitive sodium channel. However, the possible involvement of microRNAs (miRNAs) in the post-transcriptional regulation of genes associated with DDT resistance in the Drosophila melanogaster strain 91-R remains poorly understood. In this study, the majority of the resulting miRNAs discovered in small RNA libraries from 91-R and the susceptible control strain, 91-C, ranged from 16–25 nt, and contained 163 precursors and 256 mature forms of previously-known miRNAs along with 17 putative novel miRNAs. Quantitative analyses predicted the differential expression of ten miRNAs between 91-R and 91-C, and, based on Gene Ontology and pathway analysis, these ten miRNAs putatively target transcripts encoding proteins involved in detoxification mechanisms. RT-qPCR validated an inverse correlation between levels of differentially-expressed miRNAs and their putatively targeted transcripts, which implies a role of these miRNAs in the differential regulation of detoxification pathways in 91-R compared to 91-C. This study provides evidence associating the differential expression of miRNAs in response to multigenerational DDT selection in Drosophila melanogaster and provides important clues for understanding the possible roles of miRNAs in mediating insecticide resistance traits.


Introduction
While chemical insecticidal agents have been developed and widely applied to suppress pest arthropod populations in ongoing efforts to enhance the efficiency of agricultural production PLOS ONE | https://doi.org/10.1371/journal.pone.0196518 April 26, 2018 1 / 21 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 potent regulators of eukaryotic gene expression at the post-transcriptional level [36,37] that regulate cell differentiation [38], migration [39], and neuronal development [40]. Approximately 256 D. melanogaster miRNA precursors have been identified and deposited in the mir-Base (Release 21) [24] and have been shown to be involved in regulating biological processes such as development [41], immune response [42], and metabolism [43]. Recent publications implicate a role of specific miRNAs in the regulation of insecticide resistance mechanisms among insect species. Several comparative analyses have estimated significant differences in miRNA expression between chemical insecticide-resistant and -susceptible strains [44,45]. Additionally, miRNAs differentially expressed between insecticidal Bacillus thuringiensis (Bt) toxin-resistant and -susceptible strains of Ostrinia furnacalis were predicted to target potential receptor genes [46]. In Culex pipiens, miRNAs differentially expressed between deltamethrin-sensitive and -resistant strains were proposed to mediate the expression of putative cytochrome P450 target genes [47]. Specifically, Culex pipiens miR-285 and miR-278 were implicated in pyrethroid resistance through the transcriptional regulation of Cyp6n23 and Cyp6ag11 [48,49], as well as a miRNA cluster involved in regulation of Cyp9j35 and Cyp325bg3 [50] and upregulation of miR-932 that regulates transcript levels of the cuticular gene CpCPR5 [51,52].
To date, the authors have no knowledge of any investigation of the contributions of posttranscriptional gene expression regulation on DDT resistance phenotypes. In order to partially address this knowledge gap, the present study estimated significant quantitative miRNA differences between the highly DDT-resistant D. melanogaster strain 91-R compared to the DDTsusceptible strain 91-C. Additionally, correlation between levels of differentially expressed miRNAs and corresponding putative computationally predicted P450 target transcripts were made within 91-R. This study provides insights into the role of miRNAs for the regulation of metabolic resistance to DDT as well as the effects of multigenerational DDT selection on the evolution of miRNA-mediated post-transcriptional regulation in a polygenic D. melanogaster DDT resistance phenotype.

miRNA library preparation, sequencing, and annotation
The DDT-resistant 91-R and -susceptible 91-C strains were obtained from Dr. Ranjan Ganguly (University of Tennessee-Knoxville) and developed as described previously [53]. Both strains were reared on a commercially available medium (Jazz-Mix Drosophila Food, Fischer Scientific, Cat. No. AS153) under the conditions of 25 ± 1˚C, 55-70% relative humidity and a 14 h light /10 h dark cycle. 91-R has been continually selected by maintaining the flies in colony bottles with the presence of a 150 mg DDT impregnated filter paper disk, while 91-C was maintained without any exposure to DDT. Recent topical bioassays estimated that 91-R is~107-fold more resistant to DDT compared to 91-C [17].
Three biological replicates of one hundred 3-5 day-old virgin female adults were collected from both 91-R and 91-C (n = 6). In order to compare the constitutive expression of miRNAs in subsequent analyses (see below), neither population was exposed to DDT within that generation. The small RNAs (sRNAs) were immediately extracted from live collected flies from each replicate using the Qiagen miRNeasy Mini Kit according to the manufacturer's instructions (Qiagen, Valencia, CA). RNA degradation and contamination were assessed for all samples using an Agilent 2100 Bioanalyzer (Agilent Technologies, Germany), and RNA concentrations were estimated using a NanoDrop One (Thermo, Wilmington, USA). Illumina sRNA libraries were constructed from each pool, and 50-bp single-end (SE50) sequence read data were generated on an Illumina HiSeq 4000 at the Research Technology Support Facility (Michigan State University, East Lancing, MI).
All raw Illumina sequence data were imported into CLC Genomics Workbench v.9.5 (Qiagen) and all reads were processed to remove low-quality reads, poly A sequences, adapters, reads without 3' adapter, and sequences shorter than 15 nt. Using the "small RNA analysis" tool in CLC Genomics Workbench, annotation of the trimmed read data from each library (n = 6) was made by comparing against the known miRNA precursors of D. melanogaster in the miR-Base R.21 (http://www.mirbase.org/; file hairpin.fa); subsequent generation of relative miRNA counts were made within each library. Only tags matching exactly with the mature 5' or 3' regions of previously annotated D. melanogaster miRNAs (miRbase R.21; file mature.fa) were accepted and retained for the further analysis. The sRNA sequencing data with annotated tags were deposited to NCBI Short Read Archive (SRA) with the accession number SRP136631. Additionally, miRDeep2 v.0.0.8 software was used in order to predict novel miRNA candidates [54]. Illumina adapters of sRNA raw sequence reads were trimmed using cutadapt v.1.15 [55]. The trimmed reads were quality-checked and curated using personal perl script. The results were aligned to the D. melanogaster reference genome (Release 6; dmel-all-chromosome-r6. 19. fasta.gz at http://flybase.org/) using bowtie v.1.2.2 [56] and also were mapped to D. melanogaster non-coding RNA database (dmel-all-tRNA-6.1/8.fasta and dmel-all-miscRNA-6.1/8.fasta file at flybase.org) in order to filter the small conditional RNAs (scRNAs), small nucleolar RNAs (snoRNAs), small nuclear RNAs (snRNAs) and transfer RNAs (tRNAs). Novel miRNAs were then identified against known miRNA precursors (hairpin.fa) and previously annotated miR-NAs (mature.fa) using default parameters suggested by the developers [57]. The structure and minimal free energy (MFE) of all potential novel miRNAs were predicted using a RNAfold [58] with algorithms described previously [59]. The MFE -25 kcal mol -1 , the randomization test of secondary structure MFE, called randfold P-value 0.05, and miRDeep2 score ! 3 were used as the cutoff level to declare them as potential novel miRNAs [60].

Differential expression of miRNAs
Estimates of miRNA expression level analysis within each replicate library from 91-R (n = 3) and 91-C (n = 3) were performed using the "Annotate and Merge" command on the CLC Genomics Workbench v.9.5 (Qiagen). This procedure used the D. melanogaster miRBase (release 21, http://www.mirbase.org/) as the primary database and the non-coding RNA database (dmel-all-tRNA-6.1/8.fasta and dmel-all-miscRNA-6.1/8.fasta) as the secondary database for annotation. The read counts of miRNAs were first normalized using the tag per million reads (TPM) method: TPM = (number of mapped reads for each miRNA/total number of mapped reads) ×10 6 . Subsequently, the Empirical analysis of Differential Gene Expression (EDGE) algorithm [61] was used to estimate differences in read counts comparing 91-R to 91-C strains, with P-values adjusted for multiple testing calculated using the Benjamini-Hochberg false discovery rate (FDR) procedure [62]. The variance in miRNA levels between the 91-R and 91-C with a log 2 fold-change > 1.0 or < -1.0, and a FDR 0.05, were defined as significant [63].

Target prediction and functional annotation of differentially expressed miRNAs
The potential target transcripts of miRNAs predicted to be differentially expressed between 91-R and 91-C were predicted using three different types of software packages, RNAhybrid [64], DIANA [65], and ComiR [66], using the following criteria: (1) RNAhybrid: the target site prediction was restricted to the 3'-UTR region obtained from the 3'-UTR database of D.

Validation and correlated expression between miRNAs and target transcripts
Reverse transcriptase-quantitative PCR (RT-qPCR) was carried out on selected eight known miRNAs for validation of sRNA sequencing estimated differences between the 91-R and 91-C strains. Moreover, the correlation was assessed between the relative expression levels of eight differentially expressed miRNAs and their potential targeted detoxification genes: Cyp6a8, Cyp6g1, Cyp6w1, Cyp4s3, Cyp6g2, Cyp309a2, Cyp313a4, Cyp313b1, Cyp4ae1, Cyp4d2, Cyp4g15, Cyp4p3, Cyp6d5, Cyp6t3, Cyp6u1, Cyp6v1, Cyp49a1, Cyp18a1, Cyp303a1, Cyp4aa1, Cyp4e3, Cyp4g1, Cyp6a19, Cyp6a22, Cyp6a9, Cyp4d20, GstD1, GstS1, GstE6, GstE10, Esterase10, Ester-ase7, Esterase8, ABC-B7, GABA-R, mAChR-A, GluR-IB, nAcRalpha-A, and Cpr65Ec. Three biological replicates of 15 adult female flies were sampled per strain (91-R and 91-C) andidentically with treatments used in Illumina sRNA sequencing-all samples were not exposed to DDT within the generation used. Both miRNAs and total RNAs were extracted from each biological replicate using a miRNeasy Mini Kit and a RNeasy Mini Kit, with resulting extracts treated with DNAse I (Qiagen) to remove contaminating genomic DNA. The first-strand cDNA was synthesized from mature miRNA and total RNA respectively with a miScript II RT kit (Qiagen) and a SuperScript™ III Reverse Transcriptase kit (Invitrogen, Carlsbad, CA) according to the manufacturer's instructions. RT-qPCR reactions were performed with a miScript SYBR 1 Green PCR Kit (Biorad, Hercules, CA) according to the manufacturer's instructions using miRNA-specific forward primers (S1 Table) with miScript II RT kit (Qiagen) products as a template. Analogously, RT-qPCR reactions for corresponding putative target mRNA transcripts were performed using SYBR 1 Green Master Mix (Biorad) according to the manufacturer's instructions using target transcript-specific forward and reverse primers (S1 Table), with products from the SuperScript™ III Reverse Transcriptase kit used as a template. All amplification reactions were performed on a StepOnePlus Real-Time PCR system (Applied Biosystems Inc., Foster City, CA), with three technical replicates across all biological replicates. Normalized miRNA and target transcript expression levels were calculated using the 2 −ΔΔC(t) method [67] with the internal references U6 snRNA and 5S rRNA for miRNA normalization and rp49 for mRNA normalization, respectively. A one-way ANOVA was used to examine the significance of expression differences between the two samples using XLSTAT software (Addinsoft, USA). The correlation coefficient was determined by Pearson's correlation analysis between the transcript levels of the selected eight miRNAs and their corresponding putative target genes.

miRNA library preparation, sequencing, and annotation
Six miRNA libraries (three biological replicates from each of 91-R and 91-C strains) were constructed, from which sequencing generated 216.2 million total raw reads (! 30.1 million per library; Table 1). After trimming reads (i.e., removal of reads without a 3' adaptor and < 15 nt) 5.1 to 5.8 and 4.4 to 7.7 million reads were retained respectively among triplicates from 91-R and 91-C (Table 1). These were then used in subsequent analyses. Considering trimmed reads across all six libraries, the majority of the sequences (55.5%) were distributed from 16-25 nt (S1 Fig), which is the standard size of described miRNAs [68]. A class of 26-31 nt long RNA sequences accounted for 25.4% of the total reads (34,598,261) and were classified as suspected piwi-interacting RNAs (piRNAs) [69]. Of the trimmed reads, 531,843 to 992,816 miRNA tags were identified across all replicate libraries of 91-R and 91-C (Table 1), and among these between 6,050 and 9,356 matched D. melanogaster miRNA entries in the miRBase R.21 database.
Annotation of known miRNAs in the 91-R and 91-C strains identified 163 precursor and 256 mature miRNAs following alignment of trimmed reads to the D. melanogaster precursor/ mature miRNAs in the miRBase R.21 (S2 Table). None of the precursor or mature miRNAs were expressed uniquely in either strain. However, we failed to identify any of the 93 known D. melanogaster precursors that are recorded in the miRBase R.21 (S3 Table). In both strains, miR-184-3p miRNA was the most abundant among the means of non-normalized reads across triplicate libraries, with miR-8-3p, miR-276a-3p, bantam-3p, and miR-33-5p as the next most abundant miRNAs in both strains ( Table 2). Algorithms in the miRDeep2 package identified 17 potential novel miRNAs in 91-R and 91-C, of which 15 were in common between libraries derived from both strains; putatively novel-miR-3L-18860981 and novel-miR-2R-20583765 were uniquely observed in 91-R and 91-C, respectively (S4 Table). The range of estimated MFE among potential novel miRNAs was between -31.7 and -25.5 kcal mol -1 and their mature lengths from 18-25 nt. These 17 novel miRNAs were named based on the chromosome and position on which the miRNA gene is located in the D. melanogaster genome.

Differential expression analysis and RT-qPCR validation
Comparison of normalized estimates of miRNA quantity (log 2 fold-changes; S5 Table) demonstrated an overall parity across replicated libraries derived from DDT resistance strain, 91-R, as compared to the susceptible control, 91-C (counts pooled across replicates within strain; S2 Fig). Exceptions were found among ten known miRNAs that showed significant levels of differential expression between 91-R and 91-C (log 2 fold-change ! |1|; FDR 0.05). Specifically, nine miRNAs were significantly down-regulated, and only one (miR-986-5p) was significantly up-regulated in 91-R (25.3-fold) when compared with the susceptible control 91-C (Table 3). Interestingly, four out of these nine down-regulated miRNAs [miR-310, miR-311, miR-312 and miR-313] are clustered miRNAs and belong to the miR-310 family. Overall, most of the differentially expressed miRNAs were down-regulated in the DDT resistant strain 91-R, suggesting that those down-regulated miRNAs may be potentially involved in DDT resistance.

Category
Analyses of total reads data

Target transcript predictions and correlated expression between miRNAs and target transcripts
Considering only the ten miRNAs predicted to be differentially expressed between 91-R and 91-C, a total of 46,368 miRNA-target pairs were predicted by the algorithms applied by ComiR (n = 5,577), DIANA (n = 2,548), and RNAhybrid (n = 38,243; S6 Table). Overlap in output occurred among 664 transcript targets predicted by all three algorithms. Functional GO annotations were received for 603 of these 664 transcript targets, with 64.3%, 20%, and 15.7% that respectively received terms in biological process (BP), molecular function (MF), and cellular component (CC) categories (GO level 2; S7 Table). Specifically, two putative functions (development and metabolism) were highly represented in the BP category, and the MF category showed the function of 'binding' and 'catalytic activity' was most prevalent among target genes. Moreover, cell, intracellular component, cytosol, cytoplasm, and cytoskeleton were largely overrepresented in CC category (Fig 2A). A total of 2,175 biological pathways were assigned to 258 of the 664 predicted transcript targets (38.9%; Fig 2B; S8 Table), and the target Ras85D (CG9375) is associated with regulation of tissue growth and development represented 93 biological pathways. The remaining 405 transcript targets (61.1%) received no pathway annotations. Among these 2,175 biological pathways, 58 (2.7%), 54 (2.5%), and 53 (2.4%) target genes of the known differentially expressed miRNAs, respectively, were assigned to metabolism, signal transduction, and metabolic pathway (Fig 2B). Additional annotation focused on a subset of the 664 transcripts putatively targeted by differentially expressed miRNAs; specifically, those transcripts likely to be involved in xenobiotic metabolism. These putative transcripts were predicted to target phase I, II, and III detoxification pathways such as cytochrome P450s, GSTs, esterases, and ABC transporters (43 target transcripts; Fig 3). Additionally, differentially expressed miRNAs were predicted to regulate other genes associated with cuticular proteins, acetylcholine receptors, nicotinic acetylcholine receptors, and glutamate-gated chloride channels (13 target transcripts; Fig 3). The phase I, cytochrome P450 targets Cyp6a8, Cyp6g1, Cyp6g2, and Cyp6w1, were previously associated with DDT resistance as described in the Introduction. Also of note, the down-regulated miR-310-313 cluster in 91-R strain was predicted to share several P450s, sodium channel proteins, and cuticular proteins encoding target genes (Fig 3). Additionally, sixteen P450 genes (MFE -25; Cyp18a1, Cyp305a1, Cyp309a2, Cyp312a1, Cyp313b1, Cyp49a1, Cyp4ae1, Cyp4g1, Cyp4g15, Cyp4p2, Cyp4s3, Cyp6a19, Cyp6g1, Cyp6g2, Cyp6v1, and Cyp9h1) were among the predicted targets of 9 of the 17 putative novel miRNAs (S9 Table).

Discussion
Resistant phenotypes have evolved within arthropod populations that lead to survival when exposed to various chemical or biological insecticidal toxins and seemingly arise to each successively introduced novel mode of action [70]. Although direct correlations have been made between single transposon-based mutations at single genetic loci (an Accord element associated with Cyp6g1) and corresponding insecticide resistance traits [71][72][73], in many instances the causal mutations and molecular mechanisms involved in resistance are yet to be fully understood. Moreover, definitive linkages between phenotypes and corresponding adaptive mutations remain difficult to define [74], especially in instances where phenotypes arise due to the contribution of multiple genes or gene interactions, or are complexed with variance due to Differentially expressed microRNAs are associated with DDT detoxification pathways the environment [75]. For instance, DDT resistance was linked to the up-regulation of Cyp6g1 caused by the upstream integration of an Accord transposon at the DDT-R locus in D. melanogaster populations [71] but shown to be independent of DDT-R in the highly DDT resistant strain 91-R [76]. The prediction of thirteen selective sweeps that became fixed in the genome of 91-R during the course of DDT selection implicated a complex polygenic mode of evolution [15], which likely involves a combination of cis-and trans-regulatory mutations that modulate the function of stress response and neurogenic pathways [17]. Regardless of evidence that strongly implicates miRNAs as potent modulators of gene expression at the post-transcriptional level [77], the role of miRNAs in the evolution of differential gene expression in insecticide resistant phenotypes among arthropods remains suggestive or associative in many cases [44,45].
The current study identified a total of 163 precursors with 256 mature known miRNAs and 17 novel miRNAs. This study failed to identify 93 known precursor miRNAs, previously recorded in the miRBase R.21. The extraction of sRNAs from 3-5 day-old virgin adult female samples may likely have biased the number of miRNA types. Previously, miRNAs were shown to be sex-biased in D. melanogaster, where expression was preferentially associated with the reproductive functions [78]. More recent work that compared miRNA expression between mature male and female reproductive organs in Bactrocera dorsalis demonstrated that sexbiased miRNAs are likely involved in sexual differentiation [79]. Furthermore, the expression of miRNAs varied across developmental stages of Xenopus laevis [80] and across different human tissues [81]. This suggests that our use of virgin female adults might have narrowed the pool of potential miRNAs that could be encountered within the resultant sRNA libraries but is justified since adults are exposed to DDT selection within the 91-R laboratory colony. Another hypothesis may reside in the potential saturation of sRNA samples with 2S ribosomal RNA (rRNA). Specifically, D. melanogaster rRNA is composed of four individual rRNAs, 28S, 18S, 5.8S, and 2S, with the 2S rRNA being 30-nt in length [82]. Analysis of our sRNA read data revealed that 32% of non-miRNA sequences matched the D. melanogaster 2S rRNA, suggesting the possibility that the 2S rRNA component of our libraries might have affected the subsequent read depths and the inability to identify 93 known miRNAs from 91-R and 91-C strains if they Putative xenobiotic metabolism-related target genes of differentially expressed miRNAs. The network consists of differentially expressed miRNAs and their corresponding target genes. Red hexagons represent differentially expressed miRNAs; blue circles represent target genes as phase I detoxification genes; green circles represent target genes encoding phase II detoxification genes; yellow circle represents target genes as phase III detoxification gene; and purple circles represent other genes affecting the insecticide insensitivity in insects. The network was generated and visualized in Cytoscape v3.6.0. Cyp6a8, Cyp6g1, Cyp6g2, and Cyp6w1 associated with DDT resistance are in bold. https://doi.org/10.1371/journal.pone.0196518.g003 Differentially expressed microRNAs are associated with DDT detoxification pathways would have been at low copy number. Neither of the above explanations can be ruled out but require additional investigation. The current study nonetheless predicted the significant quantitative difference in ten miRNA levels between 91-R and 91-C at the adult virgin female stage.
Ten differentially expressed miRNAs and their corresponding putative target transcripts were predicted and received GO and pathway database annotations in this study. Among these ten differentially expressed miRNAs, four members of the miR-310 family, miR-310, miR-311, miR-312 and miR-313, were significantly down-regulated in 91-R when compared with the 91-C strain. The miR-310 family form a cluster on chromosomal arm 2R between the CDS of qsm and Nnf1a (positions 20,583,752 to 20,584,260) [83] and arose via duplication from the ancestral miR-91 [84]. Previous studies showed that the miR-310 family regulates genes in the D. melanogaster Toll-mediated innate immune response pathway [85], via hedgehog signaling [86], and beta-catenin that in turn regulates cell adhesion and outgrowth [87]. Furthermore, the mir-310-313 cluster has been reported to be associated with hypersensitivity to nicotine in D. melanogaster [88]. These authors observed that the expression of the miR-310-313 downregulates escargot (esg) gene expression involved in the development of sensory organs and neurons in the thoracoabdominal ganglion. The overexpression of the miR-310-313 induces an abnormal sensitivity to nicotine by abating esg transcription, then disrupting sensory organs involved in chemical perception and cuticle development. This mechanism may be similar to the DDT resistance of miR-310 family members in 91-R strain. Additionally, miR-310 family members impact synaptic functions through the regulation of neurotransmitter release at neuromuscular junctions of D. melanogaster larvae [89]. These authors observed that miR-310-313 cluster regulates neurotransmitter release at presynaptic terminals by decreasing the expression of Kinesin heavy chain-73, Khc-73, which is involved in early neural development. Interestingly, our prior research validated the significant differential regulation of genes directing neuronal development in 91-R, including Unc-115b, and CG31832 that function in neural growth and development, but detected no significant difference in expression of Khc-73 [17]. Our prior results suggest the possibility that adaptive responses to multigenerational neurotoxic DDT selection in 91-R may affect the function of neuromuscular junctions, although the role of the miR-310 family in the regulation of genes associated with neuronal functions in 91-R remains a hypothesis that needs to be tested experimentally.
Findings from this study show that the down-regulation of miR-311-3p, miR-312-3p, and miR-313-3p is correlated with the up-regulation of their respective in silico predicted cytochrome P450 target transcripts, Cyp6a8, Cyp4s3, Cyp4ae1, Cyp6g1, Cyp6g2, Cyp6t3, Cyp6v1, Cyp18a1, Cyp49a1, Cyp303a1, Cyp309a2, Cyp313a4, and Cyp313b1 in 91-R. This up-regulation of Cyp6g1 target is predicted to occur via a decreased miR-310-313 expression in 91-R and also corresponds to the cytochrome P450 initially implicated in conferring DDT resistance at the DDT-R locus via integration of the Accord transposon among D. melanogaster populations [71]. Our computational and empirical data suggests that the miR-310 family members may be involved in the posttranscriptional regulation of these cytochrome P450s, which in turn might be directly involved in DDT detoxification or stress response [17]. Thus, the involvement of miR-310-313 could explain prior evidence that even though DDT resistance in 91-R is independent of the Accord insertion [76], significant up-regulation of Cyp6g1, as well as Cyp6a8, Cyp4s3, Cyp4ae1, Cyp6g2, Cyp6t3, Cyp6v1, Cyp18a1, Cyp49a1, Cyp303a1, Cyp309a2, Cyp313a4, and Cyp313b1 occurs constitutively in the strain. Regardless of these strong correlations, additional functional analyses are required to confirm these predicted impacts of miRNA-based posttranscriptional regulation.
The function of the most highly up-regulated miRNA in 91-R strain, miR-986-5p, is currently unknown, but levels in the hemolymph of adult virgin males are known to significantly decrease over time [90]. Interestingly, the miR-986 precursor is located on the second chromosome within the third intron of the Cyp4e2 gene, which is involved in the metabolism of endogenous and exogenous compounds [91]. Our transcript target site predictions suggest that miR-986-5p could interact with transcripts of cytochrome P450s, GSTs, esterases, and superoxide dismutases (SODs). Specifically, the GSTs and SODs are a group of a multifunctional antioxidant enzymes that play an important role in mediating oxidative stress caused by reactive oxygen species (ROS) in insects [92][93][94]. Therefore, the miR-986 may be involved in the posttranslational regulation of genes that alleviate oxidative stress induced by DDT insecticide exposures. The putative targeted transcripts of miR-986-5p include Cyp6g1 and Cyp6g2. However, the RT-qPCR performed here, as well as prior RNA-seq results [17], indicate that the expression levels of putative targets Cyp6g1 and Cyp6g2 were not decreased. In contrast, both Cyp6g1 and Cyp6g2 were significantly up-regulated in 91-R, suggesting miR-986-5p may have a transcript stabilizing or enhancing effect, as was demonstrated previously [95,96]. Alternatively, algorithms used for in silico prediction of miRNA-transcript target interactions can produce variable results depending on the input database and models that are applied [97,98]. Regardless, in vivo or in vitro validation of these assumptions will be required.

Conclusion
Analyses conducted in this study focused on differentially expressed miRNAs that were predicted to regulate transcript levels of both phase I, II and III detoxification genes previously shown to be associated with the DDT resistance phenotype. Cytochrome P450s that are involved in many cellular processes, including xenobiotic detoxification, have been studied [99] for associations between miRNA levels and corresponding putatively targeted P450 transcripts. For example, a negative relationship was shown between up-regulated miRNAs miR-8534-5p and miR-375-5p and their respective predicted targeted cytochrome P450s, Cyp6b6 and Cyp4g15, in chlorantraniliprole-resistant strains of Plutella xylostella [44]. Several miRNAs down-regulated in deltamethrin-resistant mosquitoes played a role in pyrethroid resistance through the reduced targeting of Cyp325bg3, Cyp6n23, and Cyp9j35 [47,48,50], while the miR-NAs (miR-155, miR-216b, miR-499) modulate the abundance of Cyp561d2 transcripts in response to fipronil exposure [100]. Additionally, miR-285 and miR-278 differentially regulate Cyp6n23 and Cyp6ag11 in pyrethroid resistant compared to susceptible Culex pipiens [48,49].
Thus, our implication of down-regulated miR-311-3p, miR-312-3p, and miR-313-3p with the corresponding constitutive up-regulation of in silico predicted targets Cyp6a8, Cyp4s3, Cyp4ae1, Cyp6g2, Cyp6t3, Cyp6v1, Cyp18a1, Cyp49a1, Cyp303a1, Cyp309a2, Cyp313a4 in 91-R may provide yet another example of a miRNA-mediated posttranscriptional modification that contributes to an insecticide resistance trait. Additionally, 4 out of 10 miRNAs predicted to be differentially expressed between 91-R and 91-C (miR-986-5p, miR-995-3p, miR-312-3p, miR-2a-3p) were also predicted to interact with and impact the transcript level of multidrug resistance-associated protein B7 (ABC-B7). The ABC-B subfamily member, MDR49, had been previously been reported to show significant levels of differential expression between 91-R and the DDT-susceptible strain Canton-S [12] but not between 91-R and 91-C [12,17]. It was subsequently shown that the 91-R-derived 91-R-MDR49B allele provided increased DDT tolerance via transgenic expression in susceptible D. melanogaster, implicating structural, as opposed to dosage, effects on the DDT resistance trait [14]. Thus, one could hypothesize that ABC-B7 might be regulated by one or multiple miRNAs that are differentially expressed between 91-R and 91-C. In addition to further empirical study, however, it should also be kept in mind that multiple different mutational mechanisms can contribute to DDT resistance in 91-R.
In this study, we identified a set of differentially expressed known miRNAs and several novel miRNAs from DDT-resistant 91-R and -susceptible 91-C strains. Since 91-R and 91-C strains have a common initial origin (genetic background), and 91-R has been selected for survival when exposed to chronic high levels of DDT for over six decades, changes between the strains are speculated to result from either the effects of random genetic drift or from selection. The experiments described here do not allow the differentiation of any impacts from drift versus selection, but the strong correlation between differential expression of miRNAs and their corresponding in silico predicted target transcripts suggests the potential for the involvement of posttranscriptional regulation of several detoxification genes.
Moreover, the experimental procedures do not allow the formulation that down-regulated miRNAs identified in the present study could lead directly to the overexpression of detoxification genes in the 91-R strain, thus requiring further functional experiments in order to elucidate the mechanisms of miRNAs involvement in DDT resistance in 91-R. As such, the overexpressed P450s putatively targeted by differentially-regulated miRNAs in 91-R strain need to be further evaluated to identify the ones directly involved in the DDT resistance mechanism. Due to the polygenic nature of DDT resistance in 91-R (see Introduction), however, this trait may be the result of interactions or additive/non-additive effects from several distinct genetic factors. This study, for the first time, suggests that these might include impacts mediated by miRNA posttranscriptional regulation. Our results provide valuable information for exploring the mechanisms of miRNAs involved in insecticide resistance and for understanding the evolution of post-transcriptional regulation in response to DDT pressure in D. melanogaster.