Two rice lines were used in this study: Nipponbare (NIP) and BPH6-transgenic plants (BPH6G). NIP is a susceptible japonica line. BPH6G is a BPH6-transgenic line containing the BPH6 gene with its native Swarnalata promoter (IRRI Acc. No. 33964) in NIP background. The rice plants were grown in plastic cups (9 cm in diameter and 15 cm in height), in a greenhouse with 32 ± 2°C/14 h light and 26 ± 2°C/10 h dark periods. The BPH population was maintained on the susceptible rice cv. Taichung Native1 (IRRI Acc. No.00105) under controlled environmental conditions (as described above) at Wuhan University (Hu et al., 2017).

Seeds of either NIP or BPH6G were sown in plastic pots covered with nylon mesh. At the four-leaf stage, the rice seedlings were infested with third instar BPH nymphs at a rate of 10 nymphs per plant. Observation continued until susceptible control plants died, and then the rice plants were subsequently photographed and scored. At least three independent biological replicates were performed.

Measurements of honeydew excretion and BPH weight gain were performed following the method of Zheng et al. (Zheng et al., 2021). Twenty female BPH adults, which had been starved for 2 h, were weighed and introduced into a pre-weighed parafilm sachet (1.5 cm × 2.5 cm) fixed on the leaf sheath of each five-leaf stage NIP or BPH6G plant at a height of 2-3 cm above the soil. After 48 h of feeding, the surviving insects and the parafilm sachet were weighed again. The change in the insects’ weight was recorded as BPH weight gain, and the change in the sachets’ weight was recorded as honeydew excretion.

The endpoint method was utilized for sample collection throughout BPH treatment (Wu et al., 2017). Each treatment began at an individual time point and stopped at the same time. Seedlings at the four-leaf stage were infested with eight BPHs per seedling and collected after 0, 6, 12, 24, 48, 60, and 72 hours of infestation. Each treatment consisted of three biological replicates, with 15 seedlings per replicate. Leaf sheaths were mixed for non-infested controls (0 h), early infestation stage (6, 12, and 24 h), and late infestation stage (48, 60, and 72 h). BPH6G samples were annotated as R0, R_early, and R_late, and NIP samples were annotated as S0, S_early, and S_late. Leaf sheathes were excised, frozen in liquid nitrogen, and stored at -80°C until use.

RNA was isolated using an RNAiso Plus kit (TaKaRa), according to the manufacturer’s instructions. RNA quality was estimated using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific). An Illumina HiSeq 2500 was used for total RNA sequencing. The raw reads were first quality-controlled with FAST-QC by filtering low-quality reads (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Afterward, clean reads were aligned to the rice genome using Hisat2 software (Kim et al., 2015).

The mapped reads were assembled using StringTie (Pertea et al., 2015). All transcriptomes were merged to reconstruct a comprehensive transcriptome using Cuffmerge. Transcripts overlapping with known mRNAs and transcripts shorter than 200 bp were discarded. To identify lncRNAs, we utilized CPAT to predict coding transcripts, and transcripts with coding-prob score > 0.3 were removed (Wang et al., 2013). Transcript expression levels were quantified as fragments per kilobase of exon per million fragments mapped (FPKM). Differential expression P-values were calculated using the Bioconductor edgeR package (Robinson et al., 2010). We used the absolute value of log₂FC ≥ 1 and P < 0.05 as thresholds to judge the statistical significance of each differential expression result.

The sequencing data was used to predict circRNAs with the ACFS circRNA prediction pipeline (You and Conrad, 2016). Unmapped reads were obtained with BWA-MEM for circRNA identification. The head-to-tail junction was identified and the highest splicing strength score was calculated using MaxEntScan33, with a filtering criterion greater than or equal to 10. Based on the re-alignment of the unmapped reads to the circRNA candidates, reads which mapped to the circRNA back splicing junction (with an overhang of at least 6 nucleotides) were counted for each circRNA. Transcript expression levels were quantified as FPKM. Differential expression P-values were calculated using the Bioconductor edgeR package (Robinson et al., 2010). We used the absolute value of log₂FC ≥ 1 and P < 0.05 as thresholds to judge the statistical significance of each differential expression result.

Based on ceRNA theory (Salmena et al., 2011), we predicted the miRNA-mRNA, miRNA-lncRNA, and miRNA-circRNA interaction pairs using the PsTarget platform, with the standard E = 5 and UPE = 25 (http://plantgrn.noble.org/psRNATarget/analysis). Subsequently, through a combined analysis of miRNA and mRNA expression, lncRNA/circRNA-miRNA-mRNA pathways exhibiting either up-down-up or down-up-down expression modes were selected for further study. Cytoscape3.9.1 was used to display the networks (Shannon et al., 2003).

GO annotations from Gene Ontology (http://www.geneontology.org/) were downloaded (Ashburner et al., 2000). Pathway annotations were downloaded from KEGG (http://www.genome.jp/Kegg/). To identify significant GO and pathway categories, Fisher’s exact tests were applied under absolute values of P < 0.05 and FDR < 0.05 (Draghici et al., 2007).

First-strand cDNA was synthesized with a PrimeScript RT reagent kit (TaKaRa, AK2802) using 1 μg of total RNA. The cDNA was amplified by qRT-PCR using SYBR green supermix (Bio-Rad) and a CFX96 real-time system, according to the manufacturer’s instructions. All qRT-PCR primers used in this work are listed in Supplementary Table 1 . Three biological replicates were performed for each experiment. Normalized expression levels were calculated using the 2^{−ΔΔC (t)} method, with TBP as the internal reference gene.

The wild type (WT) and mutated (MUT) lncRNAs containing target miRNA binding sites were synthesized and cloned into the pGreenII 0800-miRNA reporter vector, yielding the lncRNA-WT and lncRNA-MUT vectors, respectively. The miRNA was amplified and cloned into the pCXUN overexpression vector, yielding the miRNA-OE vector. The miRNA-OE/lncRNA-WT, miRNA-OE/lncRNA-MUT, empty vector (PCXUN)/lncRNA-WT, and empty vector (PCXUN)/lncRNA-MUT pairs were co-transfected into rice protoplast cells, which were prepared following the method of Wu et al. (Wu et al., 2017). After 16 h of transfection, luciferase activities were evaluated with a dual-luciferase reporter assay system (Promega) and a SpectraMax iD5 multi-mode microplate reader (Molecular Devices). In this study, two lncRNA-miRNA interaction pairs (XLOC_042442-miR1846c and XLOC_028297-miR530) were selected to perform the assay ( Supplementary Table 1 ).

We evaluated the resistance of BPH6-transgenic plants (BPH6G) and susceptible Nipponbare plants (NIP) to BPH at the seedling stage, and the average damage severity score was calculated for each plant after infestation. By the 7th day of BPH infestation, the NIP plants died (average score of 9.0), while the BPH6G plants were still vigorously growing (average score of 2.8) ( Figures 1A, B ). To confirm these results, we measured the honeydew excretion and weight gain of BPHs allowed to feed on NIP and BPH6G plants for 48 h. BPHs on BPH6G plants excreted less honeydew and had a lower growth rate than BPHs on NIP plants ( Figures 1C, D ). These results demonstrate that BPH6G confers resistance to BPH infestation, while NIP remained susceptible.

To identify BPH-responsive non-coding RNAs in BPH6G (R) and NIP (S) plants, we obtained the FPKM values of lncRNAs and circRNAs from the whole-transcriptome RNA sequencing data after 0, 6, 12, 24, 48, 60, and 72 h of BPH infestation. The 0 h samples were taken as the non-infested controls, while the mixtures of 6, 12, and 24 h samples were taken as early-stage profiles and the mixtures of 48, 60, and 72 h samples were taken as late-stage profiles. In total, 6 treatment groups were analyzed: R0, R_early, R_late, S0, S_early, and S_late.

The raw reads from 18 rice transcriptomes were combined and filtered, resulting in the identification of 1219 lncRNAs. Most lncRNAs were less than 2,000 nt in length ( Figure 2A ), with the majority containing 2-4 exons, and 50.29% containing two exons ( Figure 2B ). The six treatment groups (R0, R_early, R_late, S0, S_early, and S_late) exhibited different lncRNA expression profiles ( Figure 2C ). The chromosomal distribution of the identified lncRNAs is shown in Figure 2D . Different chromosomes contained different numbers of lncRNAs, with chromosome 1 (Chr 1) containing the greatest number ( Figure 2D ).

The lncRNA expression levels were compared between treatments, and the differentially expressed lncRNAs (DElncRNAs) exhibiting absolute value of log₂FC ≥ 1 and P < 0.05 among the 6 treatment groups are shown in Figure 3A . The total number of DElncRNAs in group S (S_early/S_0 and S_late/S_0, 190) was higher than that of group R (R_early/R_0 and R_late/R_0, 108). During the late stage, the number of DElncRNAs in group S (S_late/S_0, 152) was 4 times greater than during the early stage (S_early/S_0, 38). The number of early and late DElncRNAs in group R was similar (R_early/R_0, 40; R_late/R_0, 68). There were 303 DElncRNAs identified in the three comparison groups among different varieties (R0/S0, R_early/S_early, and R_late/S_late) ( Figure 3A ), and the chromosomal distribution of the DElncRNAs was evaluated ( Figure 3B ). In order to identify lncRNAs related to BPH resistance in rice, we analyzed these DElncRNAs with a Venn diagram ( Figure 3C ). Overall, 60 DElncRNAs exhibited differential expression in R and S after BPH feeding (only in R_early/S_early and R_late/S_late). Thirty-two DElncRNAs were identified in R_early/R_0 which were absent in S_early/S_0 and 32 DElncRNAs were identified in R_late/R_0 which were absent in S_late/S_0. The expression profiles of these DElncRNAs were found for some modules ( Figure 3D ; Supplementary Tables 2 , 3 ). These results indicate that susceptible and resistant plants contained similar numbers of BPH-responsive lncRNAs during the early stage (38 in S_early/S_0 and 40 in R_early/R_0), although only 8 were identified in both groups ( Figure 3C ). Conversely, susceptible plants contained a greater number of BPH-responsive lncRNAs than resistant plants during the late stage ( Figures 3A, C ). That is, the two genotypes responded to BPH herbivory in entirely different manners. Subsequently, qRT-PCR was used to study the expression levels of these 6 DElncRNAs ( Figure 4A ). The data were consistent with the sequencing results, which confirmed that the results were reliable and could be used for intensive studies. Overall, the results confirmed that the numbers and types of BPH-responsive lncRNAs were different between resistant and susceptible plants, suggesting that lncRNAs are involved in the response of rice to BPH herbivory.

In plants, lncRNAs indirectly regulate mRNA expression through their regulation of miRNA (Bartel, 2009; Thomas et al., 2010). By combining the psRNA target tool results with the miRNA and mRNA high-throughput sequencing data (Tan et al., 2020), we predicted the presence of 15, 30, 20, and 83 DElncRNAs binding to 15, 35, 35, and 87 differentially expressed miRNAs (DEmiRNAs) in R_early/R0, R_late/R0, S_early/S0, and S_late/S0, respectively ( Supplementary Figure 1 ; Supplementary Tables 4 – 7 ). These DElncRNAs may affect the expression of their target differentially expressed mRNAs (DEmRNAs) (105, 242, 295, and 934 in R_early/R0, R_late/R0, S_early/S0, and S_late/S0, respectively) by sponging their corresponding DEmiRNAs ( Supplementary Tables 4 – 7 ). The DElncRNAs, DEmiRNAs, and DEmRNAs in the ceRNA network were compared using Venn diagrams ( Supplementary Figure 1 ).

To validate whether these DElncRNAs affected the expression of their corresponding DEmiRNAs, two DElncRNAs (XLOC_042442 and XLOC_028297) and their target DEmiRNAs (miR1846c and miR530) were selected for luciferase reporter assay verification. First, we constructed XLOC_042442 and XLOC_028297 WT or MUT (predicted miRNA binding sites) luciferase plasmids in the pGreenII 0800-miRNA dual luciferase reporter vector, and luciferase activity was evaluated after co-transfection of miRNA-containing (miR1846c and miR530, respectively) and luciferase plasmids. In rice protoplasts, we observed that overexpression of miR1846c and miR530 (miR1846c-OE and miR530-OE) reduced the luciferase activity significantly. At the same time, miRNA binding site mutations reversed the luciferase activity, suggesting that XLOC_042442 and XLOC_028297 function as ceRNAs by sponging miR1846c and miR530, respectively ( Figure 4B ).

The potential regulatory roles of the DElncRNAs via the ceRNA network were predicted by analyzing the functions of their target DEmRNAs through GO and KEGG pathway analysis ( Supplementary Figures 2 , 3 ). The GO annotations (P < 0.05) of the four groups of target DEmRNAs contained multiple biological processes, cellular components, and molecular functions ( Supplementary Figure 2 ). In BPH6G at the early feeding stage (R_early/R0), these GO terms were most significantly enriched in DNA replication initiation, plastid inner membrane, and syn-pimara-7,15-diene synthase activity. In NIP at the early feeding stage (S_early/S0), these GO terms were most significantly enriched in plant-type cell wall biogenesis, plasma membrane, and cellulose synthase activity. In BPH6G at the late feeding stage (R_late/R0), these GO terms were most significantly enriched in hydrogen peroxide catabolic process, integral component of membrane, and transmembrane transporter activity. In NIP at the late feeding stage (S_late/S0), these GO terms were most significantly enriched in regulation of hormone levels, plasma membrane, and ATP binding. KEGG pathway analysis was performed on the four groups of DEmRNAs, which showed that DNA replication and metabolic pathways were the most significantly enriched pathways ( Supplementary Figure 3 ).

In this study, we identified a total of 1914 circRNAs, and the genomic distribution of these circRNAs is shown in Figure 5A . Different chromosomes contained different numbers of circRNAs, with chromosome 1 (Chr 1) containing the greatest number ( Figure 5A ). In terms of length, the majority of circRNAs were between 300-600 bp ( Figure 5B ). The circRNA expression levels among the 6 treatment groups were visualized using violin plots ( Figure 5C ). GO annotation of the circRNA source genes showed that they were enriched in the following biological process terms: cellular component organization, organic substance metabolic process, and nitrogen compound metabolic process; the following cellular component terms: intracellular part, intracellular organelle, and intracellular organelle part; and the following molecular function terms: hydrolase activity, heterocyclic compound binding, and protein binding ( Figure 5D ). KEGG pathway analysis showed that fatty carbon metabolism, glycolysis/gluconeogenesis, RNA degradation, carbon fixation in photosynthetic organisms, and mRNA surveillance pathway were most significantly enriched ( Figure 5E ).

The circRNAs expression levels were compared between treatments, and we identified a total of 129 differentially expressed circRNAs (DEcircRNAs) exhibiting absolute value of log₂FC ≥ 1 and P < 0.05 among the 6 treatment groups ( Figure 6A ). Group S (S_early/S_0 and S_late/S_0, 110) contained 3 times more DEcircRNAs than group R (R_early/R_0 and R_late/R_0, 31), and the DEcircRNAs in group S were primarily down-regulated. The numbers of DEcircRNAs in the early and late stages of group R were similar (R_early/R_0, 12; R_late/R_0, 19), and less than that in group S (S_early/S_0, 54; S_late/S_0, 56). The numbers of up- and down-regulated DEcircRNAs were similar in group R (7/5 in R_early/R0, 8/11 in R_late/R0). There were 84 DEcircRNAs contained in the three comparison groups among different varieties (R0/S0, R_early/S_early, and R_late/S_late) ( Figure 6A ). The distribution and expression of DEcircRNAs were also mapped to the chromosomes of the two materials ( Figure 6B ).

In order to identify circRNAs related to BPH resistance in rice, we analyzed the DEcircRNAs between and within groups of resistant and susceptible plants before and after BPH feeding using Venn diagrams ( Figures 6C, D ). A total of 39 DEcircRNAs exhibited differential expression between R and S after BPH feeding (only in R_early/S_early and R_late/S_late). Ten DEcircRNAs were identified in R_early/R_0 which were absent in S_early/S_0, and 14 DEcircRNAs were identified in R_late/R_0 which were absent in S_late/S_0. The expression patterns of these DEcircRNAs are shown as a heat map ( Figure 6E ; Supplementary Tables 8 , 9 ). The quantitative detection results were consistent with the sequencing data, which confirmed the reliability of the sequencing results ( Figure 6F ). Overall, these results confirmed that the numbers and types of BPH-responsive circRNAs were different between resistant and susceptible plants, suggesting that circRNAs are involved in the response of rice to BPH herbivory.

By combining the psRNA target tool results with the miRNA and mRNA high-throughput sequencing data (Tan et al., 2020), we predicted the existence of 8, 10, 19, and 26 DEcircRNAs binding to 42, 61, 29, and 114 DEmiRNAs in the R_early/R0, R_late/R0, S_early/S0, and S_late/S0 comparisons, respectively ( Supplementary Tables 10 – 13 ). These DEcircRNAs may affect the expression of their target DEmRNAs (187, 401, 222, and 1090 in R_early/R0, R_late/R0, S_early/S0, and S_late/S0, respectively) by sponging their corresponding DEmiRNAs ( Supplementary Tables 10 – 13 ). The DEcircRNAs, DEmiRNAs, and DEmRNAs in the ceRNA network were compared by Venn diagrams ( Supplementary Figure 1 ).

The potential regulatory roles of the circRNAs via the ceRNA network were predicted by analyzing the functions of all the target DEmRNAs through GO and KEGG pathway analyses ( Supplementary Figures 4 , 5 ). The GO annotations (P < 0.05) of the four groups of DEmRNAs are listed in Supplementary Figure 3 and Supplementary Tables 10 – 13 , and include a multitude of biological processes, cellular components, and molecular functions. In BPH6G at the early feeding stage (R_early/R0), these GO terms were most significantly enriched in cell wall biogenesis, photosystem II antenna complex, and kinase activity. In NIP at the early feeding stage (S_early/S0), these GO terms were most significantly enriched in plant-type cell wall biogenesis, plasma membrane, and cellulose synthase (UDP-forming) activity. In BPH6G at the late feeding stage (R_late/R0), these GO terms were most significantly enriched in metabolic process, integral component of membrane, and cellulose synthase (UDP-forming) activity. In NIP at the late feeding stage (S_late/S0), these GO terms were most significantly enriched in cell wall organization, plasma membrane, and ATP binding. KEGG pathway analysis was performed on the four groups of DEmRNAs, which showed that DNA replication and metabolic pathways were the most significantly enriched pathways ( Supplementary Figure 4 ).

The lncRNA/circRNA-miRNA-mRNA interaction networks were predicted using psTarget software and consolidated using whole-transcriptome RNA sequencing data ( Figure 7A ; Supplementary Table 14 ). Specifically, the interaction network was established based on the relationship between the DElncRNAs/DEcircRNAs and DEmiRNAs and DEmRNAs. A total of 39 DElncRNAs were predicted to bind 50 DEmiRNAs and 381 target DEmRNAs, while 21 DEcircRNAs were predicted to bind 116 DEmiRNAs and 760 target DEmRNAs. The ceRNA network with 39 DElncRNAs, 21 DEcircRNAs, 133 DEmiRNAs, and 834 DEmRNAs is depicted in Figure 7A and Supplementary Table 14 .

The potential regulatory roles of the ceRNA networks were analyzed by characterizing the function of 834 target DEmRNAs through GO enrichment and KEGG pathway analyses. GO annotation showed that the network was primarily enriched in the following biological processes: cell wall organization or biogenesis, developmental growth, and single-organism cellular process; the following cellular components: intrinsic component of membrane, external encapsulating structure, and side of membrane; and the following molecular functions: transferase activity, signaling receptor activity, and carbohydrate binding ( Figure 7B ). In addition, the network was found to be enriched in the following KEGG pathways: metabolic pathways, biosynthesis of secondary metabolites, phenylpropanoid biosynthesis, starch and sucrose metabolism, and phenylalanine metabolism ( Figure 7C ). These results illustrated that both metabolism- and response-related genes are regulated by the ceRNA network in both resistant and susceptible plants in response to BPH herbivory.

Notably, the GO term “response to stimulus” (GO: 0006950) was significantly enriched in 80 DEmRNAs ( Supplementary Table 14 ). After consultation of relevant literature, several specific components of the lncRNA/circRNA-miRNA-mRNA interaction network were selected for further study: 6 DEcircRNAs, 15 DElncRNAs, 16 DEmiRNAs, and 17 DEmRNAs ( Figure 7D ). The lncRNAs XLOC_024786, XLOC_055548, and XLOC_043522, as well as the circRNAs circR9g19160 and circR11g36470, were predicted to combine with miR396f-3p, miR164c, and miR166k-5p. These three miRNA families have been described as regulators of rice innate immunity against fungi and bacteria (Raquel et al., 2018; Wang et al., 2018b; Li et al., 2019). In plants, R genes with the NB-ARC structure play vital roles in disease and pathogen resistance by recognizing specific effectors and inducing rapid and robust resistance responses (Jones and Dangl, 2006). In this study, we identified 6 NB-ARC domain-containing R proteins (MLO1-like, RPS2-like, RPM1-like, PIK2-like, RGA4-like (YP10), and RGA5-like (LP10)), which were negatively regulated by miR529b, miR399b, miR2886, miR2871b, miR396f-3p, and miR172c/miR166k-5p. These miRNAs were sponged by lncRNAs XLOC_011258, XLOC_002440, XLOC_016603/XLOC_064800, XLOC_011974/XLOC_041798, XLOC_024784, and XLOC_035677/XLOC_043522, respectively. Both CTL1/BC15 and CHIA, which encode membrane-associated chitinase-like proteins, have been reported to be involved in rice resistance (Wu et al., 2012), and are negatively regulated by the XLOC_059273/circR8g35050-miR1858a and XLOC_011258/circR7g38270-miR529b modules, respectively. Additionally, CIPK9, ATL1Q, and PDR12 have been reported to be involved in signal transduction and play significant roles in abiotic and biotic stress response (Guzmán, 2012; Kanwar et al., 2014; Nguyen et al., 2014). These genes were predicted to be the targets of miR397a, miR1846a/b/c, and miR164c, and these miRNAs were predicted to be regulated by circR7g38270, XLOC_042442, circR9g19160, circR11g36470, and XLOC_055548. Finally, PAL4 has been associated with broad-spectrum disease resistance in rice (Tonnessen et al., 2015), and is negatively regulated by the XLOC_010703-miR2090 and XLOC_020472-miR2090 modules. These examples imply that ceRNAs might play important regulatory roles in the resistance of rice to BPH herbivory.

In a previous study, an integrated analysis of the miRNAs and target mRNA genes related to BPH resistance in BPH6-transgenic plants and NIP was performed (Tan et al., 2020). A network of 34 miRNAs corresponding to 42 target genes was identified, which were potentially related to BPH resistance (Tan et al., 2020). Meanwhile, considering our ceRNA results, we obtained a core ceRNA network consisting of 6 lncRNAs, 4 circRNAs, 23 miRNAs, and 24 mRNAs ( Figure 8 ). This core ceRNA network might also play a vital role in the resistance of rice to BPH herbivory.