Skip to main content
Download PDF
  • Research article
  • Open access
  • Published:

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

BMC Genomics volume 21, Article number: 144 (2020) Cite this article

Abstract

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.

Since the report of the first BPH-resistance gene, BPH1 in 1969 [3], more than 30 ones have been identified and mapped from wild and cultivated rice germplasms [2]. 12 of them, BPH14, BPH26/2, BPH3, BPH29, BPH32, BPH18 and BPH9/1/7/10/21 were characterized by map-based cloning approaches [4,5,6,7,8,9,10]. The structure and localization of BPH-resistance proteins provides a model system for determining the molecular basis of rice-BPH interaction. BPH14 encodes an NLR (nucleotide-binding and leucine-rich repeat) protein that localizes to the cytoplasm and nucleus [4], BPH1/2/7/9/10/18/21/26 encode NLR proteins that localize to the endomembrane [10], BPH3 encodes a lectin receptor kinase that localize to the plasma membrane [6] whilst BPH29 encodes a nucleus-localized B3 domain-containing protein [7]. Recently, we cloned another BPH-resistance gene, BPH6, encoding a yet uncharacterized protein in the exocyst and interacts with OsEXO70E1, an exocyst subunit [11]. However, the BPH6-mediated molecular mechanisms against BPH remain largely undefined.

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].

MicroRNAs (miRNAs) are ~ 21-nucleotide-long regulatory RNAs produced from endonucleolytic processing of single-stranded hairpin precursors in animal and plant [18]. miRNAs specifically regulate target gene expression through binding complementary sequences to degrade mRNA or inhibit translation [19]. Plant miRNAs are involved in many development processes, including hormone signal transduction, and leaf, floral, shoot, root and vascular development [20,21,22], and play significant roles in abiotic and biotic stress responses [23,24,25,26,27,28]. miR160 is associated with local defense and systemic acquired resistance to potato late blight [24]. miR166, miR169 and miR319 participate in the regulation of rice immunity against the blast fungus Magnaporthe oryzae [25,26,27,28]. However, few miRNAs have been revealed functioning in insect response. BPH-responsive miRNAs were investigated from resistant rice in comparison with susceptible plants [29]. miR156 and miR396 negatively regulated BPH resistance through regulating Jasmonic acid (JA) and flavonoid biosynthetic pathways, respectively [30, 31].

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, high-throughput 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 miRNAs in rice.

Results

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 T2 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).

Fig. 1
figure 1

BPH resistance evaluation of the BPH6G and WT plants. 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

Full size image

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.

Small RNA library construction and sequencing

Total reads of 9,034,925 to 14,016,694 were got in S0, S_early, S_late, R0, R_early and R_late libraries, respectively (Additional file 1: Table S1). After removing all low-quality reads, poly A, incorrect adaptors and reads of < 20 nt and > 24 nt, a total of 4,503,508 to 8,547,717 clean reads remained in the 18 libraries, respectively (Additional file 1: Table S1). In the 18 libraries, the main size classes were 21 nt and 24 nt followed by 22 nt and 23 nt as previously reported for rice small RNAs (Additional file 2: Fig. S1A). Approximately 20.16–23.09% (S0), 23.68–26.35% (S_early), 19.88–24.39% (S_late), 30.66–31.04% (R0), 31.37–35.58% (R_early) and 23.57–23.67% (R_late) of the clean reads were assigned to the miRBase database (Additional file 1: Table S1). Rice miRNA is the most thoroughly studied monocot miRNA, and 738 mature miRNAs were identified in the miRbase (release 22). Accordingly, we analyzed the 738 known miRNAs in our data.

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.

Fig. 2
figure 2

Differentially expressed miRNAs in the comparisons. a, Number of miRNAs and target genes up- or down-regulated in all comparisons (fold change > 2, P < 0.05). 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

Full size image

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 miRNAs differed a little (Fig. 2c).

Identification of miRNAs related to BPH resistance

miRNA expression differences in BPH resistant and susceptible rice before BPH attack was first compared. There were 55 DEMs, including 24 up-regulated and 31 down-regulated miRNAs in R0/S0 (Fig. 2a), many of which belonged to known miRNA families including miR156, miRNA160, miR166, miR169, miR1846, miR1861 and miR319 (Additional file 3: Table S2). Members of the miRNA families were reported to be involved not only in growth, development, grain size and hormone signaling, but also in response to biotic and abiotic stress [21,22,23,24,25,26,27,28,29,30,31]. These BPH6 responsive DEMs might be involved in response to BPH.

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 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 up-regulated in R_early/S_early or R_late/S_late.

Fig. 3
figure 3

Potential BPH resistance-related miRNAs. a, Overlapping DEMs appeared in both varieties before and after BPH feeding. b, Overlapping DEMs appeared at the early and late feeding stages of both varieties and specially DEMs in the BPH6G plants. Colors represent fold-change values

Full size image

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 down-regulated, 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 previous results [11, 33], TBP was used as the reference gene for qRT-PCR analysis.

Fig. 4
figure 4

Expression profiles of mRNAs. a, FKPM values of eEF1α, GAPDH, SDHA, TBP, HSP, β-tubulin, Ubiquitin and LSD1 from RNA-sequencing data. b, Pearson correlation scatter plots of comparisons of gene expression fold-changes measured by sequencing and qRT-PCR. Rice TBP was used as an internal reference. Gene expression was quantified relative to values obtained from non-infested samples. Data represent means of three independent biological repeats. Both x and y-axes are shown in the log2 scale. Pearson’s correlation coefficient is indicated by R

Full size image

DEGs in the BPH6G and WT plants before and after BPH feeding

There were 8577 DEGs (log2FC ≥ 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).

Table 1 Differential expression genes between the BPH6-transgenic and Nipponbare plants in response to BPH feeding
Full size table

During early feeding stages, more DEGs were up-regulated 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). Up-regulated DEGs (1851) were three-fold more than down-regulated ones (657), and the number of up-regulated DEGs with FCs ≥ 5 (590) were six-fold more than down-regulated 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 R2 = 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 up-regulated genes were enriched in defense, protein modification and protein targeting to membrane. Down-regulated 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).

Fig. 5
figure 5

GO (Gene Ontology) analysis. Biological process, cellular component, and molecular function of up- (a) and down-regulated DEGs (b) in R0/S0, and the opposite expression DEGs (c) at the early or late feeding stages of the two varieties (P < 0.05). The x- and y-axis indicate the names of the clusters and the number of genes in a category, respectively

Full size image

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 guanyl-nucleotide 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 plants and down-regulated in WT, whilst 14 were down-regulated in the BPH6G plants and up-regulated in WT (Additional file 7: 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 membrane-bounded 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 oxidation-reduction process, and were enriched in the extracellular region, cell wall, membrane and nucleus (Table 2).

Table 2 Candidate BPH resistance-related genes exhibiting opposite expression in early and late feeding stages of two rice genotypes
Full size table

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). In addition, 55, 24 and 70 DEMs negatively correlated with 269, 88 and 369 target mRNAs in R0/S0, R_early/S_early and R_late/S_late, 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 up-regulated in S_early/S0 and down-regulated in R_late/R0, which negatively correlated with 20 down-regulated target genes in S_early/S0 and 4 up-regulated 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).

Table 3 The miRNA-mRNA interactions related to plant resistance
Full size table

To validate whether these miRNAs negatively regulate target expression, four miRNAs and their targets were selected for qRT-PCR verification. The results indicated that miR156b-3p negatively regulated LOC_Os02g40440 and LOC_Os03g06940 in S_early/S0, and LOC_Os07g05370 in R_late/R0, miR396c-5p negatively regulated LOC_Os02g47470, LOC_Os03g21800 and LOC_Os03g47140 in R_late/R0, and miR169g/o and miR2871a-3p negatively regulated LOC_Os07g31840. Furthermore, three target genes (LOC_Os05g38980, LOC_Os07g46560 and LOC_Os11g36180) were down-regulated by miR169i-5p.2 in S_early/S0 and S_late/S0, up-regulated in R_early/R0, while unaffected in R_late/R0 (Fig. 6a).

Fig. 6
figure 6

Negative regulation of miRNAs on their target genes. a, Contrasting expression patterns of miRNAs and their targets. Data represent the mean ± SD of three independent biological experiments. b, Fluorescence micrographs of rice protoplasts transfected with blank YFP plasmids, miRNAs and target gene plasmids. Scale bar, 50 μm. c, Western blot analysis of YFP and the target genes in rice protoplasts using anti-HA and anti-GAPDH antibodies

Full size image

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 down-regulate 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 aphid-induced 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).

Through the comparison of miRNA expression of the BPH6G and WT plants before and after BPH attack, a total of 217 known DEMs were identified (Fig. 2a-b). To identify miRNAs related to BPH response, the DEMs amongst seven comparisons: R0/S0, R_early/S_early, R_late/S_late, S_early/S0, S_late/S0, R_early/R0 and R_late/R0 were analyzed using Venn diagrams (Fig. 2b). 18 DEMs appeared in two of the comparisons (R_early/R0 and R_late/R0), whilst 61 DEMs appeared in both S_early/S0 and S_late/S0 (Fig. 2b), suggesting a lower number of miRNAs were involved in BPH defense responses in the BPH6G plants. 18 miRNAs were opposite expressed before and after BPH feeding in comparison to R0/S0, R_early/S_early and R_late/S_late, and 29 miRNAs were opposite expressed in the BPH6G and WT plants after BPH feeding (Fig. 3). Seventeen of the miRNAs, miR160c-3p/e-3p, miR166a-5p/e-5p/h-5p, miR169h/i-5p.1/i-5p.2/j/k/l, miR1859, miR1861b/j/l, miR390-3p and miR396c-3p, were present in both groups, suggesting their involvement in the defense response of rice against BPH are consistent with their respective roles in pathogen defense [24,25,26]. Previously, the BPH-responsive miRNAs were identified in a BPH15 introgression line [29]. In both the BPH15 introgression and the BPH6G lines, some miRNA exhibited similar expression trends, such as miR156b-3p, miR169h/i-3p/i-5p.1/i-5p.2/j/k/l/m/o, miR396c-5p, miR399j, miR530-5p and miR5513, suggesting a conserved and diverse resistance mechanisms against BPH mediated by BPH6 and BPH15.

miR156, miR160, miR166, miR169 and miR396 were reported to participate in rice immune response against pathogens and insects. miR160 positively regulated potato defense to late blight [24]. miR166 positively regulated rice immunity against the blast fungus via post-transcriptional control of EIN2 [26]. miR169, miR156 and miR396 negatively regulated rice immunity against the pathogens and BPH, respectively [25, 30, 31]. In the BPH6G plants, members of miR160 and miR166 family were upregulated, whilst those of miR156, miR396 and miR169 families were downregulated (Fig. 3b), implying their involvement in BPH response. In addition, miR1859 showed higher expression during heat stress treatment [22]. Members of the miR1861 family regulate starch accumulation and yield in rice [21], whilst miR390-TAS3-ARFs forms an auxin-responsive regulatory network controlling root growth [20]. These reports suggest that the miRNAs indirectly participate in BPH stress responses through altering the metabolic processes and hormone regulation. Interestingly, miR319 negatively regulated immunity to rice ragged stunt virus and blast fungus by repressing the expression of OsTCP21, leading to decrease JA-mediated defenses [27, 28], whilst members of the miR319 family were upregulated in the BPH6G plants (Fig. 3b), suggesting that miR319 might separately regulate rice immunity against BPH and pathogens.

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, up-regulated 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 up-regulated 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 H2O2 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.

Methods

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 log2FC ≥ 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. One-way 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.

Availability of data and materials

The raw sequence data of small RNA and transcriptome during this study could be found in the National Center for Biotechnology Information (NCBI) under the accession number GSE123148. This public accession is currently private and is scheduled to be released on Dec 31, 2019. I have received administrative permission to access and use these.

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE123148.

Abbreviations

AOS2:

Allene oxide synthase 2

BPH:

Brown planthopper

DEG:

Differentially expressed gene

DEM:

Differentially expressed microRNA

FPKM:

Fragments per kilobase of exon per million fragments

IPT10:

Isopentenyl-transferase 10

IRRI:

International Rice Research Institute

MiRNA:

MicroRNA

PAL:

Phenylalanine ammonia-lyase

ScRNA:

Small cytoplasmic RNA

SnoRNA:

Small nucleolar RNA

SnRNA:

Small nuclear RNA

References

  1. Cheng X, Zhu L, He G. Towards understanding of molecular interactions between rice and the brown planthopper. Mol Plant. 2013;6:621–34.

    Article  CAS  PubMed  Google Scholar 

  2. Jing S, Zhao Y, Du B, Chen R, Zhu L, He G. Genomics of interaction between the brown planthopper and rice. Curr Opin Insect Sci. 2017;19:82–7.

    Article  PubMed  Google Scholar 

  3. Pathak MD, Cheng CH, Fortuno ME. Resistance to Nephotettix impicticeps and Nilaparvata lugens in varieties of rice. Nature. 1969;223:502–4.

    Article  Google Scholar 

  4. Du B, Zhang W, Liu B, Hu J, Wei Z, Shi Z, He R, Zhu L, Chen R, Han B, et al. Identification and characterization of Bph14, a gene conferring resistance to brown plant hopper in rice. Proc Natl Acad Sci U S A. 2009;106:22163–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Tamura Y, Hattori M, Yoshioka H, Yoshioka M, Takahashi A, Wu J, Sentoku N, Yasui H. Map-based cloning and characterization of a brown planthopper resistance gene BPH26 from Oryza sativa L. ssp indica cultivar ADR52. Sci Rep. 2014;4:5872.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  6. Liu Y, Wu H, Chen H, Liu Y, He J, Kang H, Sun Z, Pan G, Wang Q, Hu J, et al. A gene cluster encoding lectin receptor kinases confers broad-spectrum and durable insect resistance in rice. Nat Biotechnol. 2015;33:301–5.

    Article  CAS  PubMed  Google Scholar 

  7. Wang Y, Cao L, Zhang Y, Cao C, Liu F, Huang F, Qiu Y, Lou X. Map-based cloning and characterization of BPH29, a B3 domain-containing recessive gene conferring brown planthopper resistance in rice. J Exp Bot. 2015;66:6035–45.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Ren J, Gao F, Wu X, Lu X, Zeng L, Lv J, Su X, Lou H, Ren G. Bph32, a novel gene encoding an unknown SCR domaincontaining protein, confers resistance against the brown planthopper in rice. Sci Rep. 2016;6:37645.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Ji H, Kim SR, Kim YH, Suh JP, Park HM, Sreenivasulu N, Misra G, Kim SM, Hechanova SL, Kim H, et al. Map-based cloning and characterization of the BPH18 gene from wild rice conferring resistance to brown planthopper (BPH) insect pest. Sci Rep. 2016;6:34376.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Zhao Y, Huang J, Wang Z, Jing S, Wang Y, Ouyang Y, Cai B, Xin XF, Liu X, Zhang C, et al. Allelic diversity in an NLR gene BPH9 enables rice to combat planthopper variation. Proc Natl Acad Sci U S A. 2016;113:12850–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Guo J, Xu C, Wu D, Zhao Y, Qiu Y, Wang X, Ouyang Y, Cai B, Liu X, Jing S, et al. Bph6 encodes an exocyst-localized protein and confers broad resistance to planthoppers in rice. Nat Genet. 2018;50:297–306.

    Article  CAS  PubMed  Google Scholar 

  12. Zhang F, Zhu L, He G. Differential gene expression in response to brown planthopper feeding in rice. J Plant Physiol. 2004;161:53–62.

    Article  CAS  PubMed  Google Scholar 

  13. Yuan H, Chen X, Zhu L, He G. Identification of genes responsive to brown planthopper Nilaparvata lugens Stål (Homoptera: Delphacidae) feeding in rice. Planta. 2005;221:105–12.

    Article  CAS  PubMed  Google Scholar 

  14. Wang Y, Wang X, Yuan H, Chen R, Zhu L, He R, He G. Responses of two contrasting genotypes of rice to brown planthopper. Mol Plant-Microbe Interact. 2008;21:122–32.

    Article  CAS  PubMed  Google Scholar 

  15. Wang Y, Guo H, Li H, Zhang H, Miao X. Identification of transcription factors potential related to brown planthopper resistance in rice via microarray expression profiling. BMC Genomics. 2012;13:687.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Li C, Luo C, Zhou Z, Wang R, Ling F, Xiao L, Lin Y, Chen H. Gene expression and plant hormone levels in two contrasting rice genotypes responding to brown planthopper infestation. BMC Plant Biol. 2017;17:57.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  17. Lv W, Du B, Shangguan X, Zhao Y, Pan Y, Zhu L, He Y, He G. BAC and RNA sequencing reveal the brown planthopper resistance gene BPH15 in a recombination cold spot that mediates a unique defense mechanism. BMC Genomics. 2014;15:674.

    Article  PubMed  PubMed Central  Google Scholar 

  18. Axtell MJ, Meyers BC. Revisiting criteria for plant MicroRNA annotation in the era of big data. Plant Cell. 2018;30:272–84.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Bartel DP. MicroRNAs: target recognition and regulatory functions. Cell. 2009;136:215–33.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Marin E, Jouannet V, Herz A, Lokerse AS, Weijers D, Vaucheret H, Nussaume L, Crespi MD, Maizel A. miR390, Arabidopsis TAS3 tasiRNAs, and their AUXIN RESPONSE FACTOR targets define an autoregulatory network quantitatively regulating lateral root growth. Plant Cell. 2010;22:1104–17.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Peng T, Sun H, Qiao M, Zhao Y, Du Y, Zhang J, Li J, Tang G, Zhao Q. Differentially expressed microRNA cohorts in seed development may contribute to poor grain filling of inferior spikelets in rice. BMC Plant Biol. 2014;14:196.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  22. Magnrauthia SK, Bhogireddy S, Agarwal S, Prasanth VV, Voleti SR, Neelamaraju S, Subrahmanyam D. Genome-wide changes in microRNA expression during short and prolonged heat stress and recovery in contrasting rice cultivars. J Exp Bot. 2017;68:2399–412.

    Article  CAS  Google Scholar 

  23. Zhang S, Yue Y, Sheng L, Wu Y, Fan G, Li A, Hu X, Shangguan M, Wei C. PASmiR: a literature-curated database for miRNA molecular regulation in plant response to abiotic stress. BMC Plant Biol. 2013;13:33.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Natarajan B, Kalsi HS, Godbole P, Malankar N, Thiagarayaselvam A, Siddappa S, Thulasiram HV, Chakrabarti SK, Banerjee AK. MiRNA160 is associated with local defense and systemic acquired resistance against Phytophthora infestans infection in potato. J Exp Bot. 2018;69:2023–36.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Li Y, Zhao SL, Li JL, Hu XH, Wang H, Cao XL, Xu XJ, Zhao ZX, Xiao ZY, Yang N, et al. Osa-miR169 negatively regulates rice immunity against the blast fungus Magnaporthe oryzae. Front Plant Sci. 2017;8:2.

    PubMed  PubMed Central  Google Scholar 

  26. Salvador-Guirao R, Hsing YI, San SB. The polycistronic miR166k-166h positively regulates rice immunity via post-transcriptional control of EIN2. Front Plant Sci. 2018;9:337.

    Article  PubMed  PubMed Central  Google Scholar 

  27. Zhang C, Ding Z, Wu K, Yang L, Li Y, Yang Z, Shi S, Liu X, Zhao S, Yang Z, et al. Suppression of jasmonic acid-mediated defense by viral-inducible microRNA319 facilitates virus infection in rice. Mol Plant. 2016;9:1302–14.

    Article  CAS  PubMed  Google Scholar 

  28. Zhang X, Bao Y, Shan D, Wang Z, Song X, Wang Z, Wang J, He L, Wu L, Zhang Z, et al. Magnaporthe oryzae induces the expression of a microRNA to suppress the immune response in rice. Plant Physiol. 2018;177:352–68.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Wu Y, Lv W, Hu L, Rao W, Zeng Y, Zhu L, He Y, He G. Identification and analysis of brown planthopper-responsive microRNAs in resistant and susceptible rice plants. Sci Rep. 2017;7:8712.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  30. Ge Y, Han J, Zhou G, Xu Y, Ding Y, Shi M, Guo C, Wu G. Silencing of miR156 confers enhanced resistance to brown planthopper in rice. Planta. 2018;248:813–26.

    Article  CAS  PubMed  Google Scholar 

  31. Dai Z, Tan J, Zhou C, Yang X, Yang F, Zhang S, Sun S, Miao X, Shi Z. The OsmiR396-OsGRF8-OsF3H-flavonoid pathway mediates resistance to the brown planthopper in rice (Oryza sativa). Plant Biotechnol J. 2019;17:1657–69.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Chen C, Ridzon DA, Broomer AJ, Zhou Z, Lee DH, Nguyen JT, Barbisin M, Xu NL, Mahuvakar VR, Andersen MR, et al. Real-time quantification of microRNAs by stem-loop RT-PCR. Nucleic Acids Res. 2005;33:e179.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  33. Hu J, Zhou J, Peng X, Xu H, Liu C, Du B, Yuan H, Zhu L, He G. The Bphi008a gene interacts with the ethylene pathway and transcriptionally regulates MAPK genes in the response of rice to brown planthopper feeding. Plant Physiol. 2011;156:856–72.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Dai X, Zhuang Z, Zhao PX. psRNATarget: a plant small RNA target analysis server (2017 release). Nucleic Acids Res. 2018;46:W49–54.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Sattar S, Tompson GA. Small RNA regulators of plant-hemipteran interactions: micromanagers with versatile roles. Front Plant Sci. 2016;7:1241.

    Article  PubMed  PubMed Central  Google Scholar 

  36. Manosalva PM, Davison RM, Liu B, Zhu X, Hulbert SH, Leung H, Leach JE. A germin-like protein gene family functions as a complex quantitative trait locus conferring broad-spectrum disease resistance in rice. Plant Physiol. 2009;149:286–96.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Liu Q, Yang J, Yan S, Zhang S, Zhao J, Wang W, Yang T, Wang X, Mao X, Dong J, et al. The germin-like protein OsGLP2-1 enhances resistance to fungal blast and bacterial blight in rice. Plant Mol Biol. 2016;92:411–23.

    Article  CAS  PubMed  Google Scholar 

  38. Jung HW, Kim KD, Hwang BK. Identification of pathogen-responsive regions in the promoter of a pepper lipid transfer protein gene (CALTPI) and the enhanced resistance of the CALTPI transgenic Arabidopsis against pathogen and environmental stresses. Planta. 2005;221:361–73.

    Article  CAS  PubMed  Google Scholar 

  39. Bethke G, Grundman RE, Sreekanta S, Truman W, Katagiri F, Glazebrook J. Arabidopsis PECTIN METHYLESTERASEs contribute to immunity against Pseudomonas syringae. Plant Physiol. 2014;164:1093–107.

    Article  CAS  PubMed  Google Scholar 

  40. Tan L, Showalter AM, Egelund J, Hernandez-Sanchez A, Doblin MS, Bacic A. Arabinogalactan-proteins and the research challenges for these enigmatic plant cell surface proteoglycans. Front Plant Sci. 2012;3:140.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Shirasu K. The HSP90-SGT1 chaperone complex for NLR immune sensors. Annu Rev Plant Biol. 2009;60:139–64.

    Article  CAS  PubMed  Google Scholar 

  42. Liu D, Chen X, Liu J, Ye J, Guo Z. The rice ERF transcription factor OsERF922 negatively regulates resistance to Magnaporthe oryzae and salt tolerance. J Exp Bot. 2012;63:3899–911.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Jacob F, Vernaldi S, Maekawa T. Evolution and conservation of plant NLR functions. Front Immunol. 2013;4:297–312.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  44. Li W, Zhu Z, Chern M, Yin J, Yang C, Ran L, Cheng M, He M, Wang K, Wang J, et al. A natural allele of a transcription factor in rice confers broad-spectrum blast resistance. Cell. 2017;170:114–26.

    Article  CAS  PubMed  Google Scholar 

  45. Gao M, Yin X, Yang W, Lam SM, Tong X, Liu J, Wang X, Li Q, Shui G, He Z. GDSL lipases modulate immunity through lipid homeostasis in rice. PLoS Pathog. 2017;13:e1006724.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  46. Macho AP, Zipfel C. Plant PRRs and the activation of innate immune signaling. Mol Cell. 2014;54:263–72.

    Article  CAS  PubMed  Google Scholar 

  47. Zhang J, Guan W, Huang C, Hu Y, Chen Y, Guo J, Zhou C, Chen R, Du B, Zhu L, et al. Combining next-generation sequencing and single-molecule sequencing to explore brown plant hopper responses to contrasting genotypes of japonica rice. BMC Genomics. 2019;20:682.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  48. Shangguan X, Zhang J, Liu B, Zhao Y, Wang H, Wang Z, Guo J, Rao W, Jing S, Guan W, et al. A mucin-like protein of planthopper is required for feeding and induces immunity response in plants. Plant Physiol. 2018;176:552–65.

    Article  CAS  PubMed  Google Scholar 

  49. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12:357–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. Anders S, Pyl PT, Huber W. HTSeq—A Python framework to work with high-throughput sequencing data. Bioinformatics. 2015;31:166–9.

    Article  CAS  PubMed  Google Scholar 

  51. Leng N, Dawson JA, Thomson JA, Ruotti V, Rissman AI, Smits BM, Haag JD, Gould MN, Stevart RM, Kendziorski C. EBSeq: an empirical Bayes hierarchical model for inference in RNA-seq experiments. Bioinformatics. 2013;29:1035–43.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Draghici S, Khatri P, Tarca AL, Amin K, Done A, Voichita C, Georgescu C, Romero R. A systems biology approach for pathway level analysis. Genome Res. 2007;17:1537–45.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. Schmittgen TD, Livak KJ. Analyzing real-time PCR data by the comparative C(T) method. Nat Protoc. 2008;3:1101–8.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

We would thank Novel Bioinformatics Ltd., Co. (Shanghai, China) for conducting the whole transcriptome sequencing. We would also like to thank reviewers for their careful reading and valuable suggestions.

Funding

This work was supported by grants from National Program on Research and Development of Transgenic Plants (2018ZX0800101B) and National Natural Science Foundation of China (31470102). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Author information

Author notes

  1. Jiaoyan Tan and Yan Wu contributed equally to this work.

Authors and Affiliations

  1. State Key Laboratory of Hybrid Rice, College of Life Sciences, Wuhan University, Wuhan, 430072, China

    Jiaoyan Tan, Yan Wu, Jianping Guo, Huimin Li, Lili Zhu, Rongzhi Chen, Guangcun He & Bo Du

  2. Institute for Biosciences and Biotechnology Research, University of Maryland, College Park, MD, 20850, USA

    Yan Wu

Authors
  1. Jiaoyan Tan
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Yan Wu
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Jianping Guo
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Huimin Li
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Lili Zhu
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Rongzhi Chen
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Guangcun He
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Bo Du
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

BD and GH designed and wrote the manuscript. BD, JT and HL participated in the experiments and data analysis. YW participated in the sample treatment, collection and sequencing. GJ developed the BPH6G line. RC and LZ prepared the experimental reagents. JT and YW contributed equally to this paper. All authors read and approved the manuscript.

Corresponding author

Correspondence to Bo Du.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Additional file 1: Table S1.

Summary of small RNA sequences. Total Reads: raw data after sequencing. Clean Reads: reads after the removal of adaptors, low-quality reads, and reads of < 20 nt and > 24 nt. Mapped Reads: clean reads mapped on the miRbase. S: WT; R: the BPH6G plants; 0, non-infested; early: early feeding stage; late: late feeding stage.

Additional file 2: Figure S1.

Size distribution and miRNAs annotation of the BPH6G and WT plants at non-infested, early and late feeding stages. A Length distribution of total reads. B Proportions of different kinds of small RNAs.

Additional file 3: Table S2.

DEMs in both varieties before and after BPH feeding.

Additional file 4: Table S3.

DEMs at the early and late feeding stages of the two varieties.

Additional file 5: Table S4.

Summary of mRNA expression libraries. Total Reads: raw data after sequencing. Clean Reads: reads after the removal of adaptors, low quality tags, and single-copy tags. Mapped Reads: clean reads mapped on the rice genome. S: WT; R: the BPH6G plants; 0, non-infested; early: early feeding stages; late: late feeding stages.

Additional file 6: Figure S2.

Hierarchical clustering analysis of DEGs of the BPH6G and WT plants after BPH feeding based on the log ratio of FPKM data. Red and green indicate upregulated and downregulated DEGs, respectively. Each row shows genes and each column represents a comparison.

Additional file 7: Figure S3.

Venn diagrams of the number of upregulated and downregulated DEGs (A), and opposite expression DEGs (B) of the BPH6G and WT plants at different feeding stages.

Additional file 8: Figure S4.

GO (Gene Ontology) analysis. Biological process, cellular component, and molecular function of up-(A) and down-regulated (B) DEGs in R_early/R0 and R_late/R0 respectively, and up- (C) and down-regulated (D) DEGs in S_early/S0 and S_late/S0 respectively (P < 0.05). The x-axis and y-axis indicate names of clusters and genes in a category, respectively.

Additional file 9: Table S5.

KEGG pathway enrichment analysis of DEGs appeared opposite expression at early or late feeding stages of two varieties.

Additional file 10: Table S6.

Integrated analysis of BPH resistance related miRNAs and their target genes.

Additional file 11: Table S7.

Primer sequences for qRT-PCR and the transformation of rice protoplasts.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Tan, J., Wu, Y., Guo, J. et al. A combined microRNA and transcriptome analyses illuminates the resistance response of rice against brown planthopper. BMC Genomics 21, 144 (2020). https://doi.org/10.1186/s12864-020-6556-6

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-020-6556-6

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us