Differentially expressed microRNAs under imidacloprid exposure and identification in Sitobion avenae

BACKGROUD MicroRNAs (miRNAs) are short single-stranded non-coding RNAs that regulate the expression of target genes, especially regulation or metabolism of endogenous or xenobiotic compounds. The de novo assembly of the transcriptomes was obtained through Illumina short-read sequencing technology in Sitobion avenae. 57 miRNAs, of which 36 were known and 21 were novel were identified. Quantitative expression levels of miRNA showed that the expression of 5 miRNAs were significant up-regulation, and the expression of 11 miRNAs were significant down-regulation in the nymph of S. avenae treated by imidacloprid compared to the control, respectively. The putative transcript target genes in S. avenae that could be regulated by these miRNAs were also carried out. The potential functions of these miRNAs in the regulation of genes involved in the metabolism, regulatory or detoxification of S. avenae were clarified based on Gene Ontology and KEGG pathway. The effects of these miRNAs identified api-miR-1000, api-miR-316, and api-miR-iab-4 on susceptibility of S. avenae to imidacloprid was determined. Modulation of the abundance of api-miR-1000, api-miR-316, and api-miR-iab-4 through the addition of inhibitors of api-miR-1000, api-miR-316, and api-miR-iab-4 to the artificial diet significantly altered the susceptibility of S. avenae to imidacloprid, which further proved that the regulatory role of these miRNAs in regulation or metabolism of insecticides. It suggested that differentially expressed microRNAs under the stress of imidacloprid could play critical regulatory role in the resistance of S. avenae to imidacloprid.


Background
Micrornas (miRNAs) are a class of endogenous small non-coding RNAs of 18~23 nucleotides in length that have been shown to be responsible for the post-transcriptional regulation of mRNAs. miRNA could suppress the translation of target mRNA molecules to silence target gene expression, through binding complementarily to 3′ untranslated regions (UTRs), coding sequences, or 5′ UTRs. [1][2][3][4][5] By imperfect complementary base pairing mainly to UTRs of target mRNAs, the miRNA guide stand directs the RISC complex, as well as to exons and5′ UTRs. 6 Though 21 to 22 bp of the guide RNA/target RNA duplex allows for loose complementarity, a near perfect match is required in the "seed" region located from positions 2 to 7 of the guide miRNA, which is considered to be related to target specifcity. 7 Studies have shown that miRNAs could regulate most of mammalian protein-coding gene, 8 which may be key mediators in a series of developmental and physiological pathways, [9][10][11][12] i.g tissue differentiation and apoptosis, cell proliferation, and morphogenesis, and the metabolism of heterogeneous compounds. 5,13 By regulating the expression of xenobiotic-metabolizing enzymes and nuclear receptors, miRNAs could mediate the detoxifcation metabolism of xenobiotics, 5 i.g human UDPglucuronosyltransferase (UGT) 1A, which participates in the metabolism of raloxifene metabolites, was proved to be post-transcriptionally regulated by mir-491-3p, 14 human P450 CYP1A1 related to the metabolism of carcinogenic and metabolites can be posttranscriptionally regulated by miR-892a. 15 Previous studies conducted in Plutella xylostella indicated that miRNAs could be involved in the resistance to chlorantraniliprole insecticides by mediating the expression levels of UGT genes. 16 UGT2B7 and UGT2B15 genes, which was proved to be post-transcriptionally regulated by mirRNA; 17 Culex pipiens suggested that miRNAs are involved in the resistance to pyrethroid insecticides by mediating the expression levels of P450 genes. 18,19 However, many known about the miRNAs are involved in regulating the detoxifcation of xenobiotics in animals and likely have essential roles in insecticide resistance, less is understood about the regulatory roles of miRNAs in the metabolism of insecticides in insects.
The grain aphid, Sitobion avenae, is a destructive pest of wheat crops distributed worldwide. The control of this pest has depended largely on the use of chemically synthesized insecticides. 20 Imidacloprid is as a chloronicotinyl insecticide that targets aphids, leafhoppers, whiteflies and other sucking pests. 21 It has been widely used for control of wheat aphid in China. Therefore, wheat aphid has evolved serious resistance to imidacloprid. 22 In this study, we performed high-throughput sequencing of short RNA libraries constructed by control and imidacloprid treatment to detect differentially expressed microRNAs, and the potential functions of these miRNAs in the regulation of genes involved in the metabolism, regulatory or detoxification of S. avenae. This will help our understanding of the regulatory roles of miRNAs in the resistance of imidacloprid in insects.

Insect culture
The strain of S. avenae has been kept in the laboratory for more 10 years, and has never been exposed to insecticides. As mentioned above, the aphids were cultured on wheat seedlings, the controlled conditions are 16-25°C, relative humidity (RH) is 60-70%, and a photoperiod is 17 h L:7 h D. 23

Imidacloprid treatments
According to the methods of Chen et al. 24 with slight modifications, leaf-dipping with aphids was performed. Imidacloprid was first dissolved in acetone, and then diluted with distilled water at an LC 30 0.5 mg L -1 concentration. Leaf-dipping with third instar nymphs aphids were immersed in the imidacloprid solution (0.5 mg L -1 ) for 5-10 s. And their imidacloprid-treated counterparts as the control (treated with distilled water) were collected at the same time points. Surviving aphids (N = 6 for each sample), including imidacloprid challenge and counterparts were collected at 24 h. The samples of aphids were flash frozen in the liquid nitrogen before preserved at -80℃. The experiment was repeated three times.

RNA isolation, miRNA library construction and Illumina sequencing
Total RNA was extracted from controls and imidacloprid treatments of S. avenae using TRIzol reagent (Invitrogen, Shanghai, China) following the manufacturer's protocol and resuspended in nuclease-free water 1% agarose gels were monitored to determine RNA degradation and contamination; RNA purity was measured using NanoDrop 2000 spectrophotometer (Thermo Scientific, Wilmington, DE); RNA concentration was checked using Qubit® RNA Assay Kit in Qubit® 2.0 Flurometer (Life Technologies, CA, USA); RNA integrity was checked using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA).
The total amount of each sample was 3 μg, which was used as input material for the small RNA library. According to the manufacturer's recommendation, the sequence library is generated by using NEBNext® Multiplex Small RNA Library Prep Set for Illumina® (NEB, USA.), and the index code is added to the attribute sequence of each sample. 3′and 5′RNA adaptors were ligated to the RNA pool using T4 RNA ligase. Then the fragments were used for reverse transcription and subsequent PCR. The final PCR products were purifed and subjected to the proprietary Solexa sequencing by synthesis method using the Illumina Genome Analyzer (San Diego, CA, USA) at the Beijing Genomics Institute (Novogene, Beijing, China).

2.
avenae miRNAs, were identified by miRDeep2 software. 25 28 The predicted miRNA target genes would be selected for further analysis. And then the target genes predicted would be aligned by the BLASTX program from NCBI (e value cut-off was 1.0E−5), as well as these gene sequences, were mapped and annotated using BLAST2GO. 29 All of the putative target genes identified were used to search the Gene Ontology (GO) database, as well as Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway for further analyses.

Comparison of miRNA expression between controls and imidacloprid treatments of S. avenae
In order to obtain the miRNAs in the response of aphids against imidacloprid, differential expression of miRNAs was performed between controls and treatments of S. avenae. The read counts of miRNA identified were analyzed using edgeR software (3.10.2), which can be obtained in Bioconductor version 3.1 (http://www.bioconductor.org/packages/release/bioc/html/edgeR.html). 30 The P-values was calculated by the Benjamini-Hochberg method. 31 P-value corrected was set at 0.05 as the threshold for significantly expressed difference by default.
Normalization of miRNA counts between libraries of controls and treatments was executed according to the total number of reads across libraries in silico. Normalized expression = Actual miRNA count/Total count of clean reads×10 6 . A false discovery rate (FDR) of the miRNAs <0.05, and a fold change ≥2 were deemed to be significant.

Quantitative real-time PCR (qPCR) validation
To validate the miRNA data obtained by the deep sequencing, 3 miRNAs were selected to confirm their expression by qPCR. Total RNA was extracted using miRNeasy Mini Kit (Qiagen, Germany) according to the manufacturer's instruction. First strand cDNA was synthesized from 2 μg of total RNA using miScript II RT kit (Qiagen) following the manufacturer's instructions. SYBR Green Master Mix (miScript SYBR Green PCR Kit, Qiagen) was used for determining miRNA expression, and performed on the ABI 7500 platform (Applied Biosystems) using an SYBR® Premix Ex Taq™ II (Tli RNaseH Plus) (Takara, Japan). And qPCR procedure was performed as previously described. 32 The optimized qPCR procedure consisted of an initial step at 94°C for 2 min, 50°C for 2 min, followed by 50 cycles of 60°C for 30 s and 94°C for 15 s followed by the melting curve.
After the cycling scheme, melting curves were got by increasing the temperature from 60 to 95°C (0.2°C s − 1 ) to denature the double-stranded DNA. The qPCR amplifications were performed in 96-well plates. Three biological replicates, with three technical replications for each, were evaluated for each sample. Analysis of the qPCR data was carried out using the 2 − ∆∆Ct method of relative quantification. 33 U6 snRNA as an endogenous control, was used to quantify miRNA expression level . The primers used for the qPCR analysis are listed in Table 1.

MiRNA inhibitor feeding in vitro, and subsequent impacts on imidacloprid susceptibility
The rearing method and the artificial diet used were as described by Gong et al. 34 with some modifications. Transparent sterilized glass tube open at both ends (4 cm in length, 2.5 cm diameter) were used as the feeding device in vitro . 25% sucrose as artificial diet was sealed between 2 layers of parafilm in a 5 cm diameter feeding device, and then healthy aphids were transfered to each device for rearing covered by the mesh to prevent them from escaping. The artificial diet was used with DEPC-treated water to ensure no RNase activity. The aphids were reared under the condition described previously. Each sample contained three replications, 100 apterous adults were used in each replication.
To evaluate the inhibitory effect of miRNAs after the feeding of miRNA inhibitors, which were mixed to the artificial diet, and the final concentration of miRNA inhibitor was 2.5 mM/L. 50 healthy apterous aphids were fed with the artificial diet, and the NC-inhibitor (a negative control) were used for the control; Following the feeding of 24 h, the aphids survival were collected for RNA extraction, 3 replicates were carried out.
To evaluate the modulation effects of the miRNAs on sensitivity of S. avenae to imidacloprid, 50 healthy apterous aphids were placed to the artificial diet that contained imidacloprid (0.1 mg/L) mixed with miRNA inhibitor at a final concentration 2.5 mM/L, the NC-inhibitor (a negative control) were used for the control. 3 replicates were carried out, and mortality was caculated at 48 h. The miRNA inhibitors used were provided by Sangon Biotech Co., Ltd (Shanghai, China).

Illumina sequencing data analysis
There were 14, 841,584, and 13, 964, 533 raw reads, identified from the control and imidacloprid treatment RNA libraries, respectively (  Table 2). The length distribution of unique filtered reads indicated a distribution in the 2 libraries (Fig. 1B), where constituted 3.75% and 3.72% of the total reads at 22 nt, respectively, where at 23 nt constituted 5.18% and 6.23% of the total reads, respectively. Where in the peak at 32 nt constituted 8.59% of the total reads of imidacloprid libraries, at 34 nt constituted 9.76% of the total reads of controls, which probably showed piRNA-like sRNAs in S. avenae. A little bit sRNAs had over a thousand reads, however most had fewer than 10 copies in the library (Fig. 1C). sRNAs in S. avenae showed a significant bias for the nucleotide U at 5′ (96.99%, and 97.51%, respectively) and 3′ ends (53.38%, and 55.32%, respectively) in both libraries (Fig. 2).

Sequence analysis of miRNA and predication.
According to the length distribution of the reads, miRNAs at 22 nt accounted for 43.86% of all miRNA species of sRNA libraries (Fig. 1D). The query of these reads of putative S.
avenae miRNAs to A. pisum genomes and all insect miRNAs in miRBase, the results showed 36 conserved sequences were identified (Table S1). In addition, the query of S. avenae filtered reads against mature miRNAs in the miRBase database, the results showed that 26 miRNA sequence families were identified (Fig. 3A). It is predicted that the miR-10 family has the largest members (n = 415), followed by let-7 (n = 401), and miR-2 (n = 180).
Insects were shared with these miRNAs, with the majority of miRNAs shared with the known A. pisum miRNAs, and a fewer extent trend for orthologs of other insects (Fig. 3B).
There were a total of 21 novel miRNAs potentially obtained, and named the prefx 'PC' (predicted candidate) based on the adopted nomenclature (Table S2). These specifc miRNAs could be mapped to S. avenae transcriptome sequences. It was predicted from Mfold that stem-loop hairpin secondary structures (<18 kcal/mole 35 ) were consist of miRNAs precursor sequences, and some predicted structures in S. avenae were shown in Fig. S1.
Differentially expressed miRNAs between imidacloprid treatments and controls were also determined according to normalized differences of Illumina read counts, 57 S. avenae miRNAs, 16 of which (28.1%) were differential expression, including 11 down-regulation and 5 up-regulation ( Table 3).
The PCR amplified products of api-miR-1000, api-miR-316, and api-miR-iab-4 indicated a single band as the expected size (60-100 bp). Subsequent analysis of qPCR results confirmed that the expressed level of Api-miR-1000 was 5.5-fold higher in imidacloprid treatments compared to controls. While the expression of Api-miR-316, and Api-miR-iab-4 were 0.09-, and 0.12-fold lower in imidacloprid treatments compared to controls, respectively. (Fig. S2). Overall, qPCR results of miRNAs were indicated similar ranging trends as that of the RNA-seq analyses.

Target prediction of miRNAs.
The 16 differentially expressed miRNAs in S. avenae, the hypothetical targets within the UTRs of transcripts from S. avenae for the miRNAs were predicted using Target Scan software (Table S3). Possible functions predicted of these hypothetical target genes were annotated by GO enrichment, which could be involved in molecular functions, biological processes and cellular components, like development processes (Table S4). There were 58 GO terms identified for target genes of miRNAs predicted according to GO level 2 ( Fig.   4A), like the genes participated in metabolic processes.
These genes including 71 metabolic pathways predicted by KEGG enrichment analysis might be affected by miRNA regulation (targeting) (Fig. 4B). These transcript targets were predicted to orthologs of related insects, which showed that differentially expressed miRNAs between the control and the imidacloprid treatment in S. avenae could bind and promote the post-transcriptional regulation, including various biological pathways like xenobiotic metabolism (Table S5).
To investigate the potential roles of miRNAs in the resistance of S. avenae to imidacloprid, the predicted target genes for the miRNAs that could participate in xenobiotic metabolism were further focused. Interestingly, target genes for some miRNAs were found to play critical roles in insects response to xenobiotic stress, including cytochrome P450s, UDPglucuronosyltransferase, glutathione S-transferase, etc. (Table 4). Some miRNAs predicted have a lot of putative target genes, and many putative target genes were regulated by multiple miRNAs. (Table S3).

Modulation of miRNAs consequently impacts the susceptibility of S. avenae to imidacloprid
The expression levels of miRNAs, api-miR-1000, api-miR-316, and api-miR-iab-4, was significantly depressed in S. avenae aphids fed miRNA inhibitors compared to the levels in S. avenae aphids fed DEPC water and NC-inhibitor as controls (Fig. 5A) inhibitors under imidacloprid exposure compared to the controls, respectively (Fig. 5B).

Discussion
It is predicted miRNAs can regulate a lot of protein-coding genes. [36][37][38] And it is demonstrated that miRNAs play critical regulatory roles in a lot of biological processes. In D. melanogaster dme-miR-1000 was located within the intron of musashi, 41 encoding an RNA-binding protein to regulate gene translation by prior expression in the nervous system; 42 Honeybee miR-1000 could be involved in musashi function of the nervous system, 43 ame-miR-1000 was as the first identified miRNA, which was expressed in a brain-selective manner of an insect. 1 cytochrome P450 and 2 glutathione S-transferase genes could be recognized by Mse-miR-316; 44 The accumulation endogenous Ubx protein could be attenuated as ectopic expression of miR-iab-4-5p, and a classical homeotic mutant phenotype could also be induced (the transformation of halteres into wings), 45 Drosophila dme-miR-iab-4 and dme-miR-iab-8 as bithorax-complex (BX-C) HOX clusters miRNAs could regulate neural patterning and reproduction through restricting HOX gene targets 46 .
In order to further understand the effect of miRNAs on the regulation of metabolism of imidacloprid, the hypothetical target genes for miRNAs identified were predicted, and many of them annotated could participate in a lot of biological processes. At the same time, several target genes from families for miRNAs were known to be involved in the metabolism of xenobiotics. 47,48 We found that several of the miRNAs targeting the genes related to xenobiotic metabolism were among the differentially expressed miRNAs, 49, 50 by combining the results of the target prediction with our differential expression analysis of miRNAs.
The mortality of aphids under imidacloprid exposure decreased significantly after their feeding api-miR-1000 inhibitor. While the mortality of aphids under imidacloprid exposure increased significantly after their feeding api-miR-316, and api-miR-iab-4 inhibitors, respectively. This showed that miRNAs differential expressied may participate in metabolism of imidacloprid, by regulating the expression of xenobitic metabolism genes.
In participated in the data analysis. All authors reviewed the manuscript.

Competing interest
The authors declare that they have no competing interests.   Note: 3ADT: reads that were removed as 3ADT was not found and length was less than 15nt and was more than 35nt. Mappable reads: reads from the raw reads that were passed by many the digital filters .  Figure 1 Characterization of miRNA sequences from Sitobion avenae Figure 2 Nucleotide bias of predicated miRNAs from Sitobion avenae Figure 3 MiRNAs characterization of Sitobion avenae.

Figure 4
Histogram presentation of GO annoation and KEGG pathway for miRNAs identified from Sitobion avenae.

Figure 5
Effects of miRNA modulation on susceptibility of Sitobion avenae to imidacloprid.

Supplementary Files
This is a list of supplementary files associated with the primary manuscript. Click to download.