A combined microRNA and transcriptome analyses illuminates the resistance response of rice against brown planthopper

Background The brown planthopper (BPH, Nilaparvata lugens Stål) is a kind of phloem-feeding pest that adversely affects rice yield. Recently, the BPH-resistance gene, BPH6, was cloned and applied in rice breeding to effectively control BPH. However, the molecular mechanisms underlying BPH6 are poorly understood. Results Here, an integrated miRNA and mRNA expression profiling analysis was performed on BPH6-transgenic (BPH6G) and Nipponbare (wild type, WT) plants after BPH infestation, and a total of 217 differentially expressed miRNAs (DEMs) and 7874 differentially expressed mRNAs (DEGs) were identified. 29 miRNAs, including members of miR160, miR166 and miR169 family were opposite expressed during early or late feeding stages between the two varieties, whilst 9 miRNAs were specifically expressed in BPH6G plants, suggesting involvement of these miRNAs in BPH6-mediated resistance to BPH. In the transcriptome analysis, 949 DEGs were opposite expressed during early or late feeding stages of the two genotypes, which were enriched in metabolic processes, cellular development, cell wall organization, cellular component movement and hormone transport, and certain primary and secondary metabolite synthesis. 24 genes were further selected as candidates for BPH resistance. Integrated analysis of the DEMs and DEGs showed that 34 miRNAs corresponding to 42 target genes were candidate miRNA-mRNA pairs for BPH resistance, 18 pairs were verified by qRT-PCR, and two pairs were confirmed by in vivo analysis. Conclusions For the first time, we reported integrated small RNA and transcriptome sequencing to illustrate resistance mechanisms against BPH in rice. Our results provide a valuable resource to ascertain changes in BPH-induced miRNA and mRNA expression profiles and enable to comprehend plant-insect interactions and find a way for efficient insect control.


Background
Rice is a primary food in China and other Asian countries (Normile 2008). The brown planthopper (BPH) is one of the most harmful insect pests of rice, which in modern rice cultivation causes severe damage and lead to large annual economic losses [1,2]. As a typical vascular-feeding insect, BPH sucks the phloem sap of rice and results in extensive dwarfing, wilting, browning and drying of the plants. Furthermore, BPH serves as a vector to transmit viral disease [1,2]. In the cultivation practice, BPH has developed resistance to most insecticides. The most economic and environment-friendly option for BPH control is to grow resistant rice varieties.
Transcript profiles contribute to our understanding of the defense mechanisms of rice against BPH. Previously, the transcriptional profiles of resistant cultivar B5 and susceptible cultivar MH63 were reported using a cDNA microarray. Expression of genes involved in an array of signaling pathways, oxidative stress, pathogen-related, and macromolecule degradation was evidently enhanced, whilst expression of those involved in the flavonoid pathway, photosynthesis and cell growth was reduced upon BPH infestation [12][13][14]. Recently, a microarray analysis of Rathu Heenati and TN1 under BPH infestation revealed that transcription factors and plant hormones played important roles in the defense response [15,16]. RNA sequencing of the BPH15 introgression line and recipient line before and after infestation by BPH identified chief defense mechanisms associated with transcription factors, hormone signaling pathway, and MAPK cascades [17].
Although BPH responsive transcriptomes profiling of miRNAs and mRNAs have been reported independently, integrated expression profiling of miRNAs and their target genes associated with the interaction of rice and BPH has not been studied. To further reveal the molecular mechanism of rice responding to BPH, highthroughput sequencing was applied to analyze the miRNA and mRNA expression profiles in BPH fed seedlings. Upon integration of these two datasets, a total of 38 miRNAs, 24 genes and 34 miRNAs corresponding to 42 target genes were identified. Our result is a valuable resource for genome-wide studies on BPH responsive genes, and the resistance mechanisms mediated by miR-NAs in rice.

Evaluation of BPH resistance
In this study, a genomic fragment containing BPH6 with its native promoter was transferred into the BPH susceptible wild type (WT), Oryza sativa subsp. japonica cv. Nipponbare, and got BPH6-transgenic plants (BPH6G). The homozygous T 2 transgenic lines were analyzed for BPH resistance using the bulk seedling test. WT plants began to wither on the 4th day and died on the 7th day after BPH infestation, but the BPH6G plants were still alive (Fig. 1a). In the BPH host choice test, the average number of BPHs settled on WT increased rapidly from 6 to 48 h, whereas those on the BPH6G lines remained relatively constant over 72 h (Fig. 1b). Moreover, the ratio of weight gain was significantly less for BPH fed on the BPH6G plants than those on WT from 12 to 72 h (P < 0.01 at 12 h) (Fig. 1c).
In our previous work, the levels of salicylic acid, JA-Ile and cis-zeatin were induced to high levels from 3 to 24 h after BPH infestation in BPH6G compared to WT [11]. Phytohormone synthesis-related genes, PAL (phenylalanine ammonia-lyase), AOS2 (allene oxide synthase 2) and IPT10 (isopentenyl-transferase 10) were selected for expression analysis in BPH6G and WT plants after BPH infestation. Expression of PAL and IPT10 increased more rapidly in the BPH6G plants from 6 to 24 h, whilst the expression levels of AOS2 increased after 48 h in both plants (Fig. 1d-f). RNA was isolated from the leaf sheathes of the BPH6G and WT plants from 0 to 72 h after BPH feeding, and divided into non-infested controls (0 h), early feeding stages (including 6, 12 and 24 h), and late feeding stages (including 48, 60 and 72 h) for high-throughput sequencing analysis.

DEMs in the BPH6G and WT plants before and after BPH feeding
After normalization of the raw sequence reads, the average normalized reads of three independent biological replicates in the libraries were chosen for further analysis. The expression levels of miRNAs were compared amongst the different groups. Using fold changes ≥2, P < 0.05 of the average value of three replicates, 231 DEMs were detected, including 119 DEMs between the different varieties and 217 DEMs between different feeding stages ( Fig. 2a-b). In the early feeding stages, there were more DEMs in WT (89) than in the BPH6G plants (61) (Fig. 2a). In the late feeding stages, the number of up-regulated DEMs (92) were four folds higher than down-regulated ones in WT (Fig. 2a), indicating that serious damage was caused by BPH.
To verify the data in miRNA sequencing, six DEMs were selected for quantitative stem-loop RT-PCR assays [32]. The results were broadly consistent to those from sequencing analysis, although expression of some miR-NAs differed a little (Fig. 2c).
To identify miRNAs related to plant resistance responses, Venn diagrams were used to show the DEMs appeared in the BPH6G plants compared to WT (R0/S0, R_early/S_ early and R_late/S_late) (Fig. 2b). There were 23 overlapping DEMs in the comparisons (Fig. 2b), 18 of which a, Image of the BPH resistance evaluation of the BPH6G and WT seedlings. G stands for BPH6G, N stands for WT; DAI, days after infestation. b, Settling of BPH nymphs on the BPH6G and WT plants in the host choice test. c, BPH weight gain ratio on the BPH6G and WT plants. d-f, Expression of phytohormone synthesis-related genes in the BPH6G and WT plants after BPH infestation. Rice TBP was used as a reference control. Genes expression was quantified relative to values obtained from non-infested WT. Data represent means ± SD of ten independent experiments (b-c) and three independent biological repeats (d-f). Asterisks indicate significant differences revealed by one-way ANOVA at P < 0.05 (*) and P < 0.01 (**), respectively showed opposite expression before and after BPH feeding (Fig. 3a, Additional file 3: Table S2). Members of the miR169 family were up-regulated before BPH feeding (R0/ S0) and down-regulated after BPH feeding (R_early/S_early or R_late/S_late). In contrast, members of miR160 and miR166 families were down-regulated in R0/S0 and upregulated in R_early/S_early or R_late/S_late. The DEMs in early and late feeding stages of the two varieties (S_early/S0, S_late/S0, R_early/R0 and R_late/R0) were analyzed by Venn diagrams (Fig. 2B). A total of 63 DEMs were expressed in R_early/R0 or R_late/R0 and 9 specifically expressed in both R_ early/R0 and R_late/R0 (Fig. 2b). Furthermore, 29 DEMs were opposite expressed in BPH6G and WT plants after BPH feeding (Fig. 3b, Additional file 4: Table S3). Among them, members of the miR169 family, miR156b-3p and miR396c-5p were downregulated, whilst members of the miR160 and miR166 families were up-regulated in BPH6-trangenic plants.
In addition, members of miR1861 and miR319, and other miRNAs appeared opposite expression in BPH6G and WT plants after BPH feeding, or were specifically expressed in both R_early/R0 and R_late/ R0 (Fig. 3b).

General mRNA expression profiles
mRNA libraries were constructed to analyze gene expression and to profile all miRNA targets that were differentially expressed in response to BPH feeding. Total reads of 95,471,364 to 111,697,630 were sequenced from 18 mRNA libraries. After deletion of low-quality reads in samples from the BPH6G plants, 84.70-89.99% of the reads were mapped to 28,988-30,383 rice genes (Additional file 5: Table S4). In the replicates from WT, 82.67-90.84% of the reads were mapped to 28,838-30,006 rice genes (Additional file 5: Table S4).
Considering that some reference genes are suppressed in host-herbivore interaction [33], we carefully selected reference genes with stable expression during BPH infestation for qRT-PCR analysis. Eight frequently used reference genes, eEF1α (Os03g08020), GAPDH (Os02g38920), SDHA (Os07g0424), TBP (Os03g45410), HSP (Os03g31300), β-tubulin (Os03g56810), Ubiquitin (Os03g03920) and LSD1 (Os12g41700) were selected to evaluate the respective FPKM values extracted from our RNA-seq data (Fig. 4a). eEF1α, GAPDH and β-tubulin were significantly reduced in S_late and R_late, and LSD1 was stable but relatively low. Combined with our b, Venn diagrams of the unique and shared DEMs. c, Stem-loop RT-PCR to verify miRNA expression in the BPH6G and WT plants. miRNA expression was normalized by U6. Data represent the mean ± SD of three independent biological repeats. Asterisks indicate significant differences revealed by one-way ANOVA at P < 0.05 (*) and P < 0.01 (**), respectively previous results [11,33], TBP was used as the reference gene for qRT-PCR analysis.

DEGs in the BPH6G and WT plants before and after BPH feeding
There were 8577 DEGs (log 2 FC ≥ 1, FDR < 0.05) detected in this, including 4608 between the different varieties and 7874 between different feeding stages (Table 1). DEGs in the BPH6G and WT plants at different feeding stages were hierarchically clustered. Amongst the four comparisons, the expression patterns of the DEGs were similar, showing consistent up-or down-regulation (Additional file 6: Fig. S2).
During early feeding stages, more DEGs were upregulated in BPH6G plants (1851) compared to WT (965) ( Table 1), and the numbers with FCs ≥ 5 were more in BPH6G plants (590) than in the WT (184). Upregulated DEGs (1851) were three-fold more than downregulated ones (657), and the number of up-regulated DEGs with FCs ≥ 5 (590) were six-fold more than downregulated ones (94) in the BPH6G plants. During late feeding stages, the number of up-regulated DEGs (1356) were similar to that of down-regulated ones (1569) in BPH6G plants. However, during early feeding stages, the down-regulated DEGs (1952) were almost two-fold more than up-regulated ones (965) in WT, indicating the response to BPH-induced wounding and physiological stresses. During late feeding stages, the number of DEGs in WT dramatically increased from 2917 to 6394, and the number with FCs ≥ 5 increased remarkably from 549 to 2854, indicating more serious damage to rice plants caused by BPH.
To verify the RNA-seq data, 30 DEGs were selected for qRT-PCR analysis. The qRT-PCR results were consistent with RNA-seq data, since the genes displayed similar fold-changes with a correlation ratio of R 2 = 0.968 (Fig. 4b).

Identification of genes related to BPH resistance
To investigate the function of BPH6, the sequencing data of BPH6G and WT plants before BPH feeding were compared. There were 3327 DEGs with FC ≥ 2, including 649 up-regulated and 2678 down-regulated ones (Table  1). These DEGs were analyzed by GO (gene ontology) enrichment to explores their functions. The upregulated genes were enriched in defense, protein modification and protein targeting to membrane. Downregulated genes were enriched in the regulation of transcription, signal transduction, cell wall organization and cell proliferation. In addition, these DEGs were enriched in plasma membrane, extracellular region, and cell wall for cellular component ( Fig. 5a-b).
Next, Venn diagrams were used to analyze the possible BPH resistance-related genes in the DEGs of the two rice genotypes. In the BPH6G plants, there were 548 and 1572 DEGs down-and up-regulated respectively after BPH feeding; while in WT, 3127 and 1521 DEGs were respectively down-and up-regulated after BPH feeding (Additional file 7: Fig. S3A). To fully understand the function of these DEGs, GO enrichment analysis were performed. When the biological processes were considered, the up-regulated genes in the BPH6G plants and the down-regulated genes in WT were both enriched in cell wall organization or biogenesis, regulation of biological process, developmental growth, anatomical structure morphogenesis and single-multicellular organism process (Additional file 8: Fig. S4A, D). Down-regulated genes in the BPH6G plants and up-regulated genes in WT were both enriched in single-organism metabolic process, primary metabolic process and response to biotic stimulus and chemical (Additional file 8: Fig. S4B, C). Genes associated with hydrolase activity, Ras guanylnucleotide exchange factor activity and protein binding were most contrasting amongst the molecular function GO terms in the two rice varieties (Additional file 8: Fig.  S4). Three cellular components of GO terms, external encapsulating structure, vesicle and intrinsic component of membrane were enriched, suggesting involvement of cell wall, vesicle and plant membrane in the response to BPH feeding (Additional file 8: Fig. S4).
To further streamline potential BPH resistance-related genes, the opposite expression DEGs during early and late feeding stages of two varieties were assessed. There were 949 DEGs in the BPH6G and WT plants after BPH feeding, of which, 935 were up-regulated in the BPH6G   Fig. S3B). GO enrichment analysis indicated that these resistance-related genes were enriched in metabolic process, cellular development, cell wall organization, cellular component movement and hormone transport for biological process, and membranebounded vesicle and cell wall for cellular component (Fig. 5c). For further information regarding the molecular and biochemical responses of rice after BPH infestation, BPH responsive DEGs were combined with KEGG processes (Kyoto Encyclopedia of Genes and Genomes). At the P < 0.05, the BPH responsive DEGs were enriched in key pathways. The up-regulated DEGs were involved in primary and secondary metabolite processes, such as limonene and pinene degradation, starch and sucrose metabolism, stilbenoid, diarylheptanoid and gingerol biosynthesis, and brassinosteroid biosynthesis. In contrast, amino and nucleotide sugar metabolism and diterpenoid biosynthesis were remarkably enriched among the down-regulated genes (Additional file 9: Table S5). Finally, 24 genes were differentially expressed in both the BPH6G and WT plants after BPH feeding, and were considered BPH resistance-related genes ( Table 2). Of these DEGs, 23 were dramatically up-regulated in the BPH6G plants and down-regulated in WT after BPH infestation. A single gene was down-regulated in the BPH6G plants and up-regulated in WT. Among them, two genes encoding germin-like proteins, two lipid transfer proteins, two cytochrome P450 family proteins and two Rop guanine nucleotide exchange factors played important roles against BPH. The majority of these   genes were enriched in response to stress, transport, cell wall organization, Rho GTPase activity and oxidationreduction process, and were enriched in the extracellular region, cell wall, membrane and nucleus (Table 2).

Integrated analysis of miRNA and mRNA expression profiles
In most cases, miRNAs negatively regulate target mRNA through translation repression or mRNA degradation    [19]. To correlate the identified miRNAs with their target genes, the psRNA target tool was used to predict miRNA targets on mRNAs using the parameters fold changes ≥2, P < 0.05 [34]. There were 89, 117, 61 and 92 DEMs that significantly and negatively correlated with 488, 1096, 235 and 498 target mRNAs in S_early/S0, S_ late/S0, R_early/R0 and R_late/R0, respectively (Fig. 2a).
To identify potential miRNA-mRNA pairs related to BPH resistance, 70 DEMs in R_early/R0 or R_late/R0 (Fig. 2b) and 29 DEMs at different feeding stages (Fig. 3b) were selected and negatively correlated with 656 target mRNAs (Additional file 10: Table S6). These miRNAs target different mRNAs during each feeding stage. For example, miR156b-3p was upregulated in S_early/S0 and down-regulated in R_late/ R0, which negatively correlated with 20 downregulated target genes in S_early/S0 and 4 upregulated ones in R_late/R0, respectively. However, the majority of these targets showed a similar trend of expression in the BPH6G and WT plants after BPH feeding (Additional file 10: Table S6). Excluding these miRNAs and their corresponding targets, 34 miRNAs corresponding to 42 target genes were differentially expressed in R_early/R0 or R_late/R0, or opposite expressed in the BPH6G and WT plants after BPH feeding, and selected as BPH resistance-related miRNA-mRNA candidates (Table 3).
miR156b-3p and miR169i-5p.2 with their targets encoding GDSL-like lipase (GDSL/LOC_Os02g40440) and Leucine Rich Repeat family protein (LRR/LOC_ Os11g36180) respectively were selected for validation in rice protoplasts. Two plasmids of each pair, one encoding pri-miRNA, and the other YFP and HA fused target, were transfected into the protoplasts. In both cases, the inflorescence signal of the blank YFP plasmid could not be weakened by the pri-miRNAs, however, that of the targets could be significantly weakened by the respective pri-miRNAs (Fig. 6b). Western blot verified the results of the YFP signal at the protein level (Fig. 6c). These results indicate that miR156b-3p and miR169i-5p.2 downregulate GDSL and LRR expression in rice cells, respectively.

Discussion
Few studies have reported the use of combined miRNA and mRNA expression profiles to analyze the responses of herbivore insects in plants, excluding studies on aphidinduced miRNA expression [35]. This study was the first to report the combined analysis of miRNA and mRNA expression profiles in BPH-infested rice, enhancing our understanding of the regulatory mechanisms of miRNA-mRNA in rice after BPH attack.
In this study, the average number and the weight gain rate of BPHs increased rapidly from 6 to 48 h and remained gently after 48 h on WT (Fig. 1b-c). In addition, there were significant differences in the expression of hormone-related genes before and after 48 h in the BPH6G plants ( Fig. 1d-f). These results demonstrate that the defense establishment and significant progression of BPH6G plants exists up to 48 h after BPH infestation. Therefore, RNA was divided into three groups, non-infested, early feeding stage (before 48 h) and late feeding stage (after 48 h).
To study BPH resistance related genes a wide range, RNA-sequencing analysis was performed in the BPH6G and WT plants under BPH infestation. Transcriptome analysis revealed notable differences in the response of the BPH6G and WT plants to BPH feeding. The inducible defense responses against BPH in BPH6G plants were more robust during early feeding stages compared with WT as a larger number of up-regulated DEGs (FCs ≥ 2) were detected in the BPH6G. In contrast, a larger number of DEGs were detected in WT during early and late feeding stages, indicating remarkable metabolic and physiological changes in WT after BPH feeding due to the absence of BPH resistance. In addition, upregulated DEGs were much higher than down-regulated ones in the BPH6G plants, suggesting that the expression of genes associated with resistance in the BPH6G plants was up-regulated.
Previously, the transcript profiles of resistant rice cultivars revealed key defense mechanisms related to transcription factors, hormone signaling, MAPK cascades and pathogen-related genes. In this study, the DEGs were analyzed during early and late feeding stages in the two varieties to reveal the BPH6-mediated defense mechanisms. The GO enrichment analysis of the DEGs of BPH6G and WT plants before BPH feeding indicate that BPH6 takes part in defense and stress, and other developmental and physiological process. There were 949 DEGs opposite expressed at early or late feeding stages between the two varieties, most of which were upregulated in the BPH6G plants, suggesting that the majority positively regulate rice immunity against BPH. These DEGs were enriched in cellular development, cell wall organization, cellular component movement and hormone transport (Fig. 5c), which were consistent with the function of BPH6, which promotes exocytosis, participates in cell wall maintenance and reinforcement, and activates hormone signaling after BPH feeding [11]. In addition, the up-regulated DEGs were involved in the primary and secondary metabolite processes, suggesting that these metabolites play important roles in rice defense against BPH. Finally, 24 genes were selected as potential candidates for BPH resistance ( Table 2). Most of the genes, excluding LOC_Os02g45420 and LOC_ Os01g72270, were highly upregulated (FC > 5) during the early stages in the BPH6G plants and highly downregulated during the late stages in WT, indicating their important roles in BPH response. The germin-like protein (GLP) gene family confers broad-spectrum resistance to pathogens and insects in plants through H 2 O 2 production due to superoxide dismutase activity at the infection site [36,37]. The overexpression of LTPs increases the resistance to pathogens and environmental stresses due to the hydrophobic protective layers of surface polymers [38]. Pectinesterase plays a regulatory role in mechanical stability and elongation of the cell wall in response to pathogen invasion in Arabidopsis [39]. Fasciclin-like arabinogalactan-proteins are implicated in plant growth and development, cell wall remodeling, hormone signaling modulation and pathogen defenses [40]. In addition, NB-ARC proteins, MYB transcription factors, ethylene response factors and HSP20 are all involved in pathogen resistance [41][42][43][44].
Integrated miRNA and mRNA expression analysis can help identify the functional miRNA-mRNA pairs related to host-insect interaction. In this study, 70 specific DEMs in the BPH6G plant (Fig. 2b) and 29 oppositely expressed miRNAs (Fig. 3b) corresponding to 656 target genes were detected under BPH attacking (Additional file 10: Table S6). However, only 34 miRNAs corresponding to 42 target genes might be potentially related to BPH response (Table 3). For example, the members of miR166 family, reported to positively regulate rice immunity against the blast fungus [26], were up-regulated in the BPH6G plants after BPH feeding (Fig. 3b). However, the targets of miR166 exhibited similar trend of expression in the BPH6G and WT plants after BPH feeding (Additional file 10: Table S6). Therefore, miR166 and its targets were excluded as BPH-related candidates. This phenomenon can be explained by the following: (1) most targets had the same expression trends in the BPH6G and WT plants after BPH feeding, (2) plants defense responses to insects include both systematic and local responses, and many targets may not be expressed at this point and (3) the accepted criteria for the DEMs and DEGs may miss key interactions. After integrated analysis of the DEMs and DEGs, several important miRNA-mRNA pairs involved in BPH stress were identified. miR156b-3p targeted to GDSL-like lipase in response to BPH (Table 3, Fig. 6). Previous studies have shown that miR156 silencing confers enhanced resistance to BPH [30], and GDSL lipases modulate immunity through lipid homeostasis [45]. Therefore, miR156b-3p may negatively regulate BPH resistance by targeting GDSL lipases. Members of the miR169 family, including miR169g/h/i-5p.1/i-5p.2/j/k/l/m/o, target some leucin rich repeat family proteins that play key roles in pattern recognition and the initiation of downstream responses [46]. In addition, members of miR1861 family target auxin response factors, miR396c-3p targets abscisic acid 8′-hydroxylase gene, and miR5513 targets ethylene-responsive transcription factor, suggesting that auxin, ABA and ethylene might all involved in the BPH response.

Conclusion
In this study, 18 libraries were constructed for the BPH6G and WT genotypes before and after BPH feeding. These libraries were amplified and sequenced, and miRNAs and mRNAs related to BPH resistance were identified. We identified members of miR160, miR166, miR169, miR1861, miR319 and miR390 families, and other miRNAs that played important roles in the BPH6-mediated resistance to BPH. DEGs potentially involved in BPH responses included genes related to metabolic process, cellular development, cell wall organization, cellular component movement and hormone transport. Additionally, 34 miRNAs corresponding to 42 target genes were identified as candidates for BPH resistance miRNA-mRNA pairs. The integrated analyses of miRNAs and genes related to BPH resistance in rice provide the basis for further research probing the functions of miRNA and targets in the BPH response, and establish a molecular basis for further studies on how plants respond to BPH infestation.

Plant and insect materials
A 7.8 kb DNA fragment containing the BPH6 gene with its native promoter, was amplified from Swarnalata (IRRI Acc. No. 33964), digested with KpnI and inserted into the binary vector pCAMBIA1300, transformed into the susceptible wild type (WT) Nipponbare (IRRI Acc. No. 136196) through Agrobacterium-mediated method, and identified by Zhang et al. [47]. A voucher specimen of the BPH6-transgenic line has been deposited in the China Center for Type Culture Collection (No. P201907). Seeds were grown in plastic cups (9 cm in diameter and 15 cm in height) with 15 plants per cup, and maintained in a greenhouse with cycles of 32 ± 2°C/ 14 h light and 26 ± 2°C/10 h dark periods.
The BPH population in this study was were kept in the laboratory and maintained on 1-month-old plants of the susceptible rice cv Taichung Native1 (IRRI Acc. No. 00105) under controlled environmental conditions as described above in Wuhan University [48].

Evaluation of rice resistance to BPH
At the four-leaf stage, the BPH6G and WT plants were infested with 8 s-instar BPH nymphs per seedling, and checked each day until all seedlings of WT died. Evaluations were carried out with three biological repeats for each line.
Host choice test was carried out as described [11], WT and BPH6G plants were grown diagonally in each bucket. At the four-leaf stage, twenty BPH nymphs at the third-instar stage were release in the buckets, and the number of nymphs settled per plant were counted at 6, 12, 24, 36, 48, 60 and 72 h after release. Ten buckets for each line were analyzed.
In BPH weight gain analysis, newly emerged brachypterous females, Parafilm sachets, 1-month-old WT and BPH6G lines were used as described by Shangguan et al. [48]. Weight increase relative to the initial weight were calculated BPH weight gain ratios. Experiments were performed five times with 10 replicates for each line.

Sample collection
The endpoint method was used to collect samples through BPH treatment [29]. Although all processes began at different time, they were stopped at the same specified time. Seedlings were infested with 8 s-instar BPH per seedling at the four-leaf stage after 0, 6, 12, 24, 48, 60, and 72 h. For analysis, three biological replicates per treatment with 15 seedlings per replicate were used. Leaf sheaths were mixed for non-infested controls (0 h), infested early (6, 12 and 24 h) and infested late (48, 60 and 72 h). Samples were referred to as R0, R_early, and R_late for the BPH6G lines, and S0, S_early, and S_late for WT. Leaf sheathes were cut and frozen in liquid nitrogen, and stored at − 80°C until use.

RNA extraction
Infested and mock leaf sheathes were used for total RNA extraction using commercial RNAiso Plus kits (TaKaRa, code no. 9109). Concentrations of RNA were checked using Qubit fluorometric quantitation (Thermo Fisher Scientific, Wilmington, DE, USA), and integrity was verified on a Bio-Analyzer 2200 (Agilent Technologies, Waldbronn, Germany).
Construction of the cDNA library and RNA mapping cDNA library for each sample was constructed using NEBNext® Ultra™ directional RNA library prep kits (NEB, code no. E7420S), and quantified on a 150 bp paired-end run by Agilent2200 and sequenced by HiSeq X (Illumina, San Diego, CA). Clean reads were obtained after the removal of adaptors, low quality reads and reads with > 5% unknown nucleotides, and mapped on rice genome (TIGR7) using the Hisat2 [49]. Gene counts were obtained by HTseq and gene expression was determined using the RPKM method [50].
miRNA library construction, sequencing and mapping miRNA libraries were prepared using Ion Total RNA-Seq Kit v2.0 (Thermo Fisher, code no. 4475936). miRNA for construction were selected according to size by polyacrylamide gel electrophoresis and processed for Proton Sequencing as per commercially available protocols. A total of 18 small RNA libraries were constructed with the BPH6G plants (R) and WT (S) infested by BPH for non-infested (R0, S0), early feeding stages (R_early, S_ early) and late feeding stages (R_late, S_late).
After deep sequencing, the raw data were evaluated in FAST-QC software (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/), including the quality distribution of nucleotides, position specific sequencing quality, GC content, proportion of PCR duplication and k-mer frequency. Raw data were processed to remove low-quality reads, adaptor sequences, contaminant reads, and reads of < 20 nt and > 24 nt. All of the sequences were aligned in the NCBI GenBank (release 227.0) and Rfam (release 13.0) database, and mapped to the rice genome to identify and remove rRNA, tRNA, scRNA, snoRNA, snRNA and small RNAs mapped to exons or introns and repeat sequences (Additional file 2: Fig.  S1B).

Differential expression analysis of miRNAs and genes
Differentially expressed miRNAs and genes were filtered by EB-Seq algorithm after significance. P-values and FDR analyses were performed at absolute values of log 2 FC ≥ 1, P < 0.05, FDR < 0.05 [51].

Target analysis
The psRNA target software (http://plantgrn.noble.org/ psRNATarget/) was used to predict miRNA targets on mRNAs based on the default parameters.

Analysis of GO (gene ontology) and KEGG pathway
GO annotations from NCBI (http://www.ncbi.nlm.nih. gov/) and GO (http://www.geneontology.org/) were downloaded. To identify DEGs pathways, the KEGG database was used. To identify significant GO and pathway categories, Fisher's exact tests were applied under absolute values of P < 0.05 and FDR < 0.05 [52].

qPCR analysis of miRNAs and mRNAs
For first strand cDNA synthesis, 2 μg total RNA were extracted using PrimeScript™ RT reagent Kits accompanied with gDNA Eraser (TaKaRa, code no. RR047A) and miRNAs were extracted using miRcute Plus miRNA First-Strand cDNA Kits (TIANGEN, code no. KR211). miRNAs were quantified by stem-loop RT-PCR [32]. Gene expression was analyzed by qPCR using SYBR green supermixes from Bio-Rad and CFX96 real-time system. Each experiment was performed in three biological replicates. The expression of miRNAs and genes were calculated through the 2 -ΔΔC (t) method [53] with internal reference genes TBP and U6, respectively. Primers are listed in Additional file 11: Table S7. Oneway ANOVA was used for statistical analyses in Microsoft Excel.

Validation of the predicted target genes of miRNAs
The role of miRNAs on the targets were investigated through counting the fluorescent cells [29]. One plasmid encoded pri-miRNA (miR156b-3p, miR169i-5p.2) was amplified from WT DNA and cloned into the binary vector pCXUN. The other containing the targets (GDSL, LRR) were amplified from WT cDNA and cloned into the binary vector pCXUN with YFP genes and HA tags. Constructs expressing miRNAs and the targets were transiently co-transfected into rice protoplasts isolated from 10-day-old WT stems. The fluorescent cells were imaged and numbered using a confocal microscope (FV10-ASW, Olympus). Protein expression was determined by Western blotting. Primers used in the experiments are listed in Additional file 11: Table S7.