The small RNA sequencing data sets generated during this study are available in the NCBI Sequence Read Archive (SRA) under accession numbers SRP136631 (<https://www.ncbi.nlm.nih.gov/sra/> \<<https://www.ncbi.nlm.nih.gov/sra/>\>).

Introduction {#sec001}
============

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 and protect human health \[[@pone.0196518.ref001]\], this intensive use of chemical insecticides has also led to the development of resistance to one or more classes of insecticides \[[@pone.0196518.ref002]--[@pone.0196518.ref004]\]. The neurotoxic organochlorine insecticide, dichlorodiphenyltrichloroethane (DDT), has been extensively used to control many agricultural insect pests and insects that vector human diseases, but was banned in most countries in the 1980s due to environmental concerns \[[@pone.0196518.ref005]\]. The genetic basis of DDT resistance in the *Drosophila melanogaster* (*D*. *melanogaster*) provides a model system for studying the evolution of insecticide resistance. Indeed, the low-level DDT resistance phenotype in *D*. *melanogaster* is thought to be associated with a single cytochrome P450, *Cyp6g1* \[[@pone.0196518.ref006],[@pone.0196518.ref007]\]. However, moderate to high-level DDT resistance is polygenic \[[@pone.0196518.ref008]\] involving modulation in DDT penetrance and excretion and multiple differentially-regulated Phase I, II, and III detoxification genes/enzymes, respectively, including P450s, glutathione S transferases (GSTs), and ATP binding cassette (ABC) transporters. Regarding the P450s, members of the *Cyp6* (*Cyp6g1*, *Cyp6g2*, *Cyp6a2*, and *Cyp6w1*) and *Cyp12* (*Cyp12d1*) subfamilies have been implicated in polygenic DDT resistance \[[@pone.0196518.ref009]--[@pone.0196518.ref013]\]. Furthermore, the overexpression and structural changes in membrane-spanning ABC transporters were shown to be involved in DDT efflux and contribute to the overall DDT resistance phenotype in *D*. *melanogaster* \[[@pone.0196518.ref014]\]. Full genome re-sequencing identified 13 different regions showing evidence of high nucleotide diversity, directional selection (selective sweeps), and thus putatively associated with the evolution of DDT resistance in the DDT-resistant strain, *91-R*, when compared with a susceptible isoline, *91-C* \[[@pone.0196518.ref015]\]. Analyses of this same whole genome re-sequencing data also identified a large panel of fixed amino acid changing mutations between *91-R* and *91-C*, many of which showed evidence of directional selection \[[@pone.0196518.ref016]\]. Most recently, a multigenic response to DDT selection was demonstrated within *91-R* where transcripts with constitutive differential-expression compared to *91-C* were enriched in cell survival, stress response, and neurological functions \[[@pone.0196518.ref017]\]. Additionally, reduced expression of the putative calcium/lipid binding domain-containing protein from gene model *CG10737* is associated with short-term (3--5 hrs) DDT knockdown resistance within the *Drosophila* Genetic Reference Panel, suggesting a role for ameliorated effects of muscle overstimulation in DDT resistance phenotypes \[[@pone.0196518.ref013]\]. Despite these lines of evidence, the systemic basis and underlying genetic control of polygenic DDT resistance remains enigmatic.

Accumulating evidence demonstrates that the post-transcriptional regulation (RNA editing and alternative splicing) may contribute to the evolution of insecticide resistant phenotypes \[[@pone.0196518.ref018],[@pone.0196518.ref019]\]. The miRNAs are a class of endogenous 18--25 nt non-coding small RNAs that, since their discovery in *Caenorhabditis elegans* over two decades ago \[[@pone.0196518.ref020]\], have been identified across several arthropod species including twelve *Drosophila* species \[[@pone.0196518.ref021]\]. Biogenesis of miRNAs occurs from transcript-derived stem-loop structures with canonical CNNC downstream elements \[[@pone.0196518.ref022]\] that are processed into a \~70 bp double-stranded precursor RNAs (pre-RNAs) within the nucleus by the ribonuclease type III enzyme Drosha \[[@pone.0196518.ref023]\]. After translocation to the cytoplasm, miRNA precursors exhibit canonical stem-loop secondary structures, which usually have two arms called the miRNA-5p arm or -3p arm \[[@pone.0196518.ref024]\]. The pre-RNAs are further processed into mature 21--22 nt miRNAs by Dicer \[[@pone.0196518.ref025],[@pone.0196518.ref026]\]. Following the degradation of the passenger strand \[[@pone.0196518.ref027]\], the resulting single-stranded guide sequence is loaded into the RNA-induced silencing complex (RISC) \[[@pone.0196518.ref028]\]. Subsequent base pairing of the miRNA guide "seed sequence" (nucleotides 2--7 at the 5\' end of guide RNAs) \[[@pone.0196518.ref029]\] with the 3\'-untranslated region (3\'-UTR) of target cellular mRNA transcripts most often targets the mRNA for degradation by RISC \[[@pone.0196518.ref030]\], either blocking translation initiation factor binding or inhibiting elongation factor progression that leads to premature termination \[[@pone.0196518.ref031]--[@pone.0196518.ref033]\]. Alternatively, a miRNA binding of a target transcript can result in translational up-regulation \[[@pone.0196518.ref034]\] within their mRNA resulting in inhibition of translation or mRNA degradation \[[@pone.0196518.ref035]\]. Overall, miRNAs are recognized as potent regulators of eukaryotic gene expression at the post-transcriptional level \[[@pone.0196518.ref036],[@pone.0196518.ref037]\] that regulate cell differentiation \[[@pone.0196518.ref038]\], migration \[[@pone.0196518.ref039]\], and neuronal development \[[@pone.0196518.ref040]\]. Approximately 256 *D*. *melanogaster* miRNA precursors have been identified and deposited in the mirBase (Release 21) \[[@pone.0196518.ref024]\] and have been shown to be involved in regulating biological processes such as development \[[@pone.0196518.ref041]\], immune response \[[@pone.0196518.ref042]\], and metabolism \[[@pone.0196518.ref043]\].

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 \[[@pone.0196518.ref044],[@pone.0196518.ref045]\]. Additionally, miRNAs differentially expressed between insecticidal *Bacillus thuringiensis* (*Bt*) toxin-resistant and -susceptible strains of *Ostrinia furnacalis* were predicted to target potential receptor genes \[[@pone.0196518.ref046]\]. In *Culex pipiens*, miRNAs differentially expressed between deltamethrin-sensitive and -resistant strains were proposed to mediate the expression of putative cytochrome P450 target genes \[[@pone.0196518.ref047]\]. Specifically, *Culex pipiens miR-285* and *miR-278* were implicated in pyrethroid resistance through the transcriptional regulation of *Cyp6n23* and *Cyp6ag11* \[[@pone.0196518.ref048],[@pone.0196518.ref049]\], as well as a miRNA cluster involved in regulation of *Cyp9j35* and *Cyp325bg3* \[[@pone.0196518.ref050]\] and upregulation of *miR-932* that regulates transcript levels of the cuticular gene *CpCPR5* \[[@pone.0196518.ref051],[@pone.0196518.ref052]\].

To date, the authors have no knowledge of any investigation of the contributions of post-transcriptional 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 DDT-susceptible 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.

Materials and methods {#sec002}
=====================

miRNA library preparation, sequencing, and annotation {#sec003}
-----------------------------------------------------

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 \[[@pone.0196518.ref053]\]. 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* \[[@pone.0196518.ref017]\].

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 miRBase 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 \[[@pone.0196518.ref054]\]. Illumina adapters of sRNA raw sequence reads were trimmed using cutadapt v.1.15 \[[@pone.0196518.ref055]\]. 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 \[[@pone.0196518.ref056]\] 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](http://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 miRNAs (mature.fa) using default parameters suggested by the developers \[[@pone.0196518.ref057]\]. The structure and minimal free energy (MFE) of all potential novel miRNAs were predicted using a RNAfold \[[@pone.0196518.ref058]\] with algorithms described previously \[[@pone.0196518.ref059]\]. 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 \[[@pone.0196518.ref060]\].

Differential expression of miRNAs {#sec004}
---------------------------------

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 \[[@pone.0196518.ref061]\] 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 \[[@pone.0196518.ref062]\]. 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 \[[@pone.0196518.ref063]\].

Target prediction and functional annotation of differentially expressed miRNAs {#sec005}
------------------------------------------------------------------------------

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 \[[@pone.0196518.ref064]\], DIANA \[[@pone.0196518.ref065]\], and ComiR \[[@pone.0196518.ref066]\], 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*. *melanogaster* (dmel-all-three_prime_UTR-r6.19.fasta at [http://www.flybase.org](http://www.flybase.org/)) with MFE ≤ -30 kcal mol^-1^; (2) DIANA: miTG score ≥ 0.8; and (3) ComiR: high threshold ≥ 0.8. Additionally, target gene ontology (GO) and corresponding pathways for all putative target transcripts were retrieved from the FlyMine database ([http://www.flymine.org](http://www.flymine.org/)), and the GO terms were classified with CateGOrizer ([https://www.animalgenome.org](https://www.animalgenome.org/)).

Validation and correlated expression between miRNAs and target transcripts {#sec006}
--------------------------------------------------------------------------

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*, *Esterase7*, *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*) and---identically 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^®^ Green PCR Kit (Biorad, Hercules, CA) according to the manufacturer's instructions using miRNA-specific forward primers ([S1 Table](#pone.0196518.s003){ref-type="supplementary-material"}) with miScript II RT kit (Qiagen) products as a template. Analogously, RT-qPCR reactions for corresponding putative target mRNA transcripts were performed using SYBR^®^ Green Master Mix (Biorad) according to the manufacturer's instructions using target transcript-specific forward and reverse primers ([S1 Table](#pone.0196518.s003){ref-type="supplementary-material"}), 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 \[[@pone.0196518.ref067]\] 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.

Results {#sec007}
=======

miRNA library preparation, sequencing, and annotation {#sec008}
-----------------------------------------------------

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](#pone.0196518.t001){ref-type="table"}). 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](#pone.0196518.t001){ref-type="table"}). 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](#pone.0196518.s001){ref-type="supplementary-material"}), which is the standard size of described miRNAs \[[@pone.0196518.ref068]\]. 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) \[[@pone.0196518.ref069]\]. 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](#pone.0196518.t001){ref-type="table"}), and among these between 6,050 and 9,356 matched *D*. *melanogaster* miRNA entries in the miRBase R.21 database.

10.1371/journal.pone.0196518.t001

###### Small RNA sequences from *91-C* and *91-R* triplicate libraries read data.

![](pone.0196518.t001){#pone.0196518.t001g}

  ------------------------------------------------------------- ---------------------------------- -------------- -------------- -------------- -------------- ------------
  **Category**                                                  **Analyses of total reads data**                                                               

  ***91-R*-1**                                                  ***91-R*-2**                       ***91-R*-3**   ***91-C*-1**   ***91-C*-2**   ***91-C*-3**   

  Raw reads                                                     35,712,203                         30,108,348     33,996,366     43,821,999     34,205,253     38,365,489

  Trimmed reads                                                 5,801,176                          5,105,542      5,439,260      7,738,822      6,065,134      4,448,327

  Annotated with ncRNA (rRNA, tRNA, snRNA, snoRNA and others)   4,399,524                          4,072,155      4,196,227      5,961,313      4,553,746      3,182,916

  *D*. *melanogaster* miRNAs                                    384,620                            252,554        339,001        455,995        336,210        356,146

  Unannotated                                                   1,017,032                          780,833        904,032        1,321,514      1,175,178      909,265

  **Data Processing/**\                                         **Analyses of unique reads**                                                                   
  **Strains**                                                                                                                                                  

  ***91-R*-1**                                                  ***91-R*-2**                       ***91-R*-3**   ***91-C*-1**   ***91-C*-2**   ***91-C*-3**   

  Unique miRNA                                                  992,816                            920,013        944,922        900,308        894,908        531,843

  Annotated with ncRNA (rRNA, tRNA, snRNA, snoRNA and others)   490,904                            510,823        503,922        344,222        375,672        174,207

  *D*. *melanogaster* miRNAs                                    7,643                              6,050          7,079          9,356          7,947          7,942

  Unannotated                                                   494,269                            403,140        433,921        546,730        511,289        349,694
  ------------------------------------------------------------- ---------------------------------- -------------- -------------- -------------- -------------- ------------

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](#pone.0196518.s004){ref-type="supplementary-material"}). 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](#pone.0196518.s005){ref-type="supplementary-material"}). 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](#pone.0196518.t002){ref-type="table"}). 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](#pone.0196518.s006){ref-type="supplementary-material"}). 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.

10.1371/journal.pone.0196518.t002

###### The most abundant reads from *91-R* and *91-C* small RNA libraries corresponding to known *Drosophila melanogaster* miRNAs in miRBase R.21.

![](pone.0196518.t002){#pone.0196518.t002g}

  ----------------------------------------------------------------------------------------------------------------------------------------------------
  miRNA           *91-R*\                                             *91-C*\                                             Mature sequence
                  Read counts[^a^](#t002fn001){ref-type="table-fn"}   Read counts[^a^](#t002fn001){ref-type="table-fn"}   
  --------------- --------------------------------------------------- --------------------------------------------------- ----------------------------
  *miR-184-3p*    27,954                                              27,402                                              `UGGACGGAGAACUGAUAAGGGC`

  *miR-8-3p*      24,289                                              25,956                                              `UAAUACUGUCAGGUAAAGAUGUC`

  *miR-276a-3p*   20,174                                              17,849                                              `UAGGAACUUCAUACCGUGCUCU`

  *bantam-3p*     7,498                                               9,036                                               `UGAGAUCAUUUUGAAAGCUGAUU`

  *miR-33-5p*     4,208                                               5,040                                               `GUGCAUUGUAGUCGCAUUGUC`

  *miR-10-5p*     3,464                                               3,954                                               `ACCCUGUAGAUCCGAAUUUGUU`

  *miR-317-3p*    2,908                                               3,292                                               `UGAACACAGCUGGUGGUAUCCAGU`

  *miR-14-3p*     2,603                                               2,611                                               `UCAGUCUUUUUCUCUCUCCUAU`

  *miR-31a-5p*    3,580                                               2,354                                               `UGGCAAGAUGUCGGCAUAGCUGA`

  *miR-312-3p*    1,056                                               2,318                                               `UAUUGCACUUGAGACGGCCUGA`

  *miR-11-3p*     1,659                                               2,275                                               `CAUCACAGUCUGAGUUCUUGC`

  *miR-311-3p*    596                                                 2,170                                               `UAUUGCACAUUCACCGGCCUGA`

  *miR-318-3p*    994                                                 1,977                                               `UCACUGGGCUUUGUUUAUCUCA`

  *miR-999-3p*    1,604                                               1,660                                               `UGUUAACUGUAAGACUGUGUCU`

  *miR-957-3p*    1,744                                               1,576                                               `UGAAACCGUCCAAAACUGAGGC`

  *miR-276b-3p*   1,846                                               1,553                                               `UAGGAACUUAAUACCGUGCUCU`

  *miR-277-3p*    2,046                                               1,489                                               `UAAAUGCACUAUCUGGUACGACA`

  *miR-995-3p*    542                                                 1,275                                               `UAGCACCACAUGAUUCGGCUU`

  *miR-305-3p*    1,000                                               1,235                                               `CGGCACAUGUUGAAGUACACUCA`

  *let-7-5p*      1,094                                               1,208                                               `UGAGGUAGUAGGUUGUAUAGU`

  *miR-986-5p*    3,069                                               139                                                 `UCUCGAAUAGCGUUGUGACUGA`

  *miR-958-3p*    1,020                                               899                                                 `UGAGAUUCUUCUAUUCUACUUU`
  ----------------------------------------------------------------------------------------------------------------------------------------------------

^a^ Non-normalized reads summed across triplicate libraries

Differential expression analysis and RT-qPCR validation {#sec009}
-------------------------------------------------------

Comparison of normalized estimates of miRNA quantity (log~2~ fold-changes; [S5 Table](#pone.0196518.s007){ref-type="supplementary-material"}) 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](#pone.0196518.s002){ref-type="supplementary-material"}). 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](#pone.0196518.t003){ref-type="table"}). 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.

10.1371/journal.pone.0196518.t003

###### Differentially expressed miRNAs between *91-C* and *91-R*.

![](pone.0196518.t003){#pone.0196518.t003g}

  ------------------------------------------------------------------------------------------------------------------------------------------------------------------------
  miRNA              *p*-value   FDR[^a^](#t003fn001){ref-type="table-fn"}   Fold-change   Log~2~\                                             miRNA sequence
                                                                                           fold-change[^b^](#t003fn002){ref-type="table-fn"}   
  ------------------ ----------- ------------------------------------------- ------------- --------------------------------------------------- ---------------------------
  *miR-986-5p*       3.0E-53     7.8E-51                                     25.3          4.66                                                `TCTCGAATAGCGTTGTGACTGA`

  *miR-2a-1-3p//*\   7.0E-04     1.3E-02                                     -2.01         -1.01                                               `TATCACAGCCAGCTTTGATGAGC`
  *miR-2a-2-3p*                                                                                                                                

  *miR-312-3p*       1.4E-07     5.3E-06                                     -2.02         -1.01                                               `TATTGCACTTGAGACGGCCTGA`

  *miR-995-3p*       2.8E-09     1.4E-07                                     -2.17         -1.12                                               `TAGCACCACATGATTCGGCTT`

  *miR-286-3p*       2.1E-03     3.2E-02                                     -2.29         -1.19                                               `TGACTAGACCGAACACTCGTGCT`

  *miR-92a-3p*       2.9E-10     1.9E-08                                     -2.34         -1.23                                               `CATTGCACTTGTCCCGGCCTAT`

  *miR-4919-3p*      1.4E-04     3.2E-03                                     -2.57         -1.36                                               `TAATCCCTGAACGACTTGCAG`

  *miR-311-3p*       1.8E-15     1.5E-13                                     -3.38         -1.76                                               `TATTGCACATTCACCGGCCTGA`

  *miR-310-3p*       9.1E-16     1.2E-13                                     -3.94         -1.98                                               `TATTGCACACTTCCCGGCCTTT`

  *miR-313-3p*       4.4E-07     1.4E-05                                     -4.89         -2.29                                               `TATTGCACTTTTCACAGCCCGA`
  ------------------------------------------------------------------------------------------------------------------------------------------------------------------------

^a^ FDR: False discovery rate. Differentially expressed miRNAs were identified at the threshold \[FDR \< 0.05 and log~2~(fold change) ≥\|1.0\|\] of *91-C*/*91-R*.

^b^ Fold change was calculated as *91-C*/*91-R*.

The RT-qPCR validation of the predicted significant quantitative differences in miRNA levels among eight known and four novel miRNAs were amplified showing that the expression levels of *miR-986-5p* were highly up-regulated in *91-R*, whereas *miR-286-3p*, *miR-4919-3p*, *miR-311-3p*, *miR-312-3p*, and *miR-313-3p* were significantly down-regulated in *91-R* ([Fig 1](#pone.0196518.g001){ref-type="fig"}). However, levels of *miR-92a-3p* and *miR-310-3p* showed no significant difference between strains (*P*-value \> 0.05). Based on Pearson's correlation coefficient test, the expression patterns of selected miRNAs showed a similar trend between the results of sRNA sequencing and RT-qPCR (R^2^ = 0.971, *P*-value \< 0.01), confirming the predicted differential expression between strains. Additionally, the expression level of one putative novel miRNA (*novel-miR-3L-10365243*) was significantly down-regulated in *91-R*, whereas three other putative novel miRNAs showed no significant differences between two strains ([Fig 1](#pone.0196518.g001){ref-type="fig"}).

![RT-qPCR validation of differentially expressed known and novel miRNAs identified between *91-R* and *91-C* strains.\
The expression levels were normalized by both *U6* and *5S rRNA*. Statistical significance was analyzed using one-way ANOVA. The asterisks represent significance, where one asterisk indicates *P* \< 0.05.](pone.0196518.g001){#pone.0196518.g001}

Target transcript predictions and correlated expression between miRNAs and target transcripts {#sec010}
---------------------------------------------------------------------------------------------

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](#pone.0196518.s008){ref-type="supplementary-material"}). 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](#pone.0196518.s009){ref-type="supplementary-material"}). 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](#pone.0196518.g002){ref-type="fig"}). A total of 2,175 biological pathways were assigned to 258 of the 664 predicted transcript targets (38.9%; [Fig 2B](#pone.0196518.g002){ref-type="fig"}; [S8 Table](#pone.0196518.s010){ref-type="supplementary-material"}), 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](#pone.0196518.g002){ref-type="fig"}).

![Gene ontology (GO) analysis and pathway annotation of the most targeted genes of the known differentially expressed miRNAs between *91-R* and *91-C* strains.\
(A) GO analysis. 664 target genes from miRNAs predicted to be differentially expressed were analyzed with FlyMine to obtain the GO terms, and the GO terms were classified with CateGOrizer and separated into three major categories. (B) Pathway annotations. 664 target genes from miRNAs predicted to be differentially expressed were submitted to FlyMine to get the pathway classification. The 30 most abundant pathways were represented.](pone.0196518.g002){#pone.0196518.g002}

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](#pone.0196518.g003){ref-type="fig"}). 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](#pone.0196518.g003){ref-type="fig"}). 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](#pone.0196518.g003){ref-type="fig"}). 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](#pone.0196518.s011){ref-type="supplementary-material"}).

![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.](pone.0196518.g003){#pone.0196518.g003}

Moreover, the RT-qPCR-estimated quantities of 26 putatively targeted P450 transcripts showed an inverse relationship with the expression level of the corresponding miRNA(s) predicted to target them in *91-R*. For example, the relative expression of three miRNAs (*miR-311-3p*, *miR-312-3p*, and *miR-313-3p*) were significantly down-regulated in *91-R* strain ([Fig 1](#pone.0196518.g001){ref-type="fig"}), while the corresponding predicted targets (*Cyp6a8*, *Cyp4s3*, *Cyp4ae1*, *Cyp6g1*, *Cyp6g2*, *Cyp6t3*, *Cyp6v1*, *Cyp18a1*, *Cyp49a1*, *Cyp303a1*, *Cyp309a2*, *Cyp313a4*, and *Cyp313b1*) were significantly up-regulated ([Fig 4](#pone.0196518.g004){ref-type="fig"}). Furthermore, other corresponding predicted detoxification targets, *Esterase8*, *Esterase10*, *GstD1*, *GstE10*, *GstS1*, *GABA-R*, *mAChR-A*, *nAcRalpha-A* and *GluR-IB* were significantly up-regulated with the down-regulation of *miR-286-3p*, *miR-2a-3p*, *miR-311-3p*, *miR-312-3p*, and *miR-313-3p* in *91-R* strain ([Fig 4](#pone.0196518.g004){ref-type="fig"}). The resulting coefficient of correlation showed that the strong negative correlation between miRNAs and putative target transcripts (R = -0.9558; *P*-value \< 0.01) suggests the possibility that higher levels of detoxification transcripts may be influenced by the reduced levels of corresponding targeting miRNAs within the DDT resistant *91-R* strain. The precise role of these detoxification genes in mediating DDT resistance remains unknown, but may include increasing cellular resistance to oxidative stress, cuticular penetrance, or other systemic responses to exposure.

![Expression level analysis of potential target detoxification genes of differentially expressed miRNAs between *91-R* and *91-C* strains.\
The expression levels were normalized by *rp49*. Statistical significance was analyzed using one-way ANOVA. The asterisks represent significance, with a single asterisk indicating *P* \< 0.05.](pone.0196518.g004){#pone.0196518.g004}

Discussion {#sec011}
==========

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 \[[@pone.0196518.ref070]\]. 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 \[[@pone.0196518.ref071]--[@pone.0196518.ref073]\], 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 \[[@pone.0196518.ref074]\], especially in instances where phenotypes arise due to the contribution of multiple genes or gene interactions, or are complexed with variance due to the environment \[[@pone.0196518.ref075]\]. 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 \[[@pone.0196518.ref071]\] but shown to be independent of DDT-R in the highly DDT resistant strain *91-R* \[[@pone.0196518.ref076]\]. 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 \[[@pone.0196518.ref015]\], which likely involves a combination of *cis*- and *trans*-regulatory mutations that modulate the function of stress response and neurogenic pathways \[[@pone.0196518.ref017]\]. Regardless of evidence that strongly implicates miRNAs as potent modulators of gene expression at the post-transcriptional level \[[@pone.0196518.ref077]\], the role of miRNAs in the evolution of differential gene expression in insecticide resistant phenotypes among arthropods remains suggestive or associative in many cases \[[@pone.0196518.ref044],[@pone.0196518.ref045]\].

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 \[[@pone.0196518.ref078]\]. More recent work that compared miRNA expression between mature male and female reproductive organs in *Bactrocera dorsalis* demonstrated that sex-biased miRNAs are likely involved in sexual differentiation \[[@pone.0196518.ref079]\]. Furthermore, the expression of miRNAs varied across developmental stages of *Xenopus laevis* \[[@pone.0196518.ref080]\] and across different human tissues \[[@pone.0196518.ref081]\]. 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 \[[@pone.0196518.ref082]\]. 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 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) \[[@pone.0196518.ref083]\] and arose via duplication from the ancestral *miR-91* \[[@pone.0196518.ref084]\]. Previous studies showed that the miR-310 family regulates genes in the *D*. *melanogaster* Toll-mediated innate immune response pathway \[[@pone.0196518.ref085]\], via hedgehog signaling \[[@pone.0196518.ref086]\], and beta-catenin that in turn regulates cell adhesion and outgrowth \[[@pone.0196518.ref087]\]. Furthermore, the *mir-310-313* cluster has been reported to be associated with hypersensitivity to nicotine in *D*. *melanogaster* \[[@pone.0196518.ref088]\]. 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 \[[@pone.0196518.ref089]\]. 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* \[[@pone.0196518.ref017]\]. 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 \[[@pone.0196518.ref071]\]. 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 \[[@pone.0196518.ref017]\]. 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 \[[@pone.0196518.ref076]\], 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 \[[@pone.0196518.ref090]\]. 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 \[[@pone.0196518.ref091]\]. 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\[[@pone.0196518.ref092]--[@pone.0196518.ref094]\]. 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 \[[@pone.0196518.ref017]\], 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 \[[@pone.0196518.ref095],[@pone.0196518.ref096]\]. 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 \[[@pone.0196518.ref097],[@pone.0196518.ref098]\]. Regardless, *in vivo* or *in vitro* validation of these assumptions will be required.

Conclusion {#sec012}
==========

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 \[[@pone.0196518.ref099]\] 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* \[[@pone.0196518.ref044]\]. Several miRNAs down-regulated in deltamethrin-resistant mosquitoes played a role in pyrethroid resistance through the reduced targeting of *Cyp325bg3*, *Cyp6n23*, and *Cyp9j35* \[[@pone.0196518.ref047],[@pone.0196518.ref048],[@pone.0196518.ref050]\], while the miRNAs (*miR-155*, *miR-216b*, *miR-499*) modulate the abundance of *Cyp561d2* transcripts in response to fipronil exposure \[[@pone.0196518.ref100]\]. Additionally, *miR-285* and *miR-278* differentially regulate *Cyp6n23* and *Cyp6ag11* in pyrethroid resistant compared to susceptible *Culex pipiens* \[[@pone.0196518.ref048],[@pone.0196518.ref049]\].

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* \[[@pone.0196518.ref012]\] but not between *91-R* and *91-C* \[[@pone.0196518.ref012],[@pone.0196518.ref017]\]. 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 \[[@pone.0196518.ref014]\]. 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](#sec001){ref-type="sec"}), 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*.

Supporting information {#sec013}
======================

###### Size distribution of small RNAs from the *91-R* and *91-C* libraries.

(TIF)

###### 

Click here for additional data file.

###### Scatter plot of differentially expressed miRNA in the *91-R* and *91-C* strains.

Each plot represents a miRNA. The X- and Y-axis show the normalized read counts of miRNAs in the two strains respectively.

(TIF)

###### 

Click here for additional data file.

###### The list of RT-qPCR primers.

(XLSX)

###### 

Click here for additional data file.

###### Known miRNAs identified in this research.

(XLSX)

###### 

Click here for additional data file.

###### Known miRNA precursors not identified in this research.

(XLSX)

###### 

Click here for additional data file.

###### Novel miRNAs predicted in this research.

(XLSX)

###### 

Click here for additional data file.

###### Differentially expressed miRNAs between *91-R* and *91-C*.

(XLSX)

###### 

Click here for additional data file.

###### Target prediction of differentially expressed miRNAs between *91-R* and *91-C* with three algorithms.

(XLSX)

###### 

Click here for additional data file.

###### Gene ontology analysis of target transcripts from differentially expressed miRNAs between *91-R* and *91-C*.

(XLSX)

###### 

Click here for additional data file.

###### The pathway annotation of target transcripts from known differentially expressed miRNAs.

(XLSX)

###### 

Click here for additional data file.

###### Target genes for novel miRNAs predicted by RNAhybrid in this research.

(XLSX)

###### 

Click here for additional data file.

[^1]: **Competing Interests:**The authors have declared that no competing interests exist.
