Comprehensive identification and characterization of lncRNAs and circRNAs reveal potential brown planthopper-responsive ceRNA networks in rice

Brown planthopper (Nilaparvata lugens Stål, BPH) is one of the most destructive pests of rice. Non-coding RNA plays an important regulatory role in various biological processes. However, comprehensive identification and characterization of long non-coding RNAs (lncRNAs) and circular RNAs (circRNAs) in BPH-infested rice have not been performed. Here, we performed a genome-wide analysis of lncRNAs and circRNAs in BPH6-transgenic (resistant, BPH6G) and Nipponbare (susceptible, NIP) rice plants before and after BPH feeding (early and late stage) via deep RNA-sequencing. A total of 310 lncRNAs and 129 circRNAs were found to be differentially expressed. To reveal the different responses of resistant and susceptible rice to BPH herbivory, the potential functions of these lncRNAs and circRNAs as competitive endogenous RNAs (ceRNAs) were predicted and investigated using Gene Ontology and Kyoto Encyclopedia of Genes and Genomes analyses. Dual-luciferase reporter assays revealed that miR1846c and miR530 were targeted by the lncRNAs XLOC_042442 and XLOC_028297, respectively. In responsive to BPH infestation, 39 lncRNAs and 21 circRNAs were predicted to combine with 133 common miRNAs and compete for miRNA binding sites with 834 mRNAs. These mRNAs predictably participated in cell wall organization or biogenesis, developmental growth, single-organism cellular process, and the response to stress. This study comprehensively identified and characterized lncRNAs and circRNAs, and integrated their potential ceRNA functions, to reveal the rice BPH-resistance network. These results lay a foundation for further study on the functions of lncRNAs and circRNAs in the rice-BPH interaction, and enriched our understanding of the BPH-resistance response in rice.

Introduction

Rice, one of the most important food crops worldwide, is host to more than 200 insect pests at different stages of its life cycle. Among these, the brown planthopper (Nilaparvata lugens Stål, BPH) is one of the most destructive. Recent surveys suggest that BPH damages rice production severely, accounting for 29.5% of the total rice crop loss due to insects and diseases and making BPH the number one pest of Chinese rice (Du et al., 2020). Therefore, the development and production of BPH-resistant rice varieties is an economical, effective, safe, and environmentally-friendly strategy to control BPH damage (Yu et al., 2022).

In order to improve the resistance of rice to BPH, the genes and signaling pathways pivotal to BPH resistance must be identified and characterized. Over the last five decades, around 40 BPH resistance genes (R genes) have been identified, and the mechanisms of rice resistance to BPH have been explored through genetic and biochemical analyses (Chen et al., 2022; Guo et al., 2023). Recently, the novel BPH R gene BPH6, encoding an exocyst-localized protein, was cloned and found to confer broad-spectrum resistance to planthoppers (Guo et al., 2018). Specifically, the BPH6 gene promotes exocytosis and participates in the maintenance and reinforcement of plant cell wall (Wu et al., 2022). Further functional characterization of BPH6 will more precisely illuminate the exact mechanism underpinning BPH6-mediated herbivory resistance.

microRNAs (miRNAs) have been demonstrated to regulate rice resistance to various pathogens and herbivores by modulating the expression of target genes at the post-transcriptional level (Kar and Raichaudhuri, 2021). For example, both miR159 and miR160 regulate the defense response against blast disease (Li et al., 2014; Chen et al., 2021), and miR159b, miR164a, and miR167d-5p are involved in the plant immune response to bacterial blight (Jia et al., 2020). These miRNAs fine-tune plant innate immunity through the integration of R gene expression, phytohormone signaling, callose deposition, and reactive oxygen species (ROS) production (Kumar et al., 2022). Two miRNAs (miR156 and miR396) have been reported to regulate rice resistance to BPH. Specifically, miR156 regulates jasmonic acid (JA) and jasmonoyl-isoleucine (JA-Ile) biosynthesis through the “miR156-OsMPK3/6-OsWRKY70” module, thus negatively regulating BPH resistance (Ge et al., 2018). miR396 also negatively regulates BPH resistance through the “miR396-OsGRF8-OsF3H-flavonoid” module (Dai et al., 2019).

Recent studies revealed that long non-coding RNAs (lncRNAs) and circular RNAs (circRNAs) may serve as competing endogenous RNAs (ceRNAs), which could be essential for regulating the circuitry of miRNAs and their target genes (Salmena et al., 2011). Specifically, lncRNAs and circRNAs can modulate the balance between miRNAs and target genes. lncRNAs are eukaryotic non-coding RNAs greater than 200 nucleotides (nt) in length (Kapranov et al., 2007; Ponjavic et al., 2007). They have been found to regulate both biotic and abiotic stress tolerance, as well as various growth and developmental processes in rice (Gao et al., 2020). Non-coding endogenous circRNA molecules are covalently closed continuous loops without 5’-3’ polarity or a polyadenylated tail and have been found to respond to abiotic and biotic stimuli and growth processes (Jeck et al., 2013). However, the specific regulatory mechanism underlying lncRNA- and circRNA-mediated BPH resistance in rice remains to be elucidated.

In this study, we performed a genome-wide analysis of lncRNAs and circRNAs in BPH6-transgenic (resistant, BPH6G) and Nipponbare (susceptible, NIP) rice plants before and after BPH herbivory via deep RNA-sequencing. We also conducted an integrated analysis of the expression profiles of circRNAs, lncRNAs, and previously identified miRNAs and mRNAs (Tan et al., 2020). lncRNA/circRNA-miRNA-mRNA ceRNA networks were generated by combining the identified and annotated target mRNAs. Our results demonstrated that lncRNAs and circRNAs act as ceRNAs to regulate BPH resistance in rice. This study provides a foundation for further research into the molecular mechanisms underlying BPH6-conferred herbivory resistance in rice.

Materials and methods

Plant and BPH materials

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

Evaluation of rice resistance to BPH

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.

Sample collection

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.

Total RNA extraction and high-throughput sequencing

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

Identification of lncRNAs and differential expression analysis

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.

Identification of circRNAs and differential expression analysis

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.

Construction and analysis of the ceRNA regulatory network

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

Analysis of Gene Ontology and Kyoto Encyclopedia of Genes and Genomes pathways

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

Quantitative real-time PCR assay

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.

Luciferase reporter assays

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

Results

lncRNA and circRNA expression in BPH6G and NIP before and after BPH feeding

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.

FIGURE 1

Figure 1 Evaluation of BPH resistance in BPH6G and NIP plants. (A) BPH resistance phenotypes of BPH6G and NIP plants. (B) BPH resistance scores of BPH6G and NIP plants. Lower scores correspond to greater resistance. Data represent the means ± SE of 3 biologically independent experiments. (C) Honeydew excretion of BPHs on BPH6G and NIP plants after 2 d of feeding. Data represent the means ± SE of 10 replicates. (D) Weight gain of BPHs after 2 d of feeding on BPH6G and NIP plants. Data represent the means ± SE of 10 replicates. Data were analyzed by ANOVA and asterisks indicate statistically significant differences. **P < 0.01.

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.

Analysis of lncRNA characteristics and response to BPH invasion

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

FIGURE 2

Figure 2 Identification and distribution of all lncRNAs. (A) Lengths of lncRNAs. (B) Number of exons. (C) FPKM distribution of lncRNAs in six groups. (D) Chromosomal distribution of all lncRNAs.

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.

FIGURE 3

Figure 3 Differentially expressed lncRNAs. (A) Number of DElncRNAs up- or down-regulated in all comparisons. (B) Genomic distribution of all DElncRNAs. The two circles (from outer to inner) represented the expression levels (log₂FC) of DElncRNAs in the resistant plants (red indicates increased expression, and green denotes decreased expression), and in the susceptible plants (orange indicates increased expression, and blue denotes decreased expression), and fold change of the expression levels, respectively. (C) Venn diagrams of the unique and shared DElncRNAs. (D) Heat map of the DElncRNAs in differential comparisons, yellow indicates N/A.

FIGURE 4

Figure 4 Detection of candidate lncRNAs expression levels and validation of target miRNAs. (A) DElncRNAs validated by qRT-PCR. (B) The luciferase reporter assay illustrates that candidate DElncRNAs could serve as a sponge for predicted target DEmiRNAs. Data are presented as the means ± SE of 3 biologically independent experiments. Data were analyzed by ANOVA and asterisks indicate statistically significant differences. *P < 0.05, **P < 0.01.

Analysis of DElncRNA function via ceRNA regulatory networks in response to BPH feeding

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

Analysis of circRNA characteristics and response to BPH

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

FIGURE 5

Figure 5 Distribution of circRNAs and GO and KEGG analyses of source genes. (A) Chromosomal distribution of circRNAs. (B) Lengths of circRNAs. (C) FPKM distribution of circRNAs in six groups. (D) GO enrichment analysis of circRNA source genes. (E) KEGG pathway enrichment analysis of circRNA source genes.

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

FIGURE 6

Figure 6 Differentially expressed circRNAs. (A) Number of DEcircRNAs up- or down-regulated in all comparisons. (B) Genomic distribution of all DEcircRNAs. The two circles (from outer to inner) represented the expression levels (log₁₀FPKM) of DEcircRNAs in the resistant plants (red indicates increased expression, and green denotes decreased expression), DEcircRNAs in the susceptible plants (orange indicates increased expression, and blue denotes decreased expression), and fold change of the expression levels, respectively. (C-D) Venn diagrams of the unique and shared DEcircRNAs. (E) Heat map of the unique DEcircRNAs in differential comparisons. (F) Validation of candidate DEcircRNAs expression levels by qRT-PCR. circR1g43630, chr1_25000624_24999798_+826-LOC_Os01g43630; circR1g38229, chr1_21426198_21424262_+1936-LOC_Os01g38229; circR3g16220, chr3_8951244_8940833_-10411-LOC_Os03g16220; circR2G55370, chr2_33928616_33914366_+14250-LOC_Os02g55370. Data were analyzed by ANOVA and asterisks indicate statistically significant differences. **P < 0.01.

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

Analysis of key ceRNA pathways responding to BPH infestation

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.

FIGURE 7

Figure 7 Regulatory network of potential BPH resistance-related lncRNAs and circRNAs. (A) CeRNA networks of potential BPH resistance-related DElncRNAs/DEcircRNAs-DEmiRNAs. Different colored circles indicate different non-coding RNA types: red represents DEcircRNAs, green represents DElncRNAs, and blue represents DEmiRNAs. (B) GO enrichment of 834 target DEmRNAs regulated by the ceRNA network. (C) KEGG enrichment of 834 target DEmRNAs regulated by the ceRNA network. (D) Core ceRNA networks of lncRNAs/circRNAs-miRNAs-mRNAs were significantly enriched in GO term “response to stimulus” (GO: 0006950). circR7g38270^a, chr7_22980701_22946863_-33838-LOC_Os07g38270; circR7g38270^b, chr7_22976110_22956434_-19676-LOC_Os07g38270; circR8g35050, chr8_22086695_22054445_-32250-LOC_Os08g35050; circR7g49270, chr7_29510790_29508498_+2292-LOC_Os07g49270; circR9g19160, chr9_11675871_11394316_+281555-LOC_Os09g19160; circR11g36470, chr11_21913121_20045605_-1867516-LOC_Os11g36470.

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.

FIGURE 8

Figure 8 CeRNA networks of lncRNA/circRNA-miRNA-mRNA related to plant resistance. Yellow indicates N/A. circR7g38270^a, chr7_22980701_22946863_-33838-LOC_Os07g38270; circR8g35050, chr8_22086695_22054445_-32250-LOC_Os08g35050; circR7g38270^b, chr7_22976110_22956434_-19676-LOC_Os07g38270; circR1g02450, chr1_795720_785269_-10451-LOC_Os01g02450.

Discussion

In rice, both lncRNAs and circRNAs have been reported to regulate plant growth and development, as well as the biotic and abiotic stress responses (Sanchita and Asif, 2020; Zhang et al., 2020). To date, several lncRNAs have been functionally characterized in rice. For example, LDMAR (long-day-specific male fertility-associated RNA) is a 1236bp photoperiod sensitive lncRNA which regulates male sterility in rice. The normal development of plant pollen under long-day conditions requires a high level of LDMAR expression (Fan et al., 2016). The expression of the lncRNA IPS1 (induced by phosphate starvation 1) is altered under nitrogen and phosphorus starvation, suggesting that IPS1 may regulate the balance of these macronutrients (Shin et al., 2018). Overexpression of the lncRNA LAIR (leucine-rich repeat receptor kinase antisense intergenic RNA) in rice increases grain yield and induces the expression of several leucine-rich repeat receptor kinase-coding genes (Wang et al., 2018a). Although there has been some research into plant circRNAs, currently-reported mechanisms of circRNA formation are largely based on bioinformatics analyses rather than convincing experimental evidence, which is still in the theoretical stage (Zhang et al., 2020). To the best of our knowledge, this is the first study to comprehensively identify and characterize BPH-responsive lncRNAs and circRNAs associated with ceRNA networks in rice. Therefore, the results presented here further our understanding of non-coding RNA regulation of the rice response to BPH infestation.

In this study, we identified BPH-responsive lncRNAs and circRNAs, and their differential expression patterns were analyzed in BPH6-transgenic and NIP rice plants. Our genome-wide analysis of lncRNA and circRNA expression showed that susceptible NIP plants contained greater numbers of both DElncRNAs and DEcircRNAs than resistant BPH6-transgenic plants, and that the majority of these DElncRNAs and DEcircRNAs were down-regulated in NIP plants. The functions of these lncRNAs and circRNAs via the ceRNA network were also predicted by GO and KEGG analyses of target mRNAs (Supplementary Figures 2–5). The GO and KEGG pathway results differed between the resistant and susceptible varieties at both the early and late feeding stages (Supplementary Figures 2–5). Our results suggest that the BPH-responsive lncRNAs and circRNAs in BPH6G and NIP were involved in different ceRNA pathways.

The interactions between different RNAs were predicted based on ceRNA theory, and multiple BPH-resistance related networks were identified and included 39 lncRNAs, 21 DEcircRNAs, 133 DEmiRNAs, and 834 DEmRNAs (Figure 7A; Supplementary Table 14). The functions of the ceRNA networks were then explored through GO and KEGG pathway analyses (Figures 7B, C). Among the ceRNA networks, 80 DEmRNAs were predicted to be involved in the response to stress (GO: 0006950) (Supplemental Table 15). A core ceRNA network of these DEmRNAs is shown in Figure 7D. Among them, OsmiR396 has been identified as a negative regulator of rice innate immunity against Magnaporthe oryzae by silencing multiple growth-regulating factors (OsGRFs) in rice (Chandran et al., 2019). Overexpression of OsmiR396 makes rice susceptible to blast, while suppression of OsmiR396 makes rice resistant to blast and improves yield (Chandran et al., 2019). In addition, overexpression of OsmiR396f resulted in enhanced immunity to Dickeya zeae by suppressing the target gene OsGRF (Li et al., 2019). A previous study showed that OsmiR396 also acts as a negative regulator of BPH resistance via the OsmiR396-OsGRF8-OsF3H-flavonoid pathway (Dai et al., 2019). Our results showed that OsmiR396f-3p, which targets the OsWRKY59 transcription factor, is regulated by lncRNA XLOC_024786. Hence, we speculate that the XLOC_024786-OsmiR396f-3p-OsWRKY59 module may play a vital role in rice resistance to BPH herbivory (Figure 7D). OsmiR164a has been reported to negatively regulate rice immunity to M. oryzae and Xanthomonas oryzae pv. oryzae by regulating the OsNAC60 transcription factor (Wang et al., 2018b; Jia et al., 2020). We also found that lncRNA XLOC_055549 could sponge miR164a, suggesting that the XLOC_055549-OsmiR164a module may be a key regulator of BPH resistance (Figure 7D). The miR166 family is highly conserved (Kumar et al., 2022). The miR166 family members OsmiR166k and OsmiR166h have been shown to function as positive regulators of defense against M. oryzae and Fusarium fujikuroi by targeting the ethylene insensitive 2 (EIN 2) gene via cross-regulation (Raquel et al., 2018). This study demonstrated that lncRNA XLOC_043522 negatively regulates the expression of miR166k-5p, thus up-regulating the expression of the R gene YR10 (Figure 7D). The ceRNA networks identified in this study are likely to play important roles in rice resistance to BPH. By integrating and analyzing the results of our ceRNA study, alongside the results of Tan et al.’s miRNA and target mRNA study, a core ceRNA network containing lncRNAs, circRNAs, miRNAs, and mRNAs was produced (Figure 8) (Tan et al., 2020). These results provide evidence of novel regulatory mechanisms underlying rice BPH resistance. Our future research will focus on elucidating of the BPH-responsive ceRNA network in rice, which will greatly increase our knowledge of plant resistance to insect damage.

The rice-BPH system is considered an excellent model for the study of plant-insect interaction and co-evolution (Jing et al., 2017). Xiao et al. provided evidence that lncRNAs might play important roles in the high fecundity and virulence adaptation of BPH (Xiao et al., 2015). A previous study has reported the differential expression of lncRNAs between two virulent BPH populations, including susceptible (TN1) and resistant (YHY15) rice planthopper varieties (Zha et al., 2022). In total, 157 differentially expressed lncRNAs and 675 differentially expressed mRNAs were found to be involved in BPH adaptation to rice resistance (Zha et al., 2022). Here, we studied the differential expression of lncRNAs and circRNAs in rice before and after BPH herbivory. Our analysis indicated that both lncRNAs and circRNAs play important roles in rice resistance to BPH. These results suggest the presence of mutual regulatory relationships between rice and BPH at the non-coding RNA level, and provide a basis for further studies of the molecular mechanism underlying co-evolution between rice and BPH.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/, GSE123148.

Author contributions

AY, BD, and YW conceived and designed the experiments. YW, WZ, and DQ performed the experiments and analyzed the data. YW wrote the paper. JG developed the BPH6G line. GL, CL, BW, SL, JC, LH, and SS contributed reagents, materials, and analysis tools. LH, LZ, and ZZ revised the paper. YW, WZ, and DQ contributed equally to this paper. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the Science and Technology Major Program of Hubei Province (2022ABA001), the Youth Project of Natural Science Foundation of Hubei Province (2022CFB747), the Youth Science Foundation Project of Hubei Academy of Agricultural Sciences (2023NKYJJ03), and funding from the Hubei Key Laboratory of Food Crop Germplasm and Genetic Improvement (2021lzjj04).

Acknowledgments

We would thank Novel Bioinformatics Ltd., Co. (Shanghai, China) for conducting the whole transcriptome sequencing. And we would also like to thank TopEdit (www.topeditsci.com) for its linguistic assistance during the preparation of this manuscript.

Conflict of interest

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

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2023.1242089/full#supplementary-material

Supplementary Figure 1 | Venn diagrams of the DElncRNAs, DEcircRNAs, DEmiRNAs, and DEmRNAs in the ceRNA network. (A–C). Venn diagrams of DElncRNAs, DEmiRNAs, and DEmRNAs in the ceRNA network of lncRNA-miRNA-mRNA. (D–F). Venn diagrams of DEcircRNAs, DEmiRNAs, and DEmRNAs in the ceRNA network of circRNA-miRNA-mRNA.

Supplementary Figure 2 | GO enrichment of the target genes of DElncRNAs via the ceRNA network at the early and late feeding stages in susceptible and resistant plants.

Supplementary Figure 3 | KEGG pathway analysis of the target genes of DElncRNAs via the ceRNA network at the early and late feeding stages in susceptible and resistant plants. The red color represents the pathways with P < 0.05.

Supplementary Figure 4 | GO enrichment of the target genes of DEcircRNAs via the ceRNA network at the early and late feeding stages in susceptible and resistant plants.

Supplementary Figure 5 | KEGG pathway analysis of the target genes of DEcircRNAs via the ceRNA network at the early and late feeding stages in susceptible and resistant plants. The red color represents the pathways with P < 0.05.

References

Ashburner, M., Ball, C. A., Blake, J. A., Botstein, D., Butler, H., Cherry, J. M., et al. (2000). Gene ontology: tool for the unification of biology. Nat. Genet. 25 (1), 25–29. doi: 10.1038/75556

PubMed Abstract | CrossRef Full Text | Google Scholar

Bartel, D. P. (2009). MicroRNAs: target recognition and regulatory functions. Cell 136 (2), 215–233. doi: 10.1016/j.cell.2009.01.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Chandran, V., Wang, H., Gao, F., Cao, X. L., Chen, Y. P., Li, G. B., et al. (2019). miR396-osGRFs module balances growth and rice blast disease-resistance. Front. Plant Sci. 9. doi: 10.3389/fpls.2018.01999

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, R., Deng, Y., Ding, Y., Guo, J., Qiu, J., Wang, B., et al. (2022). Rice functional genomics: decades’ efforts and roads ahead. Sci. China Life Sci. 65 (1), 33–92. doi: 10.1007/s11427-021-2024-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, J. F., Zhao, Z. X., Li, Y., Li, T. T., Wang, W. M. J. R. (2021). Fine-tuning roles of osa-miR159a in rice immunity against Magnaporthe oryzae and development. Rice 14 (1), 1–11. doi: 10.1186/s12284-021-00469-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Dai, Z., Tan, J., Zhou, C., Yang, X., Yang, F., Zhang, S., et al. (2019). The OsmiR396-OsGRF8-OsF3H-flavonoid pathway mediates resistance to the brown planthopper in rice (Oryza sativa). Plant Biotechnol. J. 17 (8), 1657–1669. doi: 10.1111/pbi.13091

PubMed Abstract | CrossRef Full Text | Google Scholar

Draghici, S., Khatri, P., Tarca, A. L., Amin, K., Done, A., Voichita, C., et al. (2007). A systems biology approach for pathway level analysis. Genome Res. 17 (10), 1537–1545. doi: 10.1101/gr.6202607

PubMed Abstract | CrossRef Full Text | Google Scholar

Du, B., Chen, R., Guo, J., He, G. (2020). Current understanding of the genomic, genetic, and molecular control of insect resistance in rice. Mol. Breed. 40 (2), 1–25. doi: 10.1007/s11032-020-1103-3

CrossRef Full Text | Google Scholar

Fan, Y., Yang, J., Mathioni, S. M., Yu, J., Shen, J., Yang, X., et al. (2016). PMS1T, producing phased small-interfering RNAs, regulates photoperiod-sensitive male sterility in rice. Proc. Natl. Acad. Sci. U.S.A. 113 (52), 15144–15149. doi: 10.1073/pnas.1619159114

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, C., Zheng, X., Li, H., Mlekwa, U. A., Yu, G., Jie, X. (2020). Roles of lncRNAs in rice: advances and challenges. Rice Sci. 27 (5), 384–395. doi: 10.1016/j.rsci.2020.03.003

CrossRef Full Text | Google Scholar

Ge, Y., Han, J., Zhou, G., Xu, Y., Ding, Y., Shi, M., et al. (2018). Silencing of miR156 confers enhanced resistance to brown planthopper in rice. Planta 248, 813–826. doi: 10.1007/s00425-018-2942-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, J., Wang, H., Guan, W., Guo, Q., Wang, J., Yang, J., et al. (2023). A tripartite rheostat controls self-regulated host plant resistance to insects. Nature 618, 799–807. doi: 10.1038/s41586-023-06197-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, J., Xu, C., Wu, D., Zhao, Y., Qiu, Y., Wang, X., et al. (2018). Bph6 encodes an exocyst-localized protein and confers broad resistance to planthoppers in rice. Nat. Genet. 50 (2), 297–306. doi: 10.1038/s41588-018-0039-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Guzmán, P. (2012). The prolific ATL family of RING-H2 ubiquitin ligases. Plant Signaling Behav. 7 (8), 1014–1021. doi: 10.4161/psb.20851

CrossRef Full Text | Google Scholar

Hu, L., Wu, Y., Wu, D., Rao, W., Guo, J., Ma, Y., et al. (2017). The coiled-coil and nucleotide binding domains of BROWN PLANTHOPPER RESISTANCE14 function in signaling and resistance against planthopper in rice. Plant Cell 29 (12), 3157–3185. doi: 10.1105/tpc.17.00263

PubMed Abstract | CrossRef Full Text | Google Scholar

Jeck, W. R., Sorrentino, J. A., Wang, K., Slevin, M. K., Burd, C. E., Liu, J., et al. (2013). Circular RNAs are abundant, conserved, and associated with ALU repeats. RNA 19 (2), 141–157. doi: 10.1261/rna.035667.112

PubMed Abstract | CrossRef Full Text | Google Scholar

Jia, Y., Li, C., Li, Q., Liu, P., Liu, D., Liu, Z., et al. (2020). Characteristic Dissection of Xanthomonas oryzae pv. oryzae Responsive MicroRNAs in Rice. Int. J. Mol. Sci. 21 (3), 785. doi: 10.3390/ijms21030785

PubMed Abstract | CrossRef Full Text | Google Scholar

Jing, S., Zhao, Y., Du, B., Chen, R., Zhu, L., He, G. (2017). Genomics of interaction between the brown planthopper and rice. Curr. Opin. Insect Sci. 19, 82–87. doi: 10.1016/j.cois.2017.03.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Jones, J. D., Dangl, J. L. (2006). The plant immune system. Nature 444 (7117), 323–329. doi: 10.1038/nature05286

PubMed Abstract | CrossRef Full Text | Google Scholar

Kanwar, P., Sanyal, S. K., Tokas, I., Yadav, A. K., Pandey, A., Kapoor, S., et al. (2014). Comprehensive structural, interaction and expression analysis of CBL and CIPK complement during abiotic stresses and development in rice. Cell Calcium 56 (2), 81–95. doi: 10.1016/j.ceca.2014.05.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Kapranov, P., Cheng, J., Dike, S., Nix, D. A., Duttagupta, R., Willingham, A. T., et al. (2007). RNA maps reveal new RNA classes and a possible function for pervasive transcription. Science 316 (5830), 1484–1488. doi: 10.1126/science.113834

PubMed Abstract | CrossRef Full Text | Google Scholar

Kar, M. M., Raichaudhuri, A. (2021). Role of microRNAs in mediating biotic and abiotic stress in plants. Plant Gene 26, 100277. doi: 10.1016/j.plgene.2021.100277

CrossRef Full Text | Google Scholar

Kim, D., Langmead, B., Salzberg, S. L. (2015). HISAT: a fast spliced aligner with low memory requirements. Nat. Methods 12 (4), 357–360. doi: 10.1038/nmeth.3317

PubMed Abstract | CrossRef Full Text | Google Scholar

Kumar, K., Mandal, S. N., Neelam, K., De Los Reyes, B. G. (2022). MicroRNA-mediated host defense mechanisms against pathogens and herbivores in rice: balancing gains from genetic resistance with trade-offs to productivity potential. BMC Plant Biol. 22 (1), 351. doi: 10.1186/s12870-022-03723-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, W., Jia, Y., Liu, F., Wang, F., Fan, F., Wang, J., et al. (2019). Integration analysis of small RNA and degradome sequencing reveals microRNAs responsive to Dickeya zeae in resistant rice. Int. J. Mol. Sci. 20 (1), 222. doi: 10.3390/ijms20010222

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Y., Lu, Y.-G., Shi, Y., Wu, L., Xu, Y.-J., Huang, F., et al. (2014). Multiple rice microRNAs are involved in immunity against the blast fungus Magnaporthe oryzae. Plant Physiol. 164 (2), 1077–1092. doi: 10.1104/pp.113.230052

PubMed Abstract | CrossRef Full Text | Google Scholar

Nguyen, V. N. T., Moon, S., Jung, K. H. (2014). Genome-wide expression analysis of rice ABC transporter family across spatio-temporal samples and in response to abiotic stresses. J. Plant Physiol. 171 (14), 1276–1288. doi: 10.1016/j.jplph.2014.05.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Pertea, M., Pertea, G. M., Antonescu, C. M., Chang, T.-C., Mendell, J. T., Salzberg, S. L. (2015). StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat. Biotechnol. 33 (3), 290–295. doi: 10.1038/nbt.3122

PubMed Abstract | CrossRef Full Text | Google Scholar

Ponjavic, J., Ponting, C. P., Lunter, G. (2007). Functionality or transcriptional noise? Evidence for selection within long noncoding RNAs. Genome Res. 17 (5), 556–565. doi: 10.1101/gr.6036807

PubMed Abstract | CrossRef Full Text | Google Scholar

Raquel, S. G., Hsing, Y., San, S. B. (2018). The Polycistronic miR166k-166h Positively Regulates Rice Immunity via Post-transcriptional Control of EIN2. Front. Plant Sci. 9, 337. doi: 10.3389/fpls.2018.00337

PubMed Abstract | CrossRef Full Text | Google Scholar

Robinson, M. D., McCarthy, D. J., Smyth, G. K. (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26 (1), 139–140. doi: 10.1093/bioinformatics/btp616

PubMed Abstract | CrossRef Full Text | Google Scholar

Salmena, L., Poliseno, L., Tay, Y., Kats, L., Pandolfi, P. P. (2011). A ceRNA hypothesis: the Rosetta Stone of a hidden RNA language? Cell 146 (3), 353–358. doi: 10.1016/j.cell.2011.07.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Sanchita, T. P.K., Asif, M. H. (2020). Updates on plant long non-coding RNAs (lncRNAs): the regulatory components. Plant Cell Tissue Organ Culture 140 (2), 259–269. doi: 10.1007/s11240-019-01726-z

CrossRef Full Text | Google Scholar

Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13 (11), 2498–2504. doi: 10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

Shin, S. Y., Jeong, J. S., Lim, J. Y., Kim, T., Park, J. H., Kim, J. K., et al. (2018). Transcriptomic analyses of rice (Oryza sativa) genes and non-coding RNAs under nitrogen starvation using multiple omics technologies. BMC Genomics 19 (1), 532. doi: 10.1186/s12864-018-4897-1

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Thomas, M., Lieberman, J., Lal, A. (2010). Desperately seeking microRNA targets. Nat. Struct. Mol. Biol. 17 (10), 1169–1174. doi: 10.1038/nsmb.1921

PubMed Abstract | CrossRef Full Text | Google Scholar

Tonnessen, B. W., Manosalva, P., Lang, J. M., Baraoidan, M., Bordeos, A., Mauleon, R., et al. (2015). Rice phenylalanine ammonia-lyase gene OsPAL4 is associated with broad spectrum disease resistance. Plant Mol. Biol. 87 (3), 273–286. doi: 10.1007/s11103-014-0275-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y., Luo, X., Sun, F., Hu, J., Zha, X., Su, W., et al. (2018a). Overexpressing lncRNA LAIR increases grain yield and regulates neighbouring gene cluster expression in rice. Nat. Commun. 9 (1), 3516. doi: 10.1038/s41467-018-05829-7

CrossRef Full Text | Google Scholar

Wang, L., Park, H. J., Dasari, S., Wang, S., Kocher, J.-P., Li, W. (2013). CPAT: Coding-Potential Assessment Tool using an alignment-free logistic regression model. Nucleic Acids Res. 41 (6), e74–e74. doi: 10.1093/nar/gkt006

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Z., Xia, Y., Lin, S., Wang, Y., Guo, B., Song, X., et al. (2018b). Osa-miR164a targets OsNAC60 and negatively regulates rice immunity against the blast fungus Magnaporthe oryzae. Plant J. 95 (4), 584–597. doi: 10.1111/tpj.13972

CrossRef Full Text | Google Scholar

Wu, D., Guo, J., Zhang, Q., Shi, S., Guan, W., Zhou, C., et al. (2022). Necessity of rice resistance to planthoppers for OsEXO70H3 regulating SAMSL excretion and lignin deposition in cell walls. New Phytol. 234 (3), 1031–1046. doi: 10.1111/nph.18012

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, Y., Lv, W., Hu, L., Rao, W., Zeng, Y., Zhu, L., et al. (2017). Identification and analysis of brown planthopper-responsive microRNAs in resistant and susceptible rice plants. Sci. Rep. 7 (1), 8712. doi: 10.1038/s41598-017-09143-y

CrossRef Full Text | Google Scholar

Wu, B., Zhang, B., Dai, Y., Zhang, L., Shang-Guan, K., Peng, Y., et al. (2012). Brittle culm15 encodes a membrane-associated chitinase-like protein required for cellulose biosynthesis in rice. Plant Physiol. 159 (4), 1440–1452. doi: 10.1104/pp.112.195529

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiao, H., Yuan, Z., Guo, D., Hou, B., Yin, C., Zhang, W., et al. (2015). Genome-wide identification of long noncoding RNA genes and their potential association with fecundity and virulence in rice brown planthopper, Nilaparvata lugens. BMC Genomics 16, 749. doi: 10.1186/s12864-015-1953-y

PubMed Abstract | CrossRef Full Text | Google Scholar

You, X., Conrad, T. O. (2016). Acfs: accurate circRNA identification and quantification from RNA-Seq data. Sci. Rep. 6 (1), 1–11. doi: 10.1038/srep38820

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, S., Ali, J., Zhou, S., Ren, G., Xie, H., Xu, J., et al. (2022). From Green Super Rice to green agriculture: Reaping the promise of functional genomics research. Mol. Plant 15 (1), 9–26. doi: 10.1016/j.molp.2021.12.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Zha, W., Li, S., Xu, H., Chen, J., Liu, K., Li, P., et al. (2022). Genome-wide identification of long non-coding (lncRNA) in Nilaparvata lugens’s adaptability to resistant rice. PeerJ 10, e13587. doi: 10.7717/peerj.13587

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, P., Li, S., Chen, M. (2020). Characterization and function of circular RNAs in plants. Front. Mol. Biosci. 7, 91. doi: 10.3389/fmolb.2020.00091

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, X., Xin, Y., Peng, Y., Shan, J., Zhang, N., Wu, D., et al. (2021). Lipidomic analyses reveal enhanced lipolysis in planthoppers feeding on resistant host plants. Sci. China Life Sci. 64 (9), 1502–1521. doi: 10.1007/s11427-020-1834-9

PubMed Abstract | CrossRef Full Text | Google Scholar