3.1 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 first position of the 5′ terminus is considered to be a conserved feature of mature miRNAs. Given that the first base of the 5′ end of miRNAs is known to play very important roles in the interaction between miRNAs and Argonaute complexes, the position-specific occurrence of nucleotides in the candidate miRNA sequences was analyzed. S. avenae miRNAs showed a bias towards U at the first nucleotide position; sRNAs showed a significant 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 significant at positions 1, 6, 9, 12, 16, and 22 based on analysis of the base composition of the miRNAs at each position.
3.2 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 identified (Table S2). In addition, using the filtered S. avenae reads as queries in searches against mature miRNAs identified 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 identified; these were named with the prefix ‘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/mole36) 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 amplification 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 confirmed 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.
3.3 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 confirm 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, UDP-glucuronosyltransferase 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).
3.4 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 significantly lower in S. avenae adults fed an artificial diet containing the corresponding miRNA inhibitors than in adults fed the DEPC water control or NC control, The depression efficiencies 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 significant 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).