Identi cation of differentially expressed microRNAs under imidacloprid exposure in Sitobion avenae Fabricius

Baizhong Zhang Henan Institute of Science and Technology Xu Su Henan Institute of Science and Technology Lanfen Xie Henan Institute of Science and Technology Wen-Yang Dong China Agricultural University Fan-Bin Kong Henan Institute of Science and Technology Junjie Liu Henan Instiute of Science and Technology Liuyang Lu Henan Institute of Science and Technology Ya-She Li Henan Institute of Science and Technology Shouping Zhang Henan Institute of Science and Technology Ming-Wang Shi Henan Institute of Science and Technology Xi-Ling Chen (  chenxiling456@126.com ) Henan Institute of Science and Technology Xi-Wu Gao China Agricultural University


Background
MicroRNAs (miRNAs) are between 18 and 23 nucleotides in length and have been shown to be responsible for the post-transcriptional regulation of mRNAs.2][3][4][5] UTRs of target mRNAs are mainly bound through imperfect complementary base pairing. 6The sequence of the target is an exact match to the miRNA. 7MiRNAs can mediate the detoxi cation/metabolism of xenobiotics through regulation of the expression of detoxi cation genes, including nuclear receptors and xenobiotic-metabolizing enzymes. 5For example, human UDP-glucuronosyltransferase (UGT) 1A, which metabolizes raloxifene metabolites, was shown to be regulated post-transcriptionally by mir-491-3p, 14 and human P450 CYP1A1, which is associated with the metabolism of carcinogenic and metabolites, can be regulated post-transcriptionally by miR-892a. 15revious studies conducted in Plutella xylostella indicated that miRNAs could be involved in chlorantraniliprole resistance by regulating the expression levels of UGT genes. 16UGT2B7 and UGT2B15 were shown to be regulated post-transcriptionally by miRNAs; 17 Culex pipiens miRNAs are involved in pyrethroid resistance by regulating the expression levels of P450 genes. 18,19 here are many miRNAs known to be involved in regulating enzymes that metabolize xenobiotics in animals, while the effects of miRNAs on the regulation of xenobiotic detoxi cation/metabolism in wheat aphids is less understood.
The grain aphid Sitobion avenae is the main pest of wheat and is distributed worldwide.S. avenae not only harms crops by sucking juice from plants but also spreads yellow dwarf virus disease, leading to serious economic losses.The control of this aphid has depended largely on the use of chemically synthesized pesticides, 20 such as imidacloprid, which is a neonicotinoid that targets sucking pests, such as aphids and the green leaf bug. 21It has been widely applied in wheat elds to control wheat aphid in China.As a result, S. avenae has developed serious resistance to imidacloprid. 22re, we performed high-throughput sequencing of short RNA libraries of control and imidacloprid-treated S. avenae nymphs to detect differentially expressed microRNAs, and determine the functions of miRNAs in regulating genes involved in the metabolism/detoxi cation of xenobiotics.The results of this study will help increase our understanding of the roles of miRNAs in the regulation of resistance to imidacloprid in insects.

Culture conditions
The susceptible strain of S. avenae used was maintained in a greenhouse for more than 10 years and has never been exposed to any insecticides.The culture conditions used were described previously. 23

Imidacloprid treatments
To test the susceptibility of S. avenae to imidacloprid, the leaf-dipping method was conducted as previously described by Chen et al. 24 with minor modi cations.Imidacloprid was dissolved in acetone, and then serially diluted with 0.05% (v/v) Triton X-100 in water.Wheat leaves were cut into pieces 20 mm long.The leaves with aphids were dipped into 0.05% (v/v) Triton X-100 in water with 0.5 mg/L imidacloprid (LC 10 concentration, Table S1) or without imidacloprid as a control.The leaves with treated aphids were put into a glass tube (6 cm in length, 2 cm diameter) and the open end was covered with cotton to prevent insect escape.Surviving third instar nymphs (N = 5 for each treatment) were collected at 24 h after imidacloprid or mock challenge.Aphids were ash frozen in liquid nitrogen before being stored at -80℃ for RNA extraction.The experiment was replicated three times.

RNA isolation, miRNA library construction and Illumina sequencing
Total RNA was extracted from control and imidacloprid-treated S. avenae using TRIzol reagent (Invitrogen, Shanghai, China) according to the manufacturer's protocol and resuspended in nuclease-free water.Electrophoresis on 1% agarose gels was used to determine if the extracted RNA was degraded or contaminated, and then purity, concentration and integrity of the RNA were measured using a NanoDrop 2000 spectrophotometer (Thermo Scienti c, Wilmington, DE), a Qubit® RNA Assay Kit with a Qubit® 2.0 Fluorometer (Life Technologies, CA, USA), and a RNA Nano 6000 Assay Kit with the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA), respectively.Approximately 3 μg of total RNA from each sample was used as input material for constructing the small RNA library, and the platform used for sequencing was conducted in Novogene, Beijing, China as previously described. 254 Bioinformatics analysis, miRNA prediction, and miRNA target prediction.S. avenae miRNAs, were con rmed using miRDeep2 software.26 Raw reads from the four libraries (imidacloprid treatment and control, two biological replicates each) were used as input for miRDeep2, and the data from each library were separately analyzed.The default options and settings were used to perform the miRDeep2 analysis.The raw read sequences with polyA tails and miRNAs with lengths ranging from 18 to 30 nt were selected for further analysis after trimming adaptor sequences and removing snRNA, rRNA, snoRNA and tRNA sequences.The sequences that mapped to the predicted mature Acythosiphon pisum miRNAs annotated in miRBase were identi ed as known mature miRNAs.The clean reads from the S. avenae transcriptome have been deposited in the NCBI/SRA database under accession number SRP182781. Th Miranda, RNAhybrid, and Target Scan programs were used to predict the target gene of the miRNAs according to the miRNA-binding sites.[27][28][29] Then only the binding sites commonly identi ed by the three tools were selected for subsequent analysis.The putative miRNA targets were used as queries in searches against the S. avenae transcriptome sequences.
The predicted miRNA target genes were selected for further analysis.The predicted target genes were aligned using the NCBI BLASTX program, and mapping and annotation of the gene sequences were conducted using BLAST2GO. 30The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases were used to further analyze the putative target genes of the miRNAs.

miRNAs differentially expressed between control and imidacloprid-treated S. avenae
In order to identify the aphid miRNAs that were affected by imidacloprid, the differences in the expression levels of miRNAs, as assessed by miRDeep2, between control and imidacloprid-treated S. avenae were determined.The edgeR software package (3.10.2) obtained from Bioconductor version 3.1 was used to analyze the read counts of the identi ed miRNAs (http://www.bioconductor.org/packages/release/bioc/html/edgeR.html). 31edgeR is a Bioconductor software package for analyzing the differential expression of digital gene expression data.The P-values were calculated using the Benjamini-Hochberg method. 32A corrected P-value of 0.05 was set as the default threshold for signi cance.Normalization of miRNA counts between libraries of the controls and treatments was executed according to the total number of reads across libraries.Normalized expression = Actual miRNA count/Total count of clean reads×10 6 .A fold change in miRNA expression ≥2 with a false discovery rate (FDR) <0.05 was deemed to be signi cant.

Quantitative real-time PCR (qPCR) validation
To validate the miRNA data obtained by deep sequencing, the expression levels of three miRNAs were con rmed by qPCR.First strand cDNA was synthesized from 2 μg of total RNA using the miScript II RT kit (Qiagen) following the manufacturer's instructions.The SYBR Green Master Mix (miScript SYBR Green PCR Kit, Qiagen) was used for miRNA expression assays, and qPCR was performed as previously described. 33Relative expression was calculated using the 2 −∆∆Ct method. 34U6 snRNA was used to normalize the expression levels of miRNAs as an endogenous control.The primers used are shown in Table 1.

MiRNA inhibitor feeding in vitro, and the subsequent impacts on imidacloprid susceptibility
The rearing method and the arti cial diet used were the same as previously described. 35A sterilized transparent glass tube with both ends open was used as the feeding device (4 cm in length, 2.5 cm diameter).The arti cial diet, 25% sucrose, was put between the two layers of para lm were sealed or was the para lm containing the arti cial diet used to seal one end of the feeding device, and then healthy apterous adults were transferred to the device; the ends were covered with mesh to prevent the insects from escaping.Each sample contained three replications.
To evaluate the inhibitory effect of the miRNAs on CYP6A14, 50 healthy apterous aphids were fed an arti cial diet containing miRNA inhibitors at a nal concentration of 2.5 mM/L.The NC-inhibitor (a negative control) was used as the control.Following feeding for 24 h, the surviving aphids were collected for RNA extraction.There were three replicates for each sample.
To evaluate the effects of miRNAs on the sensitivity of S. avenae to imidacloprid, 50 healthy apterous adults were fed an arti cial diet that contained imidacloprid (at the LC 50 concentration, 1.5 mg/L) with an miRNA inhibitor (2.5 mM/L).The NC-inhibitor was used as the control.Three replicates were carried out, and mortality was scored at 48 h.The miRNA inhibitors used were provided by Sangon Biotech Co., Ltd (Shanghai, China).

Illumina sequencing data analysis
A total of 14,841,584 raw reads were obtained from sequencing the control library and 13,964,533 raw reads were obtained from the imidacloprid treatment library (Table 2).The lengths of the S. avenae sRNAs ranged between 18 nt and 35 nt, comprising 5.28% of the imidacloprid reads and 5.08% of the control reads were 22 nt in size (Fig. 1A).After discarding low quality reads and adapter sequences, 6,992,451 (47.11% of total) and 6,349,363 (45.47%) unique mappable reads remained from the control and imidacloprid treatment libraries, respectively (Table 2).The length distribution of these mappable reads indicated that 3.75% of the total reads for the control and imidacloprid treatment, respectively, were 22 nt in length and 5.18% and 6.23% of the total reads, respectively, were 23 nt in length.The maximum lengths of reads from the imidacloprid and control libraries were 32 nt (8.59% of all reads) and 34 nt (9.76%), respectively, these reads are likely piRNA-like sRNAs in S. avenae.A few sRNAs had over a thousand reads; however, most had fewer than 10 copies (Fig. 1C).
The dominance of uracil (U) at the rst position of the 5′ terminus is considered to be a conserved feature of mature miRNAs.Given that the rst base of the 5′ end of miRNAs is known to play very important roles in the interaction between miRNAs and Argonaute complexes, the position-speci c occurrence of nucleotides in the candidate miRNA sequences was analyzed.S. avenae miRNAs showed a bias towards U at the rst nucleotide position; sRNAs showed a signi cant bias for U at the 5′ (96.99%, and 97.51%, respectively) and 3′ ends (53.38%, and 55.32%, respectively) in both libraries.The nucleotide U was the most abundant nucleotide at most of the positions.This higher abundance of U was signi cant at positions 1, 6, 9, 12, 16, and 22 based on analysis of the base composition of the miRNAs at each position.

Sequence and prediction of miRNAs
According to the distribution of read lengths, 22-nt miRNAs accounted for 43.86% of all species of miRNAs (Fig. 1D).The putative miRNAs in S. avenae were used as queries in searches against the miRNAs from the A. pisum genome in miRBase, and 36 conserved sequences were identi ed (Table S2).In addition, using the ltered S. avenae reads as queries in searches against mature miRNAs identi ed 26 miRNA sequence families from the miRBase database (Fig. 3A).It is predicted that the miR-10 family has the largest number of miRNAs (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, while a fewer extent miRNAs trended for other insect orthologs (Fig. 3B).
A total of 21 potentially novel miRNAs were identi ed; these were named with the pre x 'PC' (predicted candidate) based on the adopted nomenclature (Table S3).These novel miRNAs could be mapped to S. avenae transcriptome sequences.Mfold predictions revealed stem-loop hairpin secondary structures (<18 kcal/mole 36 ) that were consistent with miRNA precursor sequences.Some of the predicted structures are shown in Fig. S1.Differential expression of miRNAs between the imidacloprid treatments and controls was also determined according to normalized differences between Illumina read counts.Of the 57 S. avenae miRNAs, 16 (28.1%)were differentially expressed, 11 of which were down-regulated and 5 of which were up-regulated (Table 3).
The PCR ampli cation products of api-miR-1000, api-miR-316, and api-miR-iab-4 were single bands with the expected sizes between 60 and 100 bp.Subsequent qPCR analysis con rmed that the expression level of api-miR-1000 was 5.5-fold higher in the imidacloprid treatment samples than in the controls.The expression levels of api-miR-316, and api-miR-iab-4 were 91%, and 88% lower, respectively, in the imidacloprid treatment samples compared with the controls (Fig. S2).Overall, the results of qPCR analysis were similar to those of RNA-seq.

Prediction of miRNA targets.
The hypothetical miRNA target sequences within the UTRs of transcripts from S. avenae (including the targets of the 16 differentially expressed miRNAs) were predicted using Target Scan software (Table S4).To identify the possible functions of these hypothetical target genes, GO enrichment analysis was performed; target genes were assigned molecular function, biological process and cellular component terms (Table S5).The predicted miRNA target genes were annotated to 58 GO level 2 terms(Fig.4A), like some of these terms are related to metabolic processes.
In KEGG enrichment analysis miRNA target genes were assigned to 71 metabolic pathways, indicating that these pathways might be targeted by miRNA regulation (Fig. 4B).The orthologs of some miRNA targets are involved in biological pathways such as xenobiotic metabolism, indicating that miRNAs might be involved in the posttranscriptional regulation of these pathways (Table S6).
To con rm the potential function of miRNAs in the resistance of S. avenae to imidacloprid, the predicted target genes of the miRNAs that might participate in xenobiotic metabolism were further focused on.Interestingly, target genes for some miRNAs, such as those encoding cytochrome P450s, UDPglucuronosyltransferase and glutathione S-transferase, have been found to play critical roles in insect responses to xenobiotics (Table 4).Some miRNAs have many putative target genes, and many putative target genes are potentially regulated by multiple miRNAs (Table S3).

Modulation of miRNAs impacts the expression of predicted target genes and the susceptibility of S. avenae to imidacloprid
The expression levels of miRNAs miR-1000, miR-316, and miR-ab-4 were signi cantly lower in S. avenae adults fed an arti cial diet containing the corresponding miRNA inhibitors than in adults fed the DEPC water control or NC control, The depression e ciencies of miR-1000, miR-316, and miR-ab-4 reached 85.0%, 89.1%, and 88.2%, respectively (Fig. 5A).The transcriptional levels of CYP6A14 were dramatically increased by 2.48-fold in the miR-306-knockdown aphids.However, inhibition of the other two miRNAs had no signi cant effect on CYP6A14 transcript levels (Fig. 5B).Under imidacloprid stress, mortality decreased by 21.73% in the aphids fed miR-1000 inhibitors compared with those fed the controls.In contrast, mortality increased by 26.75% and 21.86% in the aphids fed miR-316 and miR-ab-4 inhibitors, respectively, compared with those fed the controls (Fig. 5C).

Discussion
8][39] It has also been demonstrated that miRNAs are involved in the regulation of various biological processes.In this study, miRNAs in S. avenae that might play regulatory roles in target gene responses to imidacloprid were identi ed and assessed.The identi cation and functional analysis of miRNAs differentially expressed in imidacloprid-treated S. avenae were conducted to increase our understanding of the regulation/detoxi cation of imidacloprid in insects.So far, high-throughput sequencing of small RNAs has been widely used to obtain miRNAs from a series of organisms, including S. avenae.Li et al. 40 identi ed 345 miRNAs from winged and wingless S. avenae, whereas we only identi ed 57 miRNAs in control and imidacloprid-treated S. avenae as we only selected Acyrthosiphon pisum as a reference.However, our sequencing data greatly enlarges the range of available information about the involvement of miRNAs in the regulation of metabolism genes and provide a critical basis for further studying miRNAs related to insecticide resistance in S. avenae.
Our analysis of read counts indicated that 16 miRNAs were differentially expressed between imidaclopridtreated and control S. avenae, indicating that imidacloprid could modulate the expression of miRNAs and suggesting a probable regulatory role of these miRNAs in the detoxication/metabolism of insecticides in S. avenae.Several miRNAs up/down-regulated by imidacloprid treatment, e.g., api-miR-1000, api-miR-316, and api-miR-iab-4, could be involved in the regulation of imidacloprid metabolism.The activity-dependent expression of Drosophila melanogaster miR-1000 could indicate a mechanism for allowing neuronal activity to ne-tune the strength of excitatory synaptic transmission, and regulate the expression level of vesicular glutamate transporters, and load the glutamate into synaptic vesicles. 41A target site of D. melanogaster dme-miR-1000 is located in the intron of musashi, 42 which encodes an RNA-binding protein, that regulates gene translation by prior expression in the nervous system. 43Honeybee miR-1000 is also thought to be involved in musashi function in the nervous system. 44The rst miR-1000 identi ed in insects, and is expressed in a highly selective manner in insect brains.One cytochrome P450 and two glutathione S-transferase genes could be recognized by mse-miR-316. 45The accumulation of endogenous Ubx protein could be attenuated by ectopic expression of miR-iab-4-5p, and a classical homeotic mutant phenotype could also be induced (the transformation of halteres into wings). 46rosophila dme-miR-iab-4 and dme-miR-iab-8 are miRNAs targeting bithorax-complex (BX-C) HOX bunches and participate in the regulation of reproduction and neural patterning through suppression of HOX gene targets. 47 order to further explore the role of miRNAs in the regulation of imidacloprid metabolism, the hypothetical target genes for the miRNAs identi ed were predicted, and many of them were annotated to biological processes.The target genes of several miRNAs from different families were previously shown to be involved in xenobiotic detoxication/metabolism. [48][49][50][51] By combining differential analysis with target gene prediction, we also found that several of the miRNAs targeting genes putatively involved in xenobiotic detoxication/metabolism were expressed differentially between the imidacloprid-treated and control S. avenae.S. avenae CYP6A14 could be regulated by api-miR-316 in our present study.Furthermore, the mortality of S. avenae exposed to imidacloprid signi cantly decreased after being fed an api-miR-1000 inhibitor.In contrast, the mortality of S. avenae exposed to imidacloprid signi cantly increased after being fed an api-miR-316 or api-miR-iab-4 inhibitor.The differential expression of let-7 and miR-100 might be attributable to the in uence of insecticides.A similar phenomenon was also observed in the host adaption of Myzus persicae to nicotianae, where let-7 and miR-100 participated in the posttranscriptional regulation of the expression of CYP6CY3.Our results show shows that differentially expressed miRNAs may participate in the metabolism of imidacloprid by modulating the expression of genes related to xenobiotic metabolism.
In conclusion, of 57 miRNAs identi ed in S. avenae, the transcript levels of 16 were differentially regulated by imidacloprid treatment.The prediction of target genes related to xenobiotic metabolism showed that differentially expressed miRNAs could participate in the detoxication/metabolism of imidacloprid.The data presented here represent a crucial new genomics resource for further study of small RNAs in S. avenae.

Figure 5 Effects
Figure 5

Table 1 .
Primers used in the study

Table 2 .
Distribution of miRNA reads from control imidacloprid-treated Sitobion avenae

Table 3 .
The differentially expressed miRNAs in the small RNA libraries of Sitobion avenae

Table 4
MiRNAs targeting putative genes related to xenobiotic metabolism in Sitobion avenae