Kangwenqingzhan (KW), an Indica rice variety developed by Guangdong Academy of Agricultural Sciences, with no resistance genes to BPH, is highly susceptible to brown planthopper (BPH) and KW-Bph36-NIL (NIL), a near-isogenic line, introgressed BPH-resistant gene Bph36 into KW were used as study materials. Rice seeds were germinated in a 10 cm diameter plastic cup and cultured in a greenhouse at 16/8 h day/night photoperiod under 28/25°C day/night temperature. When the rice seedlings reached 3 leaf stage, weak seedlings were removed, leaving 10 plants per cup. Then, put the transparent gauze net bag on the cup and released the 2^nd-3^rd instar larvae of BPH, ensuring 8 nymphs per seedling on average. Samples were taken at 0h (control), 6h, 24h, and 48h after BPH onset, with three replicates for each treatment. Resistant materials were named as NIL0, NIL6, NIL24, and NIL48 and susceptible samples as KW0, KW6, KW24, and KW48. BPH nymphs were released at different time points, and all materials were sampled at the same time point. Samples were taken from about 5 cm of rice stem, which was eaten by BPH, and immediately put in liquid nitrogen.

Three biological replicates were used for all RNA-seq experiments sampled from each time point. Total RNA was extracted from each sample using TRIzol Reagent (Invitrogen)/RNeasy Mini Kit (Qiagen). Total RNA was quantified and qualified by Agilent 2100 Bio-analyzer (Agilent Technologies, Palo Alto, CA, USA), Nano Drop (Thermo Fisher Scientific Inc.) and 1% agarose gel. RNA with RIN value above 7 was used for further sequencing library construction. Preserved samples were sent to UniqueGene Company (Wuhan) for transcriptome sequencing. RNA-seq library was constructed by Illumina sequencing using the removed rRNA method with Ribo-Zero™GoldKits.

Rice genomic sequence (O. sativa_323_v7.0) and the corresponding annotations were retrieved from the Phytozome database (Version 12.0) (43). After filtering low-quality raw reads with adaptor sequences, clean reads were obtained with FASTX-Toolkit (version 0.0.14) and the purified data with high quality were mapped to the reference genome using HiSAT2 (Pertea et al., 2016). At last, the mapped reads of each sample were assembled by StringTie (version 1.3.3b) in a reference-based approach (Pertea et al., 2016). The value of fragment per kb per million reads (FPKM) was calculated to estimate the expression level of genes. For biological replicates, transcripts with |log₂Ratio|≥1 and q<0.05 were designated differentially expressed between samples (6/24/48 h) for KW or NIL, determining up-regulated and down-regulated genes.

Based on transcript classification codes, “u” (unknown) intergenic transcripts were regarded as novel gene loci and used for lncRNAs identification. Five steps were adopted for the identification of bona fide lncRNAs as previously described (Wang et al., 2015): (1) transcripts with length ≥ 200 bp and detected in more than 3 samples; (2) transcripts derived from rRNA and tRNA were removed (cutoff E-value0.001); (3) transcripts encoding proteins and protein-coding domains were removed by the search against the Swiss-Prot and Pfam databases (cutoff E-value 0.001); (4) OrfPredictor was applied to predict ORFs and transcripts that encode more than 100 amino acids were removed (Min et al., 2005). (5) Transcripts were removed that did not pass the protein-coding-score test using the Coding-Non-Coding Index (CNCI) (Luo et al., 2014), Coding Potential Calculator (CPC) (Kong et al., 2007), and Pfam-scan analysis (Mistry et al., 2007). Transcripts without coding potential were retained as novel lncRNAs, which were used for further analysis.

The resulting lncRNA transcripts and known transcripts were then merged into non-redundant transcripts, which were further quantified by StringTie for each sample (Trapnell et al., 2012). Differential expression analysis for each sequenced library was performed using ballgown (Pertea et al., 2016). The corrected P value of 0.05 and abs |log2 (Fold change)| of 1 was set as the threshold for significant differential expression between samples (6/24/48 h) for KW or NIL, determining up-regulated and down-regulated DELs. Singular Enrichment Analysis from AgriGO was performed to identify significantly enriched GO terms in the gene list out of the background of the reference gene list (Du et al., 2010). GO term pathways with a false discovery rate (q-value) < 0.05 were considered as significantly altered. Pearson correlation was employed to explore the expression relationship between lncRNAs and their neighboring genes (≤ 10 Kb). To identify lncRNAs and putative target genes, the TAPIR tool was used with the default settings (Bonnet et al., 2010). The relationship between lncRNAs and genes was used to construct the interaction networks with Cytoscape (version 3.9.1) (Saito et al., 2012).

To validate RNA-seq data, we randomly selected four genes based on their expression from DEGs and tested their expression using RT-qPCR. Total RNA used earlier for sequencing was reversely transcribed using PrimeScript™ RT reagent Kit with gDNA Eraser (Perfect Real Time) (Cat. RR047). qRT-PCR reactions were performed on a LightCycler 480 using TB Green™ Premix Ex Taq™II(Tli RNaseH Plus) (Code No. RR820). The reaction conditions were as follows: 95°C for 1 min followed by 40 cycles at 95°C for 15 s, 60°C for 15s, and 72°C for 28s. All RT-qPCR reactions were performed in triplicates. The results of each reaction were analyzed based on the 2^-ΔΔCt method (Livak and Schmittgen, 2001). The detection of threshold cycle for each reaction was normalized against the expression level of the rice Actin 1 gene.

Brown planthopper (BPH) is the most destructive insect pest (IP) of rice, causing direct and indirect damage to the crops. Extensive measures like chemical and biological have been taken for the control of the pest during the past decades which alternatively resulted in increased BPH resistance and risk to our environment. Conversely, equipping host plants with resistant genetic resources is the most effective and environment-friendly approach to managing BPH. Therefore, the present study uses KW as a recurrent parent and the constructed near-isogenic lines KW-Bph36-NIL to evaluate the high-quality and complete transcriptome of rice in response to BPH. First, KW and KW-Bph36-NIL (referred to as NIL) were used to perform the resistance identification to BPH feeding. As shown in (Supplementary figure 1), KW was highly sensitive while NIL exhibited high resistance to BPH suggesting the involvement of NIL materials has a vital role in BPH-resistant responses. These results further illustrate the facilitation of the development of BPH-resistant rice varieties.

In order to validate the findings of these transcripts to determine the reliability of transcriptome in gene expression, we randomly checked six genes that were influenced by BPH treatments in the two strains. Among them, three genes (Os10g0150300; Os06g0547933; Os02g0816200) were up-regulated, and three genes (Os03g0760500; Os03g0807500; Os03g0821150) were down-regulated in terms of relative expression in both KW and NIL in our RNA-seq dataset. Therefore, we further strengthen the reliability of our results by RT-qPCR analysis. Expression patterns of all the examined genes were in agreement with RNA-seq data, further indicating the credibility of our transcriptome dataset for gene exploration ( Supplementary Figure 2 ). The sequences of primer are listed ( Supplementary Table 1 ). The results showed that the expression trends of transcripts in both analyses were consistent indicating reliability over experimental data.

Since we observed significant differences in the physiological performance of rice cultivars to BPH feeding, we attempted to further emphasize the alterations in the two cultivars. Therefore, we performed transcriptome high-throughput Illumina sequencing of treatments i.e., 0h, 6h, 24h, and 48h with 3 biological replicates for KW and NIL after BPH onset. In order to understand the dynamic responses of rice to BPH feeding in detail, we defined 6h, 24h, and 48h as the early stage, middle stage, and late stage of rice defense against BPH and 0h as control. Total reads and mapped reads for each replicate were estimated ( Supplementary Figure 3 and Supplementary Table 2 ). The quality and completeness of the transcriptome can have a substantial impact on annotation and downstream analyses and any miscalculation in the transcriptome assembly could affect the prediction analysis, phylogenetic signaling, and gene expression qualifications. Thus, we analyzed the quality of filtered clean reads and found that the Q30 of all samples was over 95%, indicating that the sequencing quality met the analytical conditions ( Supplementary Figure 4 ). A total of 2994 DEGs were annotated in rice transcriptome after BPH feeding ( Supplementary Table 3 ). Hence, it would be of great interest to find out the influence of BPH feeding on our RNA-seq dataset. Therefore, we compared the number of DEGs and other expressed genes in each cultivar against the control. The ratio of DEGs to other expressed genes was 441 (1.48%), 680 (2.27%), and 210 (0.70%) for 6h, 24h, and 48h, respectively, under BPH treatments for KW ( Figure 1A , Supplementary Table 4 ). In comparison to KW, NIL has almost the same proportion of DEGs expressed for 6h, 331 (1.09%), and 24h, 218 (0.71%), while a surprising increase was observed for the number of DEGs at 48h treatment 1958 (6.44%) ( Figure 1B and Supplementary Table 5 ). These results indicated that DEGs were less in number compared to non-altered genes and further implicated that the BPH feeding has impacted the GE level mostly at 48h treatment in NIL. We further ensured the specificity of DEGs between two cultivars by comparing specific and shared genes. Only a limited number of shared DEGs were identified which reflects different responses of rice cultivars toward BPH feeding ( Figure 1C ). When attacked by BPH, similar to the changes in KW, the number of DEGs increased in NIL. The total number of DEGs in NIL was higher than that in KW which indicated that the resistant and susceptible rice cultivars responded differently to the dynamic changes of BPH feeding.

We also compared the DEGs of KW at different time points to evaluate specific alterations at specific time points of the BPH feeding. We found 164 specific genes for 6h feeding, 350 for 24h, and 87 genes for 48h feeding respectively. There were also 48 genes shared between 3-time points of the BPH feeding ( Figure 1D ). Next was the assessment of all DEGs of KW into upregulated and downregulated genes to identify whether their regulation into transcription is affected or not. In this regard, we observed an exceeded number of downregulated genes than upregulated genes indicating suppression of some mechanisms involved in stress tolerance ( Figures 1E, F ). In the same context, DEGs of NIL at different time points were also compared to identify specific and shared genes. NIL exhibited 67 expressed genes for 6h, 22 genes for 24h, and 1638 genes for 48h feeding treatments ( Figure 1G ). Assessment of these DEGs into upregulated and downregulated transcripts indicated an unusual expression pattern as compared to KW treatments. The number of upregulated genes was high in the case of NIL for 6h and 24h treatments, while there was a significant increase in the number of downregulated genes observed for DEGs of 48h treatment ( Figures 1H, I ). These differences observed in the expression pattern of DEGs indicated some special responses of NIL to BPH feeding, especially at 48h treatment.

More DEGs changed dramatically in the NIL strain indicating its stronger resistance than KW. Based on these results, we may conclude that NIL activated different regulatory mechanisms compared to KW. This phenomenon showed that the two strains had different resistance patterns against BPH, NIL was mainly through positive regulation, while KW was more through negative regulation.

Hence we observed significant differences in the number of DEGs among the two strains; we wonder whether there will be differences among DEGs at different treatment times. For this purpose, we assessed combined DEGs expressed at different time points and found little differences among genes at various time points ( Figure 2A and Supplementary Table 6 ). But, following the expression pattern of NIL for 48h, here we also observed a significant increase in the number of DEGs for 48h treatment ( Figure 2B ), which illustrates that BPH feeding has a significant impact on the number of differential genes expression at 48h. Next, we compared the regulation patterns of genes to find out up and downregulated DEGs. Among 2316 genes expressed at 0h, 1639 were upregulated, and 677 genes were downregulated. For 6h treatment, 2623 genes were expressed consisting of 2330 upregulated and 293 downregulated genes. In addition, 2390 genes were observed for 24h treatment consisting of 2065 upregulated and 325 downregulated genes. For 48h treatment, we observed a total of 2509 expressed genes consisting of 1749 upregulated and 760 downregulated genes, respectively ( Figures 2C–F ). These comparisons indicated that the differential responses were only associated with the biotype of the rice strains, however, very few changes were observed among treatments of BPH feeding.

In addition to be used for DEGs, RNA-seq dataset is also a good source utilized for the identification of genes involved in metabolic pathways. Therefore, we assessed DEGs involved in different metabolic pathways for two strains using MapMan software and observed the same trend as the number of DEGs with significant differences involved for different metabolic pathways. Most DEGs were observed for signaling, secondary metabolites, and transcription factor in the KW strain. However, these pathways were disturbed significantly in the NIL allowing more DEGs into it ( Figures 3A, B ).

Next, we performed hierarchical clustering of DEGs observed in KW and NIL to find out differences in the expression patterns of these DEGs. Most DEGs recorded different regulation patterns for NIL compared with KW ( Figure 3C ). Meanwhile, we observed many disease-resistant and insect resistance-related proteins and transcription factors have changed expressions, such as xylanase inhibitor protein and LRR receptor. It can also be speculated that plants began to decrease their metabolic activities while increasing their focus on coping with BPH invasion. DEGs associated with phytohormones biosynthesis about ET and SA began to multiply up-regulated expression in the middle stage, and cell wall biosynthesis showed significant changes in the late-stage suggesting that plants also initiated the relevant defense mechanism against BPH.

It was worth noting that down-regulated genes decreased sharply in NIL, probably because strong resistance to insects had developed, so previously down-regulated genes returned to normal expression levels. During the late stage of BPH onset, the DEGs of KW declined sharply while the DEGs of NIL increased significantly, forming a distinct contrast illustrating different expression patterns of the two strains against BPH feeding. These results indicated that differences in the expression pattern of DEGs are due to different strains.

LncRNAs have been implicated in playing a critical role in coding gene expressions (Sebastian-delaCruz et al., 2021). To predict lncRNAs in our transcriptomic dataset, we analyzed the assembled and filtered transcripts procuring approximately 6370 expressed lncRNAs in two strains after BPH feeding, with 384 differentially expressed lncRNAs (DELs) among them ( Figure 4A ). We identify these DELs including some of them related to BPH resistance in rice using PFAM, CPC, CPAT, and CNCI, which were subjected to further analysis. These DELs ranged in length between 255 and 2970 bp with the most abundant length of 300–600 bp ( Supplementary Figure 5 ).

It was expected that lncRNAs and mRNAs could be evenly distributed on chromosomes, but mRNAs will have higher conservation than lncRNAs. The possible reason was that lncRNAs, as a regulatory factor in vivo, regulate different life activities in many ways. Further, significant differences were observed in transcript length between PCGs and lncRNAs ( Figures 4A, B ), suggesting a strong influence of lncRNAs on PCGs transcript length. We also calculated the expression level (FPKM) of the expressed PCGs and lncRNAs, depicting an average expression of 28 for the coding genes and approximately 4 for the non-coding targets ( Figure 4C ).

It has been confirmed that lncRNAs had fewer exons and were shorter in length than PCGs. The length of the expressed lncRNAs in different samples was compared with those of predicted mRNAs. An uneven length of lncRNAs was assessed in the study, however, the length distribution of mRNAs was relatively uniform ( Supplementary Figure 5 ). Moreover, lncRNAs had a smaller number of exons in contrast to PCGs. The number of exons in mRNAs ranged from 1 to 14, whereas most of the lncRNAs contained only one or two exons. The number of mRNAs containing one exon was the largest (> 25000), and then decreased with the increase in the number of exons ( Supplementary Figure 6 ).

Conservation of certain genes is essential for the maintenance of life activities and genetic stability. It is generally believed that mRNA is more conserved than lncRNA. Therefore, we assessed the influencing behavior of DEGs and DELs through hierarchical clustering and found that most DELs have positively regulated their neighboring coding genes ( Figure 4D ). Focusing on the comparison of KW and NIL at 0h, 6h, 2h, and 48h time points, some of the DELs, e.g., MSTRG.65184, MSTRG.44317, MSTRG.17612, etc. have a strong influence over the expression of neighboring coding targets ( Figure 4D ). Other examples of DELs that positively influenced their mediated target genes include MSTRG.16891 targeting Os04g0670400 and MSTRG.22175 mediating Os11g0303400 ( Figures 4E, F ).

Hence we observed that lncRNAs have significant influences over the expression of coding genes, therefore, we wanted to check whether lncRNAs have differential expression changes between treatments of the two strains. Unlike DEGs, DELs have expressed non-significant differences among treatments, except for 48h treatment where the number of expressed lncRNAs was also less compared to other treatments ( Figures 5A, B and Supplementary Table 7 ). These results indicated that besides differential background, lncRNAs have also changed the expression patterns of the coding targets. Comparisons of lncRNAs among different treatments showed upregulation of lncRNAs for all treatments except for 48h treatment where the number of downregulated lncRNAs exceeded upregulated lncRNAs supporting our results in the case of DEGs ( Figures 5C–F ).

A major theme involves is the regulatory role of lncRNAs, which trigger the expression of neighboring PCGs. To predict lncRNAs and coding genes interactions, we used Cytoscape (http://www.cytoscape.org/) as visualization tool and constructed the putative interactive network of lncRNAs targeting their presumed PCGs ( Figure 6 ). Correlation expression was the input data for the interaction networks between lncRNAs and their target coding genes. In the map, we presented four lncRNAs (MSTRG.9425, MSTRG.16026, MSTRG.28482, and MSTRG.22186) that targeted their interactive PCGs influencing their expression further confirmed the involvement of lncRNAs mediating coding genes targets. These analyses showed metabolic and signal transduction pathways mainly affected by DEGs, alternatively by DELs. According to previous reports, these pathways are closely related to the resistance of plants to insects, suggesting that lncRNAs can be involved in the regulation of rice resistance to BPH. These results showed that, compared to KW, NIL plants were more inclined to strengthen their defense response by expressing more genes involved in the defense process.

Apart from other differentiation in key pathways among the two strains, biotic and abiotic stress pathways were also strongly regulated under BPH treatments. NIL has a strong influence over KW when comparing biotic stress regulatory genes. The same regulatory trend was observed for signaling genes, cell wall-related genes, and various TFs involved in our study ( Figure 7 ). However, more secondary metabolites-related genes were perceived in KW compared to NIL, suggesting that pathogen attack has strongly influenced the metabolic pathways concerning the background of the strain ( Supplementary Figures 7A, B ). Moreover, BPH treatments have also impacted photosynthesis-related pathways among the two strains demonstrating differential responses of the strains toward pathogenic attack. Overall, these results demonstrate that the genetic background of rice genotypes has a strong influence over responses to pathogenic attacks and can be used as a strategy to overcome stresses.