The susceptible population of P. xylostella was maintained under insecticide free conditions for 10 generations in the Engineering Research Centre of Biological Control, Ministry of Education, South China Agricultural University (SCAU). The insects were kept at 60–70% relative humidity and at 25 ± 1°C under a 16:8 h light: dark cycle.

Stain IfB01 of I. fumosorosea (China Center for Type Culture Collection access number: CCTCC M 2012400) was grown on potato dextrose agar (PDA) at 26°C. The conidia were prepared as described previously (Huang et al., 2010). Healthy third instar larvae of P. xylostella were selected and treated with 1 × 10⁷ spores/ml suspension and then surviving larvae (50) were collected at 12, 18, 24, and 36 h, post-treatment. The control group larvae were treated with sterile deionized water containing 0.05% Tween-80 and the samples were collected at 0 h post-treatment.

Trizol Total RNA Isolation Kit (Takara, Japan) was used to extract total RNA from normal and treated larval samples following manufacturer's instructions. The concentrations of RNA were assessed using Nanodrop (Bio-Rad, USA) and its integrity was determined on Agilent 2100 Bioanalyzer (Agilent, USA). The small RNA libraries were constructed from each time-point of infection using a TruSeq small RNA sample preparation kit (Illumina). Briefly, RNAs were firstly ligated with 3′ adapter and after size fraction ligated to 5′ adapter. The small RNA fractions were then used for reverse transcription following PCR. The final ligation PCR products, after purification, were sequenced using Illumina Genome Analyzer (San Diego, CA, USA) at the Beijing Genomics Institute (BGI, Shenzhen, China).

To screen clean reads, raw data reads were filtered to remove low-quality, 5′ primer contaminants, without 3′ primers and insert tag, and sequences fewer than 18 nucleotides. The remaining high-quality reads were initially mapped to P. xylostella genome (GCA_000330985.1) using Bowtie software (Langmead and Salzberg, 2012), and then annotated into different classes to remove rRNA, scRNA, snoRNA, snRNA, and tRNA using Rfam database. Finally, the unannotated clean sequences were used to predict novel miRNAs using the miRDeep2 software.

The potential mRNA targets of differentially expressed miRNAs were predicted and analyzed using three different programs, such as RNAhybrid, miRanda, and TargetScan following already established criteria for target prediction (Allen et al., 2005; Schwab et al., 2005). To get more reliable results, we selected those mRNA targets which were predicted by all three programs. Additionally, functional annotation of all the predicted target genes was conducted by using Gene Ontology (GO) database and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses, with the threshold set at a corrected P-value ≤ 0.05.

Real-time quantitative PCR (RT-qPCR) is the method of choice for analyzing expression of genes and to confirm the results of RNA-Sequencing (Shakeel et al., 2017b). Thus, to confirm the results of sRNA-Seq in the current study, RT-qPCR analysis was conducted to ensure the expression levels of miRNAs displayed by Illumina sequencing results and 10 miRNAs were selected. RT-qPCR was performed on a Bio-Rad iQ2 optical system (Bio-Rad) using SsoFast EvaGreen Supermix (Bio-Rad, Hercules, CA, USA) following the instructions of the manufacturer. The U6 snRNA was used as an internal control. The reaction program was set as 95°C for 30 s, 40 cycles of 95°C for 5 s, and 55°C for 10 s with a dissociation curve generated from 65 to 95°C to ensure the purity of PCR products (Shakeel et al., 2015). Each experiment was replicated in triplicate. Finally, data analysis was performed using 2^−ΔΔCT method (Livak and Schmittgen, 2001).

To identify miRNAs in I. fumosorosea challenged P. xylostella, we constructed five small RNA libraries (Tween (TW), 12, 18, 24, and 36 h) using high-throughput Illumina sequencing platform. In total, 11,861,547; 11,872,699; 11,944,980; 11,956,814, and 11,866,077 raw reads were obtained, respectively. After low-quality sequences, adaptors, and sequences less than 18 nucleotides were discarded, 92.93, 98.17, 98.96, 98.63, and 94.62% clean reads were obtained in TW, 12, 18, 24, and 36 h, respectively, for further analysis (Table 1).

The small RNA size distribution in the five libraries showed that most of the sRNAs ranged from 18 to 30 nt, with 28 nt being the most abundant following 22, and 23 nt (Figure 1) in the five libraries. The two peaks observed at 22 and 28 nt, in the present study, represent a typical length of miRNAs and piwi-interacting RNAs, respectively. Our findings are in consistency with the typical size of miRNAs and piwi-interacting RNAs in previous reports (Wei et al., 2009; Etebari et al., 2013; Xu et al., 2015; Li et al., 2016). Among the clean reads, 85.10% sRNAs were common between 12 and 18 h, 85.31% sRNAs between 24 and 12 h, 83.28% sRNAs between 24 and 18 h, 83.92% sRNAs between 24 and 36 h, 84.55% sRNAs between 24 h and TW, 84.89% sRNAs between 36 and 12 h, 83.91% sRNAs between 36 and 18 h, 86.44% sRNAs between 36 h and TW, 85.77% sRNAs between TW and 12 h, and 84.18% sRNAs between TW and 18 h, respectively (Supplementary Figure 1).

Of the clean reads, 6,784,087, 6,965,570, 7,410,628, 7,023,256, and 6,510,880 reads from control, 12, 18, 24, and 36 h accounted for 61.54, 59.76, 62.69, 59.55, and 57.99%, respectively, and were mapped to the genome of P. xylostella (Supplementary Table 1). The annotation of sRNAs was carried out by following priority rule of rRNA etc.; (GenBank>Rfam) > known miRNA > repeat > exon > intron (Calabrese et al., 2007). The clean reads were categorized into miRNA, rRNA, snRNA, snoRNA, tRNA, and unannotated (unann). The composition of the sRNA classes in each library is displayed in Supplementary Figure 2.

After successful mapping of clean reads against P. xylostella genome, the mapped miRNA sequences were matched to miRNAs reported by Etebari and Asgari (2016). Our analysis initially identified, based on sequence similarity, in total, 191 mature miRNAs. Then, precursor sequences of these mature miRNAs were aligned to those reported by Etebari and Asgari (2016), and 102 highly confident precursor miRNAs, which produced 172 of 194 mature miRNAs. Our analysis indicated that precursor miRNA sequences of the remaining 22 conserved miRNAs were not detectable in the current assembly of P. xylostella genome. After removing those known miRNAs with read count < 10 in all libraries, remaining 116 known miRNAs with precursor sequences (Supplementary Table 2), and 15 miRNAs without precursor sequences (Supplementary Table 3) were retained for further analysis. The remaining sequences that were not matched to conserved miRNAs were used to predict novel miRNAs by using miRDeep2 program (Friedländer et al., 2012). The prediction of novel miRNAs analysis predicted 42 potential novel miRNAs from all the libraries (Supplementary Table 4) following the standard criteria of novel miRNA prediction with a miRDeep score >1, randfold P-value < 0.05, and MFE < −19 kcal/mol.

It is worth mentioning that a low copy number of miR-1, a conserved miRNA, was detected after parasitization in a previous report (Etebari et al., 2013), however, in the present study, miR-1 was the most abundant miRNA following pxy-let-7-5p, pxy-miR-184-3p, pxy-miR-10-3p, and miR-31-5p (Table 2). The abundant and common expression of these conserved miRNAs indicates that these miRNAs might play crucial roles in P. xylostella. Our results are in consistency with previous reports where a high expression of these miRNAs was observed in other insect small RNA libraries (Cai et al., 2010; Cristino et al., 2011; Liu et al., 2012). Bantam, a most abundantly expressed miRNA, plays multiple roles in insects such as apoptosis inhibition, cell proliferation and stem cell stem cell maintenance, and immunity in Drosophila melanogaster (Smibert and Lai, 2010; Fullaondo and Lee, 2012). Although a high copy number of bantam was observed in our study, however, its up-regulation after infection was less than 1-fold. Let-7, a highly conserved miRNA, has also been reported to play an important role in immunity, for example, it binds to 3' UTR of antimicrobial peptide diptericin to repress translation of this protein in D. melanogaster (Garbuzov and Tatar, 2010). Interestingly, we found that few miRNAs like miR-2755, miR-10, and miR-31 showed high expression in all treatments. A higher expression of these miRNAs after fungal treatment indicates that these miRNAs might play important roles in defending P. xylostella against pathogens.

The novel miRNA analysis identified 42 potential novel miRNAs in P. xylostella after infection (Supplementary Table 4). Among novel miRNAs, pxy-novel-miR-26 was the abundantly expressed miRNA following pxy-novel-miR-1, pxy-novel-miR-33, and pxy-novel-miR-35 (Supplementary Table 4).

The differential abundance of host miRNAs, a common observation in host-pathogen systems, changes at different infection stages following an infection (Asgari, 2011). In the present study, to find out the I. fumosorosea responsive miRNAs, a differential expression analysis was performed using the sequencing results (Figure 2). The differential expression analysis exhibited that 13, 12, 16, and 5 known miRNAs were differentially expressed in 12, 18, 24, and 36 h, respectively, compared to control (Supplementary Table 5). The top five differentially expressed known miRNAs are presented in Table 3. Furthermore, 12, 19, 13, and 11 novel miRNAs were differentially expressed, in 12, 18, 24, and 36 h, respectively, compared to control (Supplementary Table 6).

Interestingly, in the present study, we found that the expression of few conserved miRNAs like miR-2, miR-9, miR-279, miR-745, miR-7b, and miR-2767 was changed following the infection of I. fumosorosea. Our findings suggest that these miRNAs might play very important roles in P. xylostella immunity to I. fumosorosea. In accordance to our study, previous reports also suggested the important roles of these conserved miRNAs in the immunity of different insects against different pathogens, such as bacteria-injected larvae of Manduca sexta and Diadegma semiclausum parasitized P. xylostella resulted in differential expression of miR-2, miR-9, and miR-279, indicating the role of these miRNAs in immunity against bacteria and parasite, respectively (Zhang et al., 2012; Etebari et al., 2013). It is of note that miR-9 has been predicted to play an essential role in signal recognition in M. sexta, and in toll pathway in Drosophila melanogaster (Fullaondo and Lee, 2012). It is worth mentioning that the read number of most of the miRNAs dropped after infection and, overall, only 3 miRNAs (miR-282,−2796, and −34) were up-regulated while 20 miRNAs were down-regulated in all the treatments compared to control. Previously, it has been reported that when Galleria mellonella was infected with entomopathogenic fungi, Metarhizium anisopliae, at larval stage, only one miRNA (miR-210b) showed differential expression (Mukherjee and Vilcinskas, 2014), whereas, in our study, 23 miRNAs were differentially expressed, however, miR-210b was not detected in our small RNA libraries.

To validate small RNA sequencing results, 10 randomly selected miRNAs were analyzed by RT-qPCR (Figure 3). The results exhibited that the trend of the expression level of the selected miRNAs showed consistency with sequencing results except for a few miRNAs like pxy-miR-2a-3p, pxy-miR-2b-3p, and pxy-miR-274-5p.

To better understand the function of differentially expressed miRNAs, putative target genes were predicted using the genome of P. xylostella using RNAhybrid, miRanda, and TargetScan software. Our target prediction results indicated that 30,930 common spots were detected between RNAhybrid and TargetScan, 30,818 between RNAhybrid and miRanda, and 31,942 between TargetScan and miRanda. When the target prediction results of all three software were combined, 30,699 common spots were detected and were selected for further analysis (Figure 4).

The GO enrichment analysis was performed to classify the functions of miRNA target genes. The putative target genes were classified into three major categories, biological process, cellular component, and molecular function, of GO classification. Our results exhibited that cellular process, cell part, and catalytic activity were the most enriched categories in the biological process, cellular component, and molecular function, respectively, at all-time points of treatment (Supplementary Figure 3). Similar to our findings, previously, target genes of Ostrinia. furnacallis in response to Bacillus thuringiensis and Wolbachia- responsive miRNAs in Tetranychus urticae were also categorized into the cellular process, cell part, and catalytic activity (Rong et al., 2014; Xu et al., 2015).

To find out particular signaling pathways of the putative miRNA target genes, Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis was performed. The most enriched categories were transport and catabolism, signal transduction, and cancers in cellular processes, environmental information processing, and human diseases, respectively, at all-time points of infection (Supplementary Figure 4).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.