Introduction

Rice is an important food crop worldwide, and the development of improved rice varieties is crucial for ensuring global food security. Pest infestation is among the most significant stress factors affecting rice yield. The brown planthopper (BPH) Nilaparvata lugens (Stål), a monophagous insect that feeds on rice, has caused extensive damage to the rice industry in China and many other parts of Asia, leading to a marked reduction in yields1,2.

At present, the use of chemical insecticides is the fastest and most powerful method for managing BPH3. However, chemical control methods can be detrimental to human health and the environment, and have brought about insecticide-resistant BPH biotypes4,5,6. Rice has developed a complex defense system against BPH, and BPH-resistant rice varieties are mainly obtained through conventional sexual hybrid breeding; however, this requires a long breeding time and a large workload in seed selection, and it is difficult to rapidly obtain excellent resistant varieties. Therefore, utilizing insect resistance genes in rice to develop resistant rice varieties represents a promising and sustainable method for controlling rice planthopper infestation. Exploring the genes for broad-spectrum and persistent resistance to BPH is a prerequisite of gene aggregation breeding for the selection of new varieties to control BPH damage7. The effects of different resistance genes on the physiological metabolism of BPH are different, indicating that the mechanism of action is specific to rice variety, and transcriptome analysis provides an important basis for further understanding of its function8,9. Transcriptomic and metabolomic association analysis showed that rice promoted BPH resistance by inducing epigallocatechin and reducing indole acetic acid (IAA)10. In addition, transcriptome sequencing identified 14,358 DEGs and 55 potential BPH stress players related to resistance to BPH11. Furthermore, single cell sequencing showed that mesophyll cells may regulate the expression of vanillin, capsaicin and reactive oxygen species (ROS) production genes, phloem cells may regulate cell wall extension genes, and xylem cells may participate in BPH resistance by controlling the expression of chitin and pectin related genes9. The accumulation of callose in resistant rice prevented BPH infestation, and the exogenous abscisic acid (ABA) enhanced rice resistance to BPH by promoting callose formation12,13.

In recent years, there have been increasing outbreaks of the BPH in China, which is a serious threat to rice yield. Therefore, the discovery of the BPH resistance genes is growing more relevant. Mudgo, a rice germplasm that is resistant to BPH, was first discovered in 196914, followed in 1971 by the identification of two key BPH resistance genes, Bph1 and Bph2. Since then, 40 genes involved in the resistance of rice to BPH have been found in the rice genome, and many of them have been cloned15,16,17. The mechanisms of action differ among the resistance genes. For instance, Bph14, which is located on chromosome 12, encodes a typical CC-NB-LRR protein18; Bph619 and Bph3020, which are both found on chromosome 4, encode proteins containing leucine-repeat domains; Bph3, which was identified in the wild rice cultivar Rathu Heenati, is an OsLecRK1-OsLecRK3 tandem gene cluster that encodes G-type lectin receptor kinases localized to the plasma membrane of cells21. Systemic resistance to BPH was also originally present in rice without the cloned resistance genes, but during its long co-evolution with rice, BPH has overcome the effects of several originally existing resistance genes, the most notable being Bph1 and Bph222,23. This highlights the need to identify additional planthopper resistance genes for rice resistance breeding to improve the persistence and broad spectrum of resistance to BPH.

In our previous study, using a seedbox screening method, we found a rice variety (R26) that was immune to BPH24; however, the mechanism underlying the resistance to this pest remained unclear. In the current study, using RNA sequencing (RNA-seq) of the leaf sheath, we sought to delineate the differences in the transcriptomic profiles of resistant and susceptible rice varieties following feeding by BPH. Last, the function of the BPH resistance related gene (OsREM4.1) was analyzed in knockout transgenic rice. This work will provide a reference for unraveling the molecular mechanism involved in the resistance of rice to insects and provide a theoretical basis for the improved breeding of BPH-resistant varieties.

Results

Overview of the transcriptomic profiles of rice seedlings inoculated with BPH Nymphs

In previous research, using the Standard Seed Box Screening Technique, we were able to identify a new rice variety called R26 that was resistant to BPH (Fig. 1A)24. However, R26 did not contain any known resistance genes, which are identifiable with gene chips if existent. When using gene chip detection, we were unable to locate resistance genes, such as Bph6, Bph9, Bph14, Bph15, Bph18, and Bph26. Since the material (R26) could contain a new BPH-resistance gene, transcriptomics was used to examine it.

Fig. 1
figure 1

Statistical analysis of the expressed genes, transcription factors (TFs), and differentially expressed genes (DEGs) after brown planthopper (BPH) infestation. (A) The phenotype of resistance of TN1 and R26 to BPH. (B) Line charts showing the numbers of expressed genes after BPH infestation in the TN1 and R26 rice varieties. (C) Venn diagram of the genes expressed in the two varieties. (D) The TOP10 transcription factors among the expressed genes. The numbers represent the numbers of genes in the families.

We undertook a transcriptomic analysis of TN1 (sensitive, control) and R26 (resistant, case) rice seedlings 0, 0.5, 1, 3, 5, and 7 days after inoculation with BPH (Fig. 1B). After removing adapters and low-quality reads, a total of 244 Gb of clean data were obtained, with an average of 45.21 million reads per sample; more than 93.7% of the bases were assigned a Q score of 30 (Table S1). In total, 93.44–95.46% of the clean data were mapped to the reference genome. These results demonstrated that the RNA-seq data were of reliable quality. All sequencing data were uploaded to the National Center for Biotechnology Information (NCBI) (Accession ID: PRJNA1003552).

The PCA showed that the biological replicates of the two rice cultivars had good consistency, and the same material from different infection stages (time points) was clustered together, implying that the difference between the cultivars was greater than that between infection stages (Figure S1). The PCA further showed that the gene expression patterns related to resistance to BPH differed between the two rice varieties. The resistance of rice to BPH infection was seen to be a dynamic and complex process involving the coordinated activity of an intricate regulatory network comprising numerous genes. After infestation with BPH, the number of genes expressed in TN1 was lowest on day 1, while in the resistant variety, the number of expressed genes decreased to its lowest level after half a day (Fig. 1B). A total of 24,367 and 23,217 genes were identified as being expressed in TN1 and R26; of these, 21,958 were co-expressed, and 1259 were expressed only in the resistant variety (Fig. 1C). We conducted a transcription factor alignment analysis on the expressed genes and identified a total of 1891 transcription factors belonging to 56 families (Table S2). We also analyzed the transcription factors (TFs). Among these, 1891 genes were identified as TFs, which consist of 56 TF family (Table S2), and the bHLH, NACs, ERF, MYB, C2H2, and WRKY families all contained more than 100 of the expressed genes (Fig. 1D). We analyzed the DEGs between different time points following artificial inoculation with BPH nymphs and found that the number of DEGs in the two rice cultivars first in-creased, peaking after 1 day, and subsequently decreased (Figure S2). Interestingly, the number of DEGs in TN1 was significantly higher than that in the resistant variety. Additionally, compared with day 0, TN1 exhibited the fewest DEGs on day 3, while the number of DEGs in the resistant variety decreased over time from the day of inoculation (Figure S2).

DEGs in R26 and TN1 plants before and after BPH feeding

To capture genes associated with resistance to BPH, we generated Venn diagrams of differentially expressed genes at different time points to identify genes that were induced in the resistant varieties but were differently expressed in the susceptible varieties. In this analysis, genes were identified across five time points of BPH infection. Compared with the control group, 426 and 159 genes were down-regulated and up-regulated in TN1 (Fig. 2A and B) respectively, and 159 and 35 genes were down-regulated and up-regulated in R26 (Fig. 2C and D). At clusters 7 and 8, some of the up-regulated genes in R26 belonged to NB-LRR family of disease-resistance genes, such as BGIOSGA033552 and BGIOSGA034264, and lignin-synthesis-related genes were also highly expressed in R26 after 1 day of infestation with BPH (BGIOSGA000512, BGIOSGA012707, and BGIOSGA036496) (Fig. 2E). In cluster 10, an auxin-responsive protein IAA27 (BGIOSGA035021) was highly expressed (Fig. 2E).

Fig. 2
figure 2

An overview of the transcriptomic analysis. (A) Venn diagram of down-regulated genes in comparison with those in TN1. (B) Venn diagram of up-regulated genes in comparison with those in TN1. (C) Venn diagram of down-regulated genes in comparison with those in R26. (D) Venn diagram of up-regulated genes in comparison with those in R26. (E) Differential gene expression analysis showing up- and down-regulated genes across all sixteen clusters.

Identification of the core genes by using WGCNA

A weighted gene co-expression network analysis (WGCNA) was performed based on the genes displaying differential abundance between TN1 and R26 over the 7 days. Beta (power) = 21 was selected as the optimal soft threshold for constructing a scale-free network (Figure S3A) and displaying network connectivity with different soft thresholds (Figure S3B) to further establish a weighted co-expression network model, classify genes, and divide them into different modules. This resulted in a total of 12 modules, each represented by a different color in Fig. 3A. The turquoise module contained the greatest number of genes (245), and the gray module contained the smallest number (5). A correlation analysis was conducted between different BPH inoculation times and modules, and a correlation heat map was drawn (Fig. 3B). The heat map analysis indicated that the blue module was significantly and positively correlated with CK0d (TN1 control) (0.71, 1e-06), and the green module and yellow modules were significantly correlated with CK3d (0.66, 1e-05) and CK1d (0.66, 1e-05), respectively. The magenta module was significantly correlated with CK7d (0.71, 1e-06). The red module and pink modules were significantly correlated with M0.5d (R26 infestation for 12 h) (0.63, 4e-05 and 0.52, 0.001, respectively). As can be seen in Fig. 3B, the expression patterns of the blue, green, magenta, yellow, and red modules showed the same trends in TN1 and R26, and there was a higher correlation with TN1. The pink module showed directly opposing trends in the TN1 and R26 varieties: it was positively correlated with R26 and negatively correlated with TN1, suggesting that it may be related to resistance to BPH. The trend of gene expression in the pink module was similar to the trends in cluster 1 and cluster 3 in the MaSigpro analysis (Figure S4 and Table S3).

Fig. 3
figure 3

Gene module analysis using a weighted gene co-expression network analysis (WGCNA). (A) Hierarchical clustering diagram of the modules. (B) WGCNA of the genes with dominant expression in the two rice varieties after brown planthopper (BPH) infestation. Red indicates a positive correlation and blue indicates a negative correlation. The darker the color, the stronger the correlation. The numbers in parentheses represent significant P-values; the smaller the value, the greater the significance. (C) KEGG pathway analysis of the unique DEGs in the pink module. (D) Co-expression regulatory network analysis of the pink module.

In order to determine the probable regulatory roles of the DEGs, and thus the possible functions of the target modules, we also carried out a KEGG pathway enrichment analysis. It was discovered that 17 pathways implicated the DEGs (Fig. 3C). Eight of the 17 pathways, such as flavonoid biosynthesis, the MAPK signaling pathway (plant), biosynthesis of other secondary metabolites, signal transduction, plant–pathogen interaction, environmental adaptation, phenylpropanoid biosynthesis, and environmental infor-mation processing, were related to the plant immune response. To gain a deeper comprehension of the resistance of R26 to BPH, a functional protein association network was generated by entering all identified genes from the pink module into the STRING database. Central to the network and associated with lignin synthesis were genes related to PAL and 4CL3 (Fig. 3D). The central genes were functionally annotated via alignment with the rice database (Table S4). We analyzed the expression of these central genes at each dpi, which reached the peak at 12 h in R26(Figer S5). These phenomena indicate that these central genes may be the key genes for resistance to BPH.

In addition, RT-qPCR analysis of nine genes (BGIOSGA016973, BGIOSGA035600, BGIOSGA025904, BGIOSGA027799, BGIOSGA012776, BGIOSGA035957, BGIOSGA008199, BGIOSGA008200, and BGIOSGA016055) (Fig. 4) showed that the expression trend for the RT-qPCR results was consistent with that for the transcriptome data.

Fig. 4
figure 4

Verification of gene expression levels by using RT-qPCR. (AI) The bar chart shows the expression levels of nine genes from TN1 and R26 as determined through RT-qPCR; the heat map shows the expression levels of the same genes as determined through RNA-seq analysis. The bar is the standard error of each group of processing data (n = 3).

Resistance of G01-8 mutants to the BPH

In order to detect whether the remorin protein plays a role in the resistance of rice to BPH, we mutated OsREM4.1 by using CRISPR/cas9 in ZH11 and found that we obtained a mutant with an additional base pair at the target position. This mutant plant was named g01-8 (Fig. 5A). To verify the resistance of remorin to BPH, a seedling resistance analysis was performed using wild-type ZH11 plants as the control and TN1 (susceptible) and Rathu Heenati (RHT, resistant) plants as additional controls. After 7 d of BPH infes-tation, the TN1 plants died and the RHT plants survived well, while half of the ZH11 plants survived and the G01-8 plants died, thus demonstrating the effect of the OsREM4.1 gene on resistance to BPH (Fig. 5B, C). The callose content was tested as well, the results showing that the callose content was also lowest in G01-8 plants (Fig. 5D).

Fig. 5
figure 5

Verification of remorin-4.1-mediated resistance to BPH. (A) Point detection of the mutation position. (B) Photo showing the difference in ZH11 and G01-8 rice plants’ BPH resistance. (C) BPH-resistance scores of ZH11 and G01-8 plants (negatively related to BPH resistance) obtained from observations 7 days after BPH infestation (n = 3 independent experiments, each with 10 rice plants). (D) Callose content in leaf sheath of different rice plants. The bar is the standard error of each group of processing data (n = 3).

Schematic diagram of key potential mechanisms for regulating resistance to BPH

We created a schematic diagram to summarize the key mechanisms of the resistance of R26 rice to BPH (Fig. 6). When BPHs fed on rice seedlings, many biological processes were significantly affected, such as plant signaling pathways, the MAPK signaling pathway (plant), biosynthesis of other secondary metabolites, and plant–pathogen interactions. Similar to previous research results, we suggest that when the BPH feed on rice plants, they first activated plant signaling pathways, such as SA, JA, and MAPK immune cascade reactions25,26. These signaling pathways cause changes in plants’ secondary metabolic pathways, leading to the strengthening of rice cell walls. In agreement with previous findings, more secondary-metabolite-related genes were detected in R26, indicating that pathogen attacks strongly affect metabolic pathways27. BPH feeding can lead to calcium ion flux, trigger H2O2 accumulation, and cause protein blockage, as well as callose deposition on the sieve28,29. However, in R26, in addition to initiating H2O2 accumulation and causing cell death, increased remorin expression may also lead to decreased intercellular permeability, resulting in callose accumulation and increased resistance to the BPH.

Fig. 6
figure 6

Potential mechanisms of the resistance of R26 to BPH. When BPH nymphs feed on rice seedlings, resistant rice regulates a variety of biological pathways for defense, including plant hormone signal transduction, the MAPK signaling pathway, and biosynthesis of other secondary metabolites. Proteins labeled in red are genes affected and upregulated by BPH infestation.

Discussion

BPH is one of the most prevalent pests endangering rice production. Cultivating rice varieties that harbor insect-resistance genes is considered the most economical, effective, environmentally friendly, and sustainable method for controlling BPH10. Screening resistance genes from different germplasm resources and identifying rice resistance genes against BPH are the keys to breeding insect-resistant rice varieties. Under normal circumstances, differences in gene expression are related to breed inheritance. When BPH infects rice, different resistance genes would cause different resistance reactions. We previously identified a rice variety with persistent resistance to BPH that may contain a novel resistance gene to the BPH by genomic association analysis (GWAS)24. In this study, for the first time, RNA expression patterns in the leaf sheath of insect-sensitive and insect-resistant rice seedlings during BPH infestation were investigated using a time-course analysis. To obtain more accurate expression data related to resistance to BPH in rice, we only selected the leaf sheaths of rice seedlings, as this is where BPH feeding is concentrated. Our data showed that, after infestation with BPH, the number of genes expressed in TN1 was higher than that expressed in the resistant variety (Fig. 1B), and a greater number of genes were specifically expressed in TN1 than in R26 (Fig. 1C). After 1 day of BPH feeding, the genetic changes were the largest, which was different from the results of a previous study30, possibly because the materials contained different resistance genes. These results suggest that when rice was not exposed to BPH, there were significant differences in the number of genes expressed between resistant and susceptible varieties, which are determined by rice genotype, and after BPH infection, the DEGs between resistant and susceptible varieties may be the reason why R26 was resistant to BPH. TFs play a regulatory role in the growth and development, secondary metabolism and stress resistance of plants. In previous research, WRKY, MYB and bZIP were significantly upregulated in resistant plant in the middle stage to increase resistance against BPH11,27, in our study, we also identified WRKYs, bZIPs and MYBs, and half of these TFs were downregulated compared with susceptible varieties, which may not be the main resistance mode in R26.

Genes are a determining factor of the phenotype, and transcriptome sequencing could offer valuable insights into the functionality of rice resistance genes to BPH infestation. Previous studies have also utilized transcriptome analysis to elucidate the mechanisms underlying resistance in different rice varieties against BPH. As BPH primarily feed on rice leaf sheaths and cause significant damage to the crop at seedling stage, this study employed transcriptomic analysis on the leaf sheaths of “TN1” and “R26” with BPH infestation to unveil their resistance mechanisms. Basic resistance exists in both susceptible and resistant rice varieties. Feeding by herbivores can promote the expression of the PLD gene and the release of volatiles from green leaves, while resistance to BPH can be restored through the application of volatiles to green leaves31,32. In addition, BPH infestation was also reported to activate MAPK-, ethylene-, and SA-related signaling pathways18,25. Here, the DEGs between the TN1 and R26 varieties were mainly associated with “phenylalanine, tyrosine, and tryptophan biosynthesis”, “biosynthesis of other secondary metabolites”, and “MAPK signaling pathway (plant)” in the modules assessed in the WGCNA. In addition to inducing the production of metabolites for resistance, pests also induce the production of defensive proteins, which are usually enzymes, such as phenylalanine ammonialyase (PAL), peroxidase (PO), and polyphenol oxidase (PPO). PAL is an important enzyme that acts as a link between plant birth metabolism and secondary metabolism, and plays an important role in plants. It can regulate the resistance to BPH by regulating the biosynthesis and accumulation of SA and lignin33. After BPHs fed on rice plants, the PAL activity in the resistant rice variety was significantly higher than that in susceptible rice34, which was consistent with our experimental results of increases in the gene expression of defense-related secondary metabolites such as PAL resistant rice, and the inhibition of BPH infestation by the high expression of protease inhibitors. In addition, the 4CL was the main downstream target of the PAL pathway, catalyzing the final step to the formation of various cinnamate CoA sulphoids and controlling the pathway into the metabolism of different phenylpropanoid compounds. 4CL used cinnamic acid and its hydroxyl or methoxy derivatives such as 4-coumaric acid, ferulic acid, erucic acid as substrate to produce the corresponding coenzyme A ester. 4CL was also one of the rate-limiting enzymes in lignin synthesis, which it regulated. In this study, through WGCNA analysis, we found that PAL, 4CL3, and 4CL5 were the hub genes of R26 resistance to BPH infection, which may be caused by the BPH infection resulting in increased vanillin synthesis in resistant rice. plants containing vanillin have strong insecticidal and repellent activities35,36,37, making R26 resistant to BPH. Additionally, we found that after 12 h of BPH infestation, the expression of these genes in the resistant rice variety R26 was significantly increased, implying that they were involved in the process of resistance to the pest. In summary, the RNA-seq results indicated that the immune differences between the TN1 and R26 varieties stimulated by BPH infection, such as those involving the MAPK, lignin synthesis, and metabolism-related pathways, are variety-specific. The responses of resistance genes are faster in resistant rice varieties than in susceptible ones, and this helps protect plants against BPH attacks and keeps them alive longer.

Plant hormones play important roles in plant resistance to insects and pathogens. Studies have shown that auxin is crucial to plant growth and development, as well as in plant–pathogen interactions38. Rice containing Bph14 showed an increase in the expression levels of salicylic acid (SA) related genes (EDS1, PAD, and ICS1) after feeding by BPH, while the expression levels of jasmonic acid (JA) and ethylene signaling pathway genes decreased18. After feeding on materials containing the Bph6 gene by BPH, genes related to salicylic acid (SA), jasmonic acid (JA), and gibberellin (CKs) were upregulated, and the response speed was faster than that of susceptible varieties19. In our study, we also found that in resistant varieties, hormones quickly respond to feeding by the BPH, while in susceptible varieties, the response speed is slower. In this study, we found that plant hormone signal transduction pathways were enriched in the brown module and identified IAA1 as the hub gene. IAA1 (BGIOSGA002262) expression levels were lowest in the resistant variety 1 day after exposure to BPH; this result was similar to previous single-cell sequencing results9. Moreover, it was found that auxin can inhibit PR1 expression, IAA1 may be related to insect resistance in resistant varieties, and the antagonism between IAA1 and SA may represent a balance between plant growth and defense39,40. Studies have shown that ABA is a negative regulator associated with plant disease resistance. However, ABA also has a positive effect on resistance to plant diseases. ABA enhances plant resistance to pathogens through mechanisms that include callose deposition and the thickening of callose in the cell wall to form the corpus callosum41. The remorin protein is the main component of the plasma membrane and intercellular desmosis, and participates in the interaction of plant hormones to improve plant disease resistance. It can regulate its temporal and spatial expression in the plasma membrane through phosphorylation, and it can induce callose to regulate cytodesmosis closure, which can effectively limit the intercellular movement and proliferation of PVX virus42. By knocking out the rice OsREM4.1 gene, we found that it was more susceptible to BPH than the wild type was (Fig. 6). it is possible that OsREM4.1 mutation causes callus failure to accumulate and reduces rice resistance to BPH. The results indicate that OsREM4.1 plays a role in rice resistance to BPH by regulating callose accumulation, and may provide attachment medium for other anti-insect compounds; the molecular mechanism needs further study. The results of this study will be helpful for the early screening and molecular marker-assisted breeding of rice varieties resistant to BPH.

Methods

Plant materials and insect populations

The rice varieties used in this study were the susceptible variety TN1 and the insect-resistant variety R26 identified in our previous study24. Seeds were provided by the Rice Research Institute, Fujian Academy of Agricultural Sciences (Fuzhou, China). The rice seeds were soaked at 26 °C for 48 h, exposed, and sown in seedling trays (40 × 30 × 8 cm) containing paddy soil. The seedling trays were evenly divided into 10 rows with 12 seedlings per row, two varieties in alternate rows, and three biological replicates per variety. The seedlings were grown in a greenhouse for 7 days (28℃, 14 h/10 h light/dark photoperiod). The BPH-resistance genes were detected at the Wuhan Greenfafa Institute of Novel Genechip R&D Co., Ltd (Greenfafa, Wuhan, China) using an Illumina SLYm1R SNP chip containing 31,753 SNPs, and the chromosomal locations of each probe on the SLYm1r SNP chip were determined by using the Nipponbare MSU7.0 as a reference genome.

Plant infestation with BPH and sample collection

BPH specimens were provided by the Rice Research Institute, Fujian Academy of Agricultural Sciences (Fuzhou, China). BPH feeding was performed as previously described24. The BPH specimens were reared on TN1 in the greenhouse at a temperature of 26 ± 2 °C and a relative humidity of 65 ± 5% with a 14 h/10 h light/dark photoperiod (artificial lighting). Three trays containing 10-day-old seedlings were covered with insect-proof nets, and the seedlings of both varieties were artificially inoculated, with 8–10 nymphs of 1–2 instar BPHs per seedling. Leaf sheaths fed by the BPHs that showed consistent growth were collected at 0, 0.5, 1, 3, 5, and 7 days after infestation. A total of 36 samples were collected, including three biological replicates. The stems were immediately placed in liquid nitrogen and stored at − 80 °C for subsequent RNA extraction.

RNA extraction and transcriptome sequencing

Total RNA was isolated by using Trizol reagent (Invitrogen Life Technologies), after which the RNA concentration, quality, and integrity were determined by a NanoDrop spectrophotometer (Thermo Scientific). Three micrograms of RNA was used as the input material for RNA sample preparation. The prepared RNA libraries were sequenced on the NovaSeq 6000 platform (Illumina) by Shanghai Personal Biotechnology Co., Ltd. Adapters and low-quality reads were filtered out by using Cutadapt (v1.15), yielding high-quality, clean data for further analysis. The clean data were aligned with the rice reference genome (Oryza sativa indica, ASM465v1; http://plants.ensembl.org/Oryza_indica/Info/Index) using Hisat2 (v2.0.5)43. The original expression of the genes was compared with the read count for each gene by using HTSeq (0.9.1). Gene expression was then normalized based on fragments per kilobase of exon per million fragments mapped (FPKM) reads. Genes with FPKM > 1 were considered to be expressed. Differentially expressed genes (DEGs) were identified using DESeq (version 1.30.0) with |log2FoldChange|>1 and a P-value of < 0.05 as the threshold. MaSigPro (version 1.64.0)44 was employed for time-course analysis. The gene expression matrix for all of the samples was input into MaSigPro with the susceptible rice variety (TN1) as the control check. Then, the design matrix was constructed based on the two rice varieties and six developmental time points. The degree of polynomial regression was set to 9, and the forward elimination algorithm was used for stepwise regression, with R2 > 0.9 as the threshold; the significantly fitting genes were differentially expressed over time between the two varieties. Hierarchical clustering (hclust) was performed based on the correlation coefficient distance.

All of the genes were mapped to terms in the Gene Ontology (GO) database, and the number of differentially enriched genes for each term was calculated. The topGO package was used to analyze the GO term enrichment for the DEGs, and the hypergeometric distribution method was used to calculate the P-value (enrichment P-value: <0.05). The main biological functions of the differentially abundant genes were determined based on the identified significantly enriched GO terms associated with the DEGs. The TBtools (version 1.120)45 software was used for a KEGG pathway enrichment analysis (P < 0.05) of the DEGs.

Weighted gene co-expression network analysis

Gene co-expression networks were analyzed by using the WGCNA (version 1.72) package in R46. The co-expression modules were obtained by the automatic network construction and module detection function (blockwiseModules) with the following default parameters: soft threshold power = 10, TOMtype = signed, mergeCutHeight = 0.25, and minModuleSize = 30. Protein interactions were predicted using the STRING database. The hub genes were selected and visualized with Cytoscape (version 3.9.1).

Construction of the REM4.1 mutant and BPH resistance evaluation

The OsREM4.1 mutant was created using CRISPR/Cas9 editing technology, and the target sequence for the mutation was 5’-CGTCGTCGTCGTACCACCGGCGG-3’. Additionally, ZH11 rice was utilized to create the mutant. R: TTGTTGACCTTGGCGACCTC and F: CAAACGGCGGCTAGTGGTAG were the mutant detection primers. Additionally, the previously established approach to evaluating the phenotypic resistance of rice that was variably resistant to BPH was used24. The callose content was determined with plant callose ELISA kit.

Quantitative real-time PCR analysis and data statistical analysis

Total RNA was extracted as described for RNA sequencing (Sect. 2.3). The isolated RNA was reverse-transcribed to cDNA using the RT First Strand cDNA Synthesis Kit (Servicebio, G3330). qPCR was conducted on a CFX Connect (Bio-Rad). Each run had three biological replicates and three technical replicates. Primers were designed using Primer Premier 5 (Table S5). The relative expression levels of each gene were calculated by using the 2-ΔΔCt method, with glyceraldehyde-3-phosphate dehydrogenase (GAPDH) as the reference gene.

The statistical calculations and histogram drawing of experimental data were performed using Graphad Prism 7. Statistical analyses were performed using one-way ANOVA or Student’s t tests, and p-values < 0.05 were considered to indicate statistical significance.