The NILs containing either BPH14 (B14), BPH15 (B15), or both BPH14/BPH15 (B1415) genes were developed using inbred indica rice variety ‘Wushansimiao’, as the RP (Wang et al., 2019). Seeds were planted in plastic cups (15 cm high by 9 cm wide) at a density of 15 plants per cup, and greenhouse-grown under a 10 h dark (26 ± 2°C)/14 h light (32 ± 2°C) cycle. The BPHs were maintained at Wuhan University, China, on ‘Taichung Native1’ (TN1; susceptible cultivar, IRRI Acc. No.00105) under environmental conditions identical to those of the rice plants.

BPH nymphs (third instar) were introduced at a rate of 8 BPH per seedling to four-leaf stage B14, B15, B1415, and RP seedlings. As described previously (Huang et al., 2001), seedlings were ascribed a resistance score during examination. The average damage severity score (0, 1, 3, 5, 7, or 9) was calculated for each plant after infestation.

Pre-weighed parafilm sachets were used to confine starved (2 h) third instar BPH nymphs and fastened to the leaf sheathes of one-month-old B14, B15, B1415, and RP plants (Pathak et al., 1982). The sachets were removed and emptied of BPH insects after 2 d of active herbivory. All sachets were weighed post BPH removal, and the weight difference before and after 2 d of herbivory was recorded as the amount of honeydew excretion.

Both BPH treatment and sample collection were accomplished according to the endpoint method (Wu et al., 2017). All treatments ended at the same time, despite beginning at different times. After 0, 3, 6, 12, 24, 48, and 72 h, four-leaf stage B14, B15, B1415, and RP seedlings were infested at a rate of 8 BPH nymphs (third instar) per seedling. Each experiment consisted of three biological replicates per treatment, with each replicate containing 15 seedlings. Leaf sheath samples were designated as either the non-infested group (0 h), early infestation group (3, 6, and 12 h), or late infestation group (24, 48, and 72 h). The experimental sample designations were as follows: B14_0, B14_early, and B14_late for the B14 lines; B15_0, B15_early, and B15_late for the B15 lines; B1415_0, B1415_early, B1415_late for the B1415 lines; and RP_0, RP_early, and RP_late for the RP lines. All samples were frozen with liquid N₂ and stored at -80°C prior to analyses.

Total RNA was collected from leaf sheathes with Trizol (Invitrogen). Quality was established with a Bioanalyzer 2200 (Aligent). All samples were stored at -80°C prior to analyses.

A TruSeq Stranded mRNA Library Prep Kit (Illumina) was utilized for preparation of the cDNA libraries, according to the standard protocol. Briefly, oligo (dT) magnetic beads were utilized to purify poly-A mRNA from 1 μg total RNA, which was then fragmented (200-600 bp) for 6 min with divalent cations (85°C). Both first- and second-strand cDNA synthesis were carried out using the cleaved RNA fragments. dUTP mix was utilized for second-strand cDNA synthesis, allowing for second strand separation. The cDNA fragments were then ligated with indexed adapters, A-tailed, and end-repaired. To remove the second-strand cDNA, the ligated cDNA was purified and subjected to uracil DNA glycosylase. The cDNA libraries were created by using PCR to enrich the purified first-strand cDNA. An Agilent 2200 was used for library quality control, and the libraries were sequenced using NovaSeq 6000 on a 150 bp paired-end run.

Adapter sequences and low-quality reads were removed in order to acquire clean reads. Hisat2 was utilized to align the clean reads with the reference genome (IRGSP1.0, Ensembl) (Kim et al., 2015). Gene counts were acquired with HTseq. Gene expression was quantified according to the fragments per kilo base of exon per million fragments mapped (FPKM) (Anders et al., 2015).

DEGs were filtered using the DESeq2 algorithm (Love et al., 2014). Statistically significant DEGs were determined according to P-value (< 0.05), fold change (FC; log₂FC > 1 or log₂FC < -1), and FDR (< 0.05) (Benjamini et al., 2001). Here, DEGs are defined as transcripts exhibiting a P-value < 0.05 and at least a 2-fold change in FPKM (log₂FC > 1 or log₂FC < -1).

GO evaluation was carried out to elucidate the biological importance of the identified DEGs (Ashburner et al., 2000), using GO annotations downloaded from the Gene Ontology (http://www.geneontology.org/), UniProt (http://www.uniprot.org/), and NCBI (http://www.ncbi.nlm.nih.gov/) databases. Statistically significant GO categories were determined with the Fisher’s exact test (P-value < 0.05).

KEGG pathway evaluation was carried out to determine the biological pathways associated with the identified DEGs according to the KEGG database. Statistically significant KEGG pathways were determined with the Fisher’s exact test (P-value < 0.05) (Draghici et al., 2007).

A PrimeScript RT Reagent Kit containing gDNA Eraser (RR047A, TaKaRa) was used to convert total RNA into first-strand cDNA. qRT-PCR was accomplished on a CFX96 real-time system (Bio-Rad) with SYBR Green Real-Time PCR Master Mix (QPK-201, Toyobo). All primers are listed in Supplementary Table 1 . Gene expression was evaluated by relative quantification, with TBP as the endogenous reference (Livak and Schmittgen, 2001).

The NB domain of BPH14 and the OsPOX8.1 coding sequence were amplified from ‘B5’ and ‘Wushansimiao’ cDNAs, and then respectively cloned into the ZeBaTA-based pCXUN expression vector with a Myc tag at the c-terminus (Chen et al., 2009). All primers are listed in Supplementary Table 1 . The aforementioned constructs were transiently transfected into 10-day old rice stem protoplasts as described previously (Chen et al., 2006).

Transfected protoplasts were extracted using a protein extraction buffer containing 5 mM MgCl₂, 100 mM Tris-HCl (pH 7.5), 0.5% (w/v) Triton X-100, and 1 mM EDTA, with 1 mM PMSF and 2 mM DTT included just prior to the assay. Total soluble proteins were collected from rice protoplasts (5×10⁶ cells per sample) using 200 μL of extraction buffer. SDS-PAGE was carried out to separate 10 μL of the extract. The extract was diluted (1:1000) with dilution buffer (3% [w/v] BSA, 150 mM NaCl, 0.1% [w/v] Tween 20, 20 mM Tris-HCl [pH 7.4]) and used for anti-Myc antibody (M192-3, MBL) immunoblotting, and subsequently incubated with 5% (w/v) skim milk-diluted (1:10,000) secondary antibody conjugated to horseradish peroxidase (115-035-003, Jackson). Detection was carried out using Tanon high-sig ECL protein gel blotting substrate.

All experiments consisted of three biological replicates, except where stated otherwise. Equivalent results were obtained using three independent biological experiments. Statistically significant differences were identified using Student’s t-tests at P value < 0.05.

Evaluation of ROS production in rice protoplasts (as shown in Figure 1C ) was carried out with a modified chemiluminescence method (Zhang et al., 2007). Briefly, the protoplasts were transfected for 16-22 h and then quantified and diluted to 1×10⁵ cells/200 μL with W5. To the diluted protoplasts was added 20 μM of the luminol derivative 8-amino-5-chloro-7-phenylpyrido [3,4-d] pyridazine-1,4 (2H,3H) dione (L-012) (Wako) and 20 μg/mL horseradish peroxidase (Sigma-Aldrich). Luminescence was captured using a SpectraMax iD5 multi-mode microplate reader (Molecular Devices).

In this study, four-leaf stage NILs containing either the BPH14 (B14), BPH15 (B15), or BPH14/BPH15 genes (B1415), as well as their RP, were infested with BPH. RP plants began to wither after 4 d of BPH herbivory (average score of 4.7), and wilted completely after 7 d (average score of 8.2). However, the B14, B15 and B1415 plants showed no visible damage (average scores of 3.3, 2.0, and 1.5, respectively) and survived until the end of the experiment (average scores of 4.6, 5.6, and 3.3, respectively) ( Figures 2A–D ).

To investigate the antibiosis effects of the NIL and RP plants, we measured the quantity of BPH-secreted honeydew. Overall, BPH feeding on RP and NIL plants produced very little honeydew from 3 to 6 h after infestation. Interestingly, the most significant differences in honeydew production were observed at 12 h after infestation, with the amount of honeydew production remaining relatively constant on NIL plants (from a minimum of 0.16 mg at 6 h to a maximum of 0.36 mg at 12 h) and increasing sharply on RP plants (from a minimum of 0.21 mg at 6 h to a maximum of 14.4 mg at 12 h). After 12 h, honeydew production increased on both NIL and RP plants, and remained high from 24 to 72 h after infestation ( Figure 2E ).

To identify DEGs functionally associated with defense against BPH in NILs pyramiding the BPH14 and BPH15 genes, RNA was extracted from the leaf sheaths of B14, B15, B1415, and RP plants after infestation (0-72 h). Samples were grouped as non-infested (0 h), early feeding stage (3, 6 and 12 h), or late feeding stage (24, 48 and 72 h) for RNA-seq.

Differences in BPH-responsive gene expression were analyzed using mRNA libraries. From 36 mRNA libraries, a total of 31,739,768 to 50,072,060 reads were sequenced. After removing low quality sequences, 82.38%-89.76%, 82.96%-90.53%, 81.69%-90.48%, and 87.83%-90.93% of the reads were mapped to 25,781,469-42,511,732 (RP), 32,708,227-43,047,757 (B14), 29,297,676-41,357,463 (B15), and 33,097,294-44,922,512 (B1415) rice genes, respectively ( Supplementary Table 2 ).

Subsequent to normalization, the average normalized reads from three independent biological replicates were selected for further studies. In total, 14,492 DEGs were identified among 17 comparisons, including nine comparisons among the different varieties (B14_0/RP_0, B14_early/RP_early, B14_late/RP_late, B15_0/RP_0, B15_early/RP_early, B15_late/RP_late, B1415_0/RP_0, B1415_early/RP_early, B1415_late/RP_late) and eight comparisons among the different feeding stages (RP_early/RP_0, RP_late/RP_0, B14_early/B14_0, B14_late/B14_0, B15_early/B15_0, B15_late/B15_0, B1415_early/B1415_0, B1415_late/B1415_0) ( Figure 3 and Supplementary Table 3 ).

A total of 690, 1,388, and 1,400 DEGs were identified in the B14_0/RP_0, B14_early/RP_early, and B14_late/RP_late comparisons, respectively; a total of 6,150, 3,235, and 709 DEGs were identified in the B15_0/RP_0, B15_early/RP_early, and B15_late/RP_late comparisons, respectively; and a total of 960, 1,247, and 520 DEGs were identified in the B1415_0/RP_0, B1415_early/RP_early, and B1415_late/RP_late comparisons, respectively ( Figure 3 ). In addition, 6,165 DEGs were identified in RP plants (1,644 in RP_early/RP_0 and 4,521 in RP_late/RP_0), 7,959 DEGs were identified in B14 plants (889 in B14_early/B14_0 and 7,070 in B14_late/B14_0), 7,784 DEGs were identified in B15 plants (941 in B15_early/B15_0 and 6,843 in B15_late/B15_0), and 4,234 DEGs were identified in B1415 plants (1,224 in B1415_early/B1415_0 and 3,010 in B1415_late/B1415_0). These results illustrate that the DEGs were responsive to BPH feeding, with a higher response in B14 and B15 plants compared to RP, and a lower response in B1415 plants than in RP plants ( Figure 3 ). These results suggest the presence of different modes of regulation at the early and late herbivory stages among the different rice plants.

During interactions between host plants and herbivores, reference gene expression is often suppressed (Hu et al., 2011). Normalization candidates were chosen after identifying which of the following common rice reference genes were the most stably-expressed: RPS27α (Os01g0328400), ACTIN1 (Os03g0718100), β-tubulin (Os03g0780600), eEF1α (Os03g0177500), GAPDH (Os02g0601300), SDHA (Os07g0134800), HSP (Os03g0426900), LSD1 (Os12g0611000), TBP (Os03g0657000), and Ubiquitin (Os03g0131300). Each was evaluated using FPKM values extracted from the RNA-seq data. Overall, both RPS27α and ACTIN1 expressions were significantly reduced after BPH herbivory in all groups. Compared with other candidates, TBP exhibited the most stable and appropriate expression level and was chosen as the endogenous reference gene for qRT-PCR validation assays ( Figure 4A ). The expression levels of eight DEGs were determined by qRT-PCR utilizing gene-specific primers ( Supplementary Table 1 ) for RNA-seq verification, and we found that the data were in agreement ( Figure 4B ).

To discover BPH resistance-associated genes, DEGs appearing in the comparisons of the resistant NIL vs. RP plants were analyzed by Venn diagrams, respectively ( Figures 5A–C ). We identified 150 overlapping DEGs in the B14_0/RP_0, B15_0/RP_0, and B1415_0/RP_0 comparisons ( Figure 5A ); 267 overlapping DEGs in the B14_early/RP_early, B15_early/RP_early, and B1415_early/RP_early comparisons ( Figure 5B ); and 218 overlapping DEGs in the B14_late/RP_late, B15_late/RP_late, and B1415_late/RP_late comparisons ( Figure 5C ). By combining these overlapping results, we obtained 531 DEGs common to B14, B15, and B1415 plants before and after BPH feeding, which may be involved in BPH resistance ( Supplementary Table 4 ).

To functionally categorize these 531 DEGs, we analyzed their associated GO terms and KEGG pathways. The DEGs were mainly enriched in the defense response, phosphorylation, and salt stress response GO biological processes; in the ATP binding, nucleotide binding, and kinase activity GO molecular functions; and the integral component of membrane, vacuole, and vacuolar membrane GO cellular components ( Figure 5D ). For KEGG analysis, the BPH-responsive DEGs were found to be primarily enriched in alpha-linolenic acid metabolism, amino sugar and nucleotide sugar metabolism, fatty acid metabolism, fatty acid degradation, and monoterpenoid biosynthesis. ( Figure 5E ).

Finally, we comprehensively evaluated both the expression patterns of and the relevant research reports pertaining to the identified DEGs, and ultimately landed on 11 BPH resistance-related genes. Among these, nine were significantly up-regulated in the resistant NIL plants compared with RP plants before and after BPH herbivory, while two (Os02g0599500 and Os10g0180800) were down-regulated in the resistant NIL plants compared with RP plants before and after BPH herbivory ( Figure 5F ).

The DEGs of both resistant NIL and RP plants at the early and late stages of herbivory were compared with those at the non-infested stage using Venn diagrams ( Figures 6A, B ). A total of 31 DEGs were specifically expressed in B14_early/B14_0, B15_early/B15_0, and B1415_early/B1415_0, while a total of 228 DEGs were specifically expressed in B14_late/B14_0, B15_late/B15_0, and B1415_late/B1415_0 ( Figures 6A, B ). The common DEGs were further pooled, and 258 DEGs were found to be shared only in resistant NILs during the early or late stages of BPH feeding ( Supplementary Table 5 ).

To functionally categorize these 258 DEGs, we analyzed their associated GO and KEGG pathways. The DEGs were mainly enriched in the response to oxidative stress, chloroplast organization, and response to karrikin GO biological processes; the metal ion binding, transferase activity, and glucosyltransferase activity GO molecular functions; and the cytosol, cell wall, and extracellular region GO cellular components ( Figure 6C ). For KEGG analysis, the BPH-responsive DEGs were found to be primarily enriched in betalain biosynthesis, biosynthesis of secondary metabolites, metabolic pathways, phagosome, and phenylpropanoid biosynthesis ( Figure 6D ).

Finally, we comprehensively evaluated both the expression patterns of and the relevant research reports pertaining to the identified DEGs, and ultimately landed on 10 BPH resistance-related genes. Most of these candidates, excluding Os06g0341300, were rapidly up- or down-regulated during the early BPH feeding stage, specifically in resistant NIL plants, with significant differences remaining during the late herbivory stages ( Figure 6E ).

Defense against BPH and other pathogens often involves the generation of ROS (Hu et al., 2017). Specifically, the rapid accumulation of ROS serves as a signal that coordinates an astonishing diversity of defense processes, while also being directly toxic to intruders (Gechev et al., 2006). The OsPOX8.1 gene, encoding a class III peroxidase, is highly up-regulated in response to blast and bacterial blight, where it is involved in the generation of ROS (Yin et al., 2000; Xiao et al., 2015). Here, the BPH-responsive candidate gene OsPOX8.1 was found to belong to the GO category “response to oxidative stress” ( Figure 6C ). Through qRT-PCR validation, we found that OsPOX8.1 was significantly up-regulated by BPH herbivory only in B14, B15, and B1415 plants, and the degree of up-regulation was higher in B1415 plants than in B14 or B15 plants. Consistent with the RNA-seq results, OsPOX8.1 expression was responsive from the early through the late feeding stages ( Figure 1A ).

To verify whether OsPOX8.1 regulates ROS levels, rice protoplasts were first transformed with OsPOX8.1. An empty vector construct (control, CK) and an auto-activated construct of the NB domain of BPH14 (BPH14-NB) were utilized as negative and positive controls, respectively (Hu et al., 2017). According to the immunoblotting experiments, each of the transformed proteins exhibited expected expression patterns ( Figure 1B ). ROS production in the protoplast lines was measured histochemically using the chemiluminescence method. The protoplasts transformed with OsPOX8.1 exhibited ROS accumulation, which was significantly stronger than that of protoplasts transformed with CK, but weaker than that of the protoplasts transformed with BPH14-NB ( Figure 1C ). These results indicated that OsPOX8.1 enhanced ROS production in rice protoplasts.