Combined miRNA and mRNA sequencing reveals the defensive strategies of resistant YHY15 rice against differentially virulent brown planthoppers

Introduction The brown planthopper (BPH) poses a significant threat to rice production in Asia. The use of resistant rice varieties has been effective in managing this pest. However, the adaptability of BPH to resistant rice varieties has led to the emergence of virulent populations, such as biotype Y BPH. YHY15 rice, which carries the BPH resistance gene Bph15, exhibits notable resistance to biotype 1 BPH but is susceptible to biotype Y BPH. Limited information exists regarding how resistant rice plants defend against BPH populations with varying levels of virulence. Methods In this study, we integrated miRNA and mRNA expression profiling analyses to study the differential responses of YHY15 rice to both avirulent (biotype 1) and virulent (biotype Y) BPH. Results YHY15 rice demonstrated a rapid response to biotype Y BPH infestation, with significant transcriptional changes occurring within 6 hours. The biotype Y-responsive genes were notably enriched in photosynthetic processes. Accordingly, biotype Y BPH infestation induced more intense transcriptional responses, affecting miRNA expression, defenserelated metabolic pathways, phytohormone signaling, and multiple transcription factors. Additionally, callose deposition was enhanced in biotype Y BPH-infested rice seedlings. Discussion These findings provide comprehensive insights into the defense mechanisms of resistant rice plants against virulent BPH, and may potentially guide the development of insect-resistant rice varieties.


Introduction
Rice (Oryza sativa L.) was domesticated approximately 10,000 years ago in the lower Yangtze Valley in China.From there, it spreads across Asia, Africa, Europe, and the Americas, and now serves as a staple crop for more than half of the world's population (Cheng et al., 2013;Gutaker et al., 2020).However, rice production suffers from numerous pests and pathogens.Among them, the brown planthopper (BPH; Nilaparvata lugens Stål) is considered extremely destructive (Bottrell and Schoenly, 2012;Cheng et al., 2013;Jing et al., 2017).BPH feeds on phloem sap and causes dwarfing, wilting, browning, drying, and ultimately death in severe cases (Wang et al., 2008).In addition, BPH serves as a vector for viral diseases, resulting in significant yield shortfalls and economic losses (Cheng et al., 2013).
Over the course of an evolutionary arms race between these two species, rice has developed sophisticated defensive mechanisms against BPH (Cheng et al., 2013), including both basic defenses and resistance (R) gene-mediated defenses.To date, approximately 40 major BPH resistance genes have been identified in cultivated and wild rice species (Wang et al., 2023).Molecular cloning and functional characterization of BPH resistance genes have clarified the molecular mechanisms of rice resistance to BPH (Zheng et al., 2021).One mechanism involves the occlusion of sieve tubes with callose, which is a common plant defense against sap-sucking insects.This mechanism is the most effective in rice varieties carrying BPH resistance genes such as B5 (carrying Bph14 and Bph15), RI35 (carrying Bph14), and YHY15 (carrying Bph15) (Hao et al., 2008).In susceptible varieties such as Taichung Native1 (TN1), b-1,3-glucanases (which are only weakly induced in resistant plants) decompose the deposited callose and thereby facilitate continuous feeding by BPH (Hao et al., 2008).Bph14 encodes a typical CC-NB-LRR protein which interacts with the transcription factors (TFs) OsWRKY46 and OsWRKY72 to increase the expression of the receptor-like cytoplasmic kinase gene RLCK281 and the callose synthase gene by binding to their promoters (Du et al., 2009;Hu et al., 2017).Bph15 is located on the short arm of chromosome 4 and is composed of a gene cluster encoding three plant lectin receptor-like kinase proteins (LecRLKs) (Yang et al., 2004;Lv et al., 2014;Xiao et al., 2016).However, the precise mechanism by which Bph15 confers resistance to BPH is unknown.
Phytohormones are thought to play pivotal roles in the interaction of rice plants and BPH herbivory, including salicylic acid (SA), ethylene (ET), jasmonic acid (JA), cytokinin (CK), brassinosteroid (BR), and abscisic acid (ABA).The SA pathway contributes to the immune response against piercing-sucking insects and is involved in Bph6-, Bph9-, Bph14-, and Bph29mediated resistance in rice (Du et al., 2009;Wang et al., 2015;Zhao et al., 2016;Hu et al., 2017;Guo et al., 2018).Antagonistically, ET negatively regulates BPH resistance in rice (Lu et al., 2011(Lu et al., , 2014;;Ma et al., 2020).However, the role of JA in rice resistance to BPH remains controversial as the silencing of different genes related to JA biosynthesis and signaling results in diverse impacts on BPH resistance.For example, silencing 9-lipoxygenase (OsLOX9/OsHI-LOX, a JA biosynthesis-related gene) enhances BPH resistance, suggesting that JA negatively regulates BPH resistance (Zhou et al., 2009).On the other hand, silencing coronatine insensitive 1 (OsCOI1, a JA receptor gene) has no effect on BPH resistance (Ye et al., 2012).According to studies in allene oxide cyclase (AOC, a JA biosynthesis-related enzyme)-and Myelocytomatosis protein 2 (MYC2, a bHLH TF in the JA pathway)-knockout mutant rice, JA appears to be a positive regulator of BPH resistance (Xu et al., 2021).The exogenous application of JA similarly suggests this (Guo et al., 2018).Recent research also suggests that CK may positively regulate BPH resistance in a JA-dependent manner (Zhang et al., 2022).Conversely, BR promotes BPH susceptibility by modulating SA and JA signaling (Pan et al., 2018).ABA enhances BPH resistance by promoting callose formation (Dinh et al., 2013;Liu et al., 2017) and synergizes with JA to stimulate the expression of TFs in BPH-infested rice (Li et al., 2022).Furthermore, a coordinated CK, SA, and JA signaling network has been found to be activated in Bph6-near isogenic lines (NILs) (Guo et al., 2018).Taken together, these findings suggest that the various phytohormones play diverse roles in the BPH defense response, and that there is complex crosstalk between them.
BPH biotypes with increased virulence have emerged in response to pressures imposed by these defense mechanisms, which are capable of overcoming resistance conferred by major resistance genes.Biotype 1 BPH, which exhibits low virulence on resistant rice varieties, is widely distributed across southeast Asia and primarily parasitizes susceptible varieties such as TN1 (Alam & Cohen, 1998).Rearing biotype 1 BPH on resistant rice variety YHY15 (carrying Bph15) for several years resulted in the development of highly-virulent biotype Y BPH, which are able to overcome resistance conferred by Bph15 (Jing et al., 2011).Similarly, after force-feeding 40 generations of local BPH using resistant IR56 rice (carrying Bph3), the resulting BPH population (IR56-BPH) was able to overcome Bph3-conferred resistance (Zheng et al., 2016).In addition, the Mudgo BPH population has been reported to cause substantial damage to Mudgo rice plants (carrying Bph1) (Ji et al., 2013;Wan et al., 2019).In response to selective pressures imposed by these resistance genes, BPH populations accumulate adaptations over generations which eventually allow them to overcome such resistance mechanisms.In effect, this process constitutes a loss of resistance in formerlyresistant rice varieties against them.However, the effectiveness of resistance strategies used by different rice varieties against BPH populations with varying levels of virulence requires further investigation.Such combinations of avirulent/virulent BPH and resistant rice provide ideal models for studying resistance adaptation mechanisms.
Just as plants have evolved intricate defense mechanisms to protect themselves against herbivorous insects, insects have in turn developed strategies to overcome plant defenses (Liu Q. et al., 2021).Interactions between plants and insects involve an array of molecular, biochemical, and physiological processes occurring at multiple levels.Multi-omics analyses integrate genomics, transcriptomics, proteomics, and metabolomics, as well as other "omics" approaches, in order to clarify the intricate signaling pathways, molecular responses, and biochemical processes involved in the dynamic interplay between plants and insects (Wang et al., 2020;Shi et al., 2023).Transcriptional profiling has aided our understanding of the defense mechanisms utilized by rice against BPH.For example, research suggests that BPH infestation results in the upregulation of genes involved in signaling, oxidative stress, pathogen-related response, and macromolecule degradation, as well as the downregulation of genes associated with flavonoid biosynthesis, photosynthesis, and cell growth (Zhang et al., 2004;Yuan et al., 2005;Wang et al., 2008).In addition, microarray analyses of BPH-infested Rathu Heenati and TN1 rice underscore the importance of TFs and phytohormones in the defense response (Wang et al., 2012;Li et al., 2017).
MicroRNAs (miRNAs) are small (approximately 21 nt in length) regulatory RNAs produced through endonucleolytic processing of hairpin precursors (Axtell and Meyers, 2018).miRNAs regulate gene expression by binding to complementary sequences in mRNA molecules, resulting in degradation and/or translational inhibition (Bartel, 2009).In plants, miRNAs are involved in various processes such as phytohormone signaling; abiotic and biotic stress response (Zhang et al., 2013(Zhang et al., , 2016;;Li et al., 2017;Natarajan et al., 2018;Salvador-Guirao et al., 2018;Zhang et al., 2018;Liu X. et al., 2021); and leaf, flower, shoot, root, and vascular tissue development (Marin et al., 2010;Peng et al., 2014;Mangrauthia et al., 2017).However, only a few miRNAs have been found to play roles in the insect-plant interaction (Jing et al., 2023).For example, in an investigation of resistant and susceptible rice varieties, the BPH-responsive miRNAs miR156 and miR396 were found to negatively regulate BPH resistance by regulating JA and flavonoid biosynthesis, respectively (Wu et al., 2017;Ge et al., 2018;Dai et al., 2019).Studies of BPH-responsive mRNA and miRNA transcriptomes have uncovered certain universal responses of rice to BPH infestation (Li et al., 2017;Wu et al., 2017;Nanda et al., 2020).Additionally, integrated expression profiling of miRNAs and target genes associated with the BPH-rice interaction has been conducted (Tan et al., 2020).Such combined miRNA and mRNA analyses will be key to unraveling the transcriptional responses of rice to BPH infestation.
In this study, we employed high-throughput sequencing to analyze the mRNA and miRNA expression profiles of YHY15 rice seedlings infested with either biotype 1 BPH (avirulent population) or biotype Y BPH (virulent population).Furthermore, we combined sequence analysis and physiological assays to reveal the underlying resistance mechanisms of rice against BPH.The findings presented in this work will provide a valuable resource for further genomewide investigations of BPH-responsive genes, as well as studies of Bph15-mediated resistance.Moreover, these findings improve our understanding of the intricate interactions between rice and BPH, and may be used in the development of effective BPH management strategies.

Plant and insect materials
In this study, we used two rice varieties (Taichung Native1 [TN1] and YHY15) and two BPH populations (biotype 1 and biotype Y).TN1 is a susceptible rice cultivar, while YHY15 is a recombinant inbred line (RIL) derived from the RI93 × TN1 F2 population carrying resistance gene Bph15 (Yang et al., 2004).Biotype 1 BPH (avirulent) originated from Wuhan University (China) and are reared on TN1.Biotype Y BPH (virulent) were developed by rearing biotype 1 BPH on YHY15 plants beginning in January 2007 (Jing et al., 2011).All rice plants were grown from seeds sown in sponges (6 cm diameter, 2 cm height), with 8 rice plants per cup.Rice plants were reared in a controlled-environment incubator maintained at 30 ± 2 °C during daytime hours (16 h, 06:00-22:00) and 28 ± 2 °C during nighttime hours (8 h, 22:00-06:00).Rice plants were grown for approximately 2-5 weeks following sowing, depending on experimental needs.The BPH populations were reared at Xinyang Normal University (China) under the following conditions: 26 ± 1 °C, 16 h light/8 h dark cycle.Third instar BPH nymphs were used for the infestation experiments.

Evaluation of BPH resistance in rice
At the two-leaf (2-week-old) stage, YHY15 rice plants were infested with third instar biotype 1 or biotype Y BPH nymphs at a rate of 15 nymphs per seedling.The growth status of each plant was photographically recorded daily until all biotype Y-infested se edlings died.Each expe rimen t consiste d of three biological replicates.

Measurement of BPH weight gain and honeydew excretion
BPH weight gain and honeydew excretion were measured as described previously (Shi et al., 2021).Briefly, newly-emerged adult female BPH were weighed using an electronic balance (Mettler Toledo, MS105DU, Switzerland) and subsequently separated into pre-weighed parafilm sachets (2 × 2.5 cm) fixed to the leaf sheaths of 4-week-old rice plants.After 48 h, the insects were carefully removed from the sachets, and both the insect and the honeydew in each sachet were separately weighed.Weight gain was calculated by comparing each insect's weight before and after feeding, and the weight gain ratio was calculated by dividing the weight gain by the initial weight.Both the weight gain and honeydew excretion assays were conducted using at least 37 replicates.

BPH infestation and sample collection
Three-week-old YHY15 rice seedlings were infested with thirdinstar biotype 1 (RT) and biotype Y (RY) BPH at a rate of 15 nymphs per seedling.Each treatment consisted of three biological replicates, with six seedlings per replicate.For RNA-seq and miRNA-seq analysis, leaf sheaths were collected from noninfested controls (0 h), and during early (6 h) and late (48 h) infestation.According to their time of collection, samples of noninfested rice plants were labeled 'R0', samples of biotype 1-infested rice plants were labeled as 'RT6' or 'RT48', and samples of biotype Y-infested rice plants were labeled as 'RY6' or 'RY48'.Each sampled leaf sheath blade was excised, frozen in liquid nitrogen, and stored at -80 °C for further use.
2.5 mRNA transcriptome sequencing and analysis 2.5.1 RNA extraction, quantification, and qualification Total RNA was isolated from rice samples using RNA Trizol reagent (Life Technologies, NY, USA), following the manufacturer's instructions.RNA degradation and contamination were evaluated using 1% agarose gels.RNA purity was quantified using a NanoPhotometer spectrophotometer (Implen, CA, USA).RNA integrity was assessed using an RNA Nano 6000 Assay Kit for the Bioanalyzer 2100 system (Agilent Technologies, CA, USA).

Library construction, quality control, and sequencing
RNA libraries were generated using 1 µg of RNA per sample and the NEBNext Ultra RNA Library Prep Kit for Illumina.Index codes were incorporated to distinguish between samples.Briefly, mRNA was purified using poly-T oligo-attached magnetic beads, followed by fragmentation using divalent cations.First-strand cDNA synthesis was carried out using random hexamer primers and M-MuLV Reverse Transcriptase.The second strand was synthesized using DNA Polymerase I and RNase H. Overhangs were blunted using exonucleases or polymerases and NEBNext adaptors were ligated following adenylation.The library fragments (250~300 bp) were purified using an AMPure XP system.The USER Enzyme was applied to size-selected, adaptorligated cDNA prior to PCR.PCR was carried out using Phusion High-Fidelity DNA polymerase, Universal PCR primers, and Index (X) Primer.The PCR products were purified using an AMPure XP system, and library quality was evaluated using an Agilent Bioanalyzer 2100 system.Index-coded samples were clustered with a cBot Cluster Generation System using a TruSeq PE Cluster Kit v3-cBot-HS (Illumina), according to the manufacturer's instructions.After clustering, the libraries were sequenced on an Illumina NovaSeq platform, generating 150 bp paired-end reads.

Data analysis
Raw data (fastq format) were first processed using in-house perl scripts.In this step, low-quality reads, adapter sequences, and poly-Ns were removed.Subsequently, the Q20, Q30, and GC content of the clean reads were calculated.All downstream analyses utilized only clean, high-quality data.Next, the reference genome and gene model annotation files were downloaded.An index of the reference genome was constructed, and clean paired-end reads were aligned to the reference genome, using Hisat2 (v2.0.5).The mapped reads from each sample were assembled with StringTie (v1.3.3b) using a reference-based approach.featureCounts (v1.5.0-p3) was used to count the number of reads mapped to each gene.The Transcripts Per Kilobase Million (TPM) of each gene was quantified based on the gene length and the number of reads mapped to the gene.
Differentially expressed genes (DEGs) were identified using the R (v1.16.1) package DESeq2, with three biological replicates per treatment, according to the following criteria: P-value < 0.05, false discovery rate (FDR) < 5, and absolute value of log 2 fold change (FC) ≥ 1.The DEGs underwent additional screening through soft clustering using the Mfuzz package, employing a fuzzy c-means algorithm, as previously reported (Kumar and Futschik, 2007).DEGs exhibiting similar expression patterns were categorized into 20 clusters, and the genes within these clusters were subjected to Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses.GO annotations were downloaded from NCBI (http://www.ncbi.nlm.nih.gov/) and GO (http:// www.geneontology.org/).The KEGG database was used to identify BPH-responsive pathways.Fisher's exact tests were applied to identify significant GO and KEGG categories according to the absolute values of P < 0.05 and FDR < 0.05.

RNA extraction, quantification, and qualification
RNA extraction, quantification, and qualification were conducted as described in section 2.5.1.

Library construction, quality control, and sequencing
Briefly, 3' and 5' adaptors were ligated to the 3' and 5' ends of small RNAs, respectively.Next, first strand cDNA was synthesized after hybridization with the reverse transcription primer.The double-stranded cDNA library was generated through PCR enrichment.After purification and size selection, libraries with 18~40 bp insertions were selected for Illumina sequencing with SE50.The library was quantified with Qubit and real-time PCR and the library size distribution was evaluated with a Bioanalyzer.Quantified libraries were pooled and sequenced on an Illumina platform, according to the effective library concentration and amount of data required.

Data analysis
Raw data (fastq format) were first processed using in-house perl and python scripts.In this step, low-quality reads; reads containing poly-Ns, 5' adapter sequences, or poly-As/Ts/Gs/Cs; and reads missing 3' adapter sequences or insert tags were removed.Subsequently, the Q20, Q30, and GC content of the clean reads were calculated.All downstream analyses utilized only clean, highquality data.Small RNA tags were mapped to the reference sequence using Bowtie (Langmead et al., 2009), either without mismatches or with only one mismatch, to analyze their expression and distribution.Mapped small RNA tags were used to identify known miRNAs, with miRBase (v22.0)used as a reference.miRNAs were identified using a modified version of mirdeep2 (Friedlander et al., 2012) and srna-tools-cli was used to draw the secondary structures.Custom scripts were used to obtain miRNA counts as well as to determine base bias at the first position of identified miRNAs of a certain length and at each position of all identified miRNAs.To remove tags originating from proteincoding genes, repeat sequences, rRNA, tRNA, snRNA, and snoRNA, small RNA tags were mapped to RepeatMasker, the Rfam database, or species-specific data.DEGs were identified using DESeq2.P-values were adjusted using the Benjamini & Hochberg method.A corrected P-value of 0.05 was selected as the threshold for determining significantly differential expression.

Analysis of transcriptional signatures of phytohormone responses
To identify the transcriptional signatures of BPH-responsive phytohormone responses, Hormonometer was used to compare gene expression in rice with gene expression in phytohormonetreated Arabidopsis thaliana (Volodarsky et al., 2009;Liu et al., 2016).Specifically, only orthologous genes detected in the RNA-seq analysis which were related to the Arabidopsis thaliana probe set identifiers were selected for further analyses.

qRT-PCR validation of DEGs
First-strand cDNA was synthesized using a PrimeScript RT Reagent Kit with gDNA Eraser (Takara, Japan), according to the manufacturer's instructions.qRT-PCR assays of candidate genes were conducted using a PrimeScript ™ RT reagent Kit with gDNA Eraser (Perfect Real Time) (Takara, RR047A, China).qRT-PCR was carried out on a CFX96 Real-Time System (Bio-Rad, CA, USA) according to the following protocol: 95 °C for 5 min, followed by 40 cycles at 95 °C for 5 s, 60 °C for 30 s, and 72 °C for 30 s. Relative gene expression was calculated with the 2−11Ct method, using PP2A as the reference gene.All primer sequences are listed in Supplementary Table 1.

Callose staining and evaluation
Callose staining and evaluation were performed as described in a previous study (Hu et al., 2017).Briefly, fresh sheaths collected from two-leaf stage rice seedlings infested with either biotype 1 or Y BPH for 48 h were fixed in an ethanol:acetic acid (3:1 v/v) solution for 5 h.The fixative was changed frequently to ensure thorough fixing and clearing.The samples were then rehydrated successively in 70% ethanol for 2 h, in 50% ethanol for 2 h, and in water overnight.After rinsing three times with water, the samples were treated with 10% NaOH for 1 h to make the tissues transparent.After rinsing four times with water, the samples were incubated in 150 mM K 2 HPO 4 (pH 9.5) containing 0.01% aniline blue for 4 h.Finally, the samples were mounted on a slide and callose deposits were observed with a positive fluorescence microscope (Nikon, Eclipse 80i, Japan) under the UV channel.

Statistical analysis
Statistical analyses were conducted using R (v4.0.4) and SPSS (v22.0)(IBM SPSS, Somers, NY, USA).Two-sided Student's t-tests were used to determine statistically significant differences between groups.All bioinformatics analyses were conducted using R packages.

YHY15 exhibits differential resistance to biotype 1 and biotype Y BPH
YHY15 rice, which contains the BPH resistance gene Bph15, exhibits robust resistance against avirulent BPH biotypes (i.e., biotype 1) (Yang et al., 2004).In this study, YHY15 seedlings were subjected to infestation by either biotype 1 or biotype Y BPH.The results suggest that YHY15 exhibited distinct responses under the two infestation scenarios.In response to infestation with biotype Y BPH, YHY15 seedlings exhibited signs of withering at 4 days and eventually wilted completely by 7 days.In contrast, seedlings infested with biotype 1 BPH remained healthy and continued to grow vigorously (Figure 1A).To evaluate the ability of YHY15 rice to affect the biology of BPH insects, we measured the weight gain and amount of honeydew produced by the two BPH populations.As expected, biotype Y exhibited a significantly higher weight gain ratio (mean = 41.4%)than biotype 1 (mean = −2.2%)(Figure 1B).In addition, biotype Y produced significantly more honeydew (mean = 16.49mg) than biotype 1 (mean = 0.57 mg) (Figure 1C).These findings indicate that YHY15 rice plants are highly resistant to biotype 1 BPH but are susceptible to biotype Y BPH, as previously reported (Yang et al., 2004;Jing et al., 2011;Guan et al., 2022).Overall, these results provide strong evidence for the contrasting responses of YHY15 rice to BPH with different levels of virulence, and highlight the efficacy of the Bph15 gene in conferring resistance against specific BPH strains.

Overview of the miRNA-and RNAseq results
To better understand the mechanisms underlying the differential resistance displayed by YHY15 rice to virulent and avirulent BPH, we performed RNA and miRNA sequencing analyses using leaf sheaths from YHY15 seedlings infested with either biotype 1 (RT) or Y (RY) BPH for either 6 (RT6, RY6) or 48 h (RT48, RY48).Un-infested rice plants were used as controls and named as R0.After constructing and sequencing the mRNA/ miRNA libraries, high-quality raw sequence reads were normalized and subjected to further analysis.
Out of the 738 identified miRNAs, miR396, miR167, miR166, miR162, miR159, miR156, miR820, miR408, miR1425, miR1862, miR444, and miR827 exhibited the highest relative abundance (Supplementary Table 2).Principal component analysis (PCA) of the miRNA-seq data revealed significant variation between the R0 group and the RY6 or RY48 group, indicating that biotype Y BPH infestation led to distinct fluctuations in the miRNA profiles (Figure 2A).Among the 50 differentially expressed miRNAs (DEMs) detected in R0 vs RY6, 13 were upregulated and 37 were downregulated (Figure 2B; Supplementary Table 3).As infection progressed, more DEMs (64) were identified in R0 vs RY48, 16 of which were upregulated and 48 were downregulated (Figure 2B; Supplementary Table 3).In total, 38 DEMs were shared between R0 vs RY6 and R0 vs RY48, with 12 unique DEMs detected at the earlier time point and 26 at the later time point (Figure 2C).These results suggest that many more miRNAs are involved in the interaction between YHY15 and biotype Y BPH than in the interaction between YHY15 and biotype 1 BPH.Among these, miR156, miR5076, miR1856, miR398, miR5072, miR5079, miR408, miR2873, and miR169 exhibited significant variation in amplitude (Table 1).Many of these miRNAs have been reported to play specific roles in plant developmental processes and biotic and abiotic stress responses (Sharma et al., 2015;Zhao et al., 2017;Lin et al., 2018;Liu et al., 2019;Liebsch and Palatnik, 2020;Li et al., 2020;Liu X. et al., 2021;Gao et al., 2022;Pachamuthu and Hari Sundar, 2022;Zhao et al., 2022).Their high abundance and variable expression patterns suggest that they may contribute to BPH resistance in YHY15 seedlings, and therefore subsequent analyses were conducted on these miRNAs.
In comparison, PCA of the RNA-seq data revealed distinct separation among all groups (Figure 2D; Supplementary Table 4).Through comparison of RNA expression levels across different groups, we identified 1027 and 1154 differential expression analysis (DEGs) in R0 vs RT6 and R0 vs RT48, respectively.Notably, a larger number of DEGs were observed in R0 vs RY6 and R0 vs RY48 (6817 and 6273, respectively), suggesting that infestation with biotype Y BPH results in more severe effects than infestation with biotype 1 BPH (Figure 2E).A Venn diagram was used to compare the expression patterns of DEGs in avirulent/ virulent BPH-infested rice at the early (6 h) and late (48 h) stages of infestation, and the identified DEGs were subjected to GO and KEGG enrichment analyses (Figure 2F; Supplementary Figure 1; Supplementary Table 5).We identified 730 overlapping DEGs in R0 vs RT6 and R0 vs RY6 at the early infestation stage, which were enriched in cytochrome P450 (KEGG), protein modificationrelated biological processes (BP, GO), extracellular region (CC, GO), and monooxygenase activity (MF, GO) (Figure 2F; Supplementary Figures 1A, B).Specifically, genes in R0 vs RY6 were enriched in ribosome, photosynthesis proteins, carbon fixation in photosynthetic organisms, and DNA replication proteins (Figure 2F; Supplementary Figures 2A, B).However, genes in R0 vs RT6 did not exhibit these enrichments.We identified 827 overlapping DEGs in R0 vs RT48 and R0 vs RY48 at the late infestation stage, which were enriched in photosynthesis proteins, porphyrin metabolism, zeatin biosynthesis, and cytochrome P450, suggesting disturbance to the photosynthetic system (KEGG) (Figure 2F; Supplementary Figures 3A, B).Additionally, genes in R0 vs RY48 were enriched in ribosome-related pathways, DNA replication proteins, photosynthesis proteins, and carbon fixation in photosynthetic organisms (Figure 2F; Supplementary Figures 4A,   B).These results indicate that infestation with biotype Y BPH results in damage to the photosynthetic system earlier than infestation with biotype 1 BPH, which is consistent with the phenotypic observations (Figure 1A).

Transcriptional responses in YHY15 rice to BPH infestation
To elucidate the genome-wide responses of YHY15 rice to BPH infestation, the transcriptome was analyzed over the entire infection time course.According to Mfuzz analysis, distinct temporal patterns were clustered into 20 groups (Figure 3; Supplementary Table 6).Notably, clusters 1 to 4 exhibited significantly lower levels of transcription in RY6/48 compared to R0, whereas this trend was less pronounced in RT6/48 (Figure 3; Supplementary Table 6).The genes in these clusters were subjected to GO and KEGG enrichment analyses to explore their functions.Genes in cluster 1 were enriched in photosynthesis proteins, carbon fixation in photosynthetic organisms, terpenoid backbone biosynthesis, and porphyrin metabolism (Figure 4A; Supplementary Table 7).Accordingly,   GO analysis confirmed that many of these genes were related to the photosynthesis biological process, as well as photosynthetic cellular components such as thylakoid, photosynthetic membrane, and photosystem (Supplementary Figure 5A; Supplementary Table 7).These results suggest that infestation with biotype Y BPH results in significantly greater physical impairment than does infestation with biotype 1 BPH, which was consistent with the growth observations (Figure 1A).
In clusters 5 to 10, gene expression levels increased significantly in the RY groups, but only slightly in the RT groups (Figure 3; Supplementary Table 6).Genes in cluster 6 (Figure 4B; Supplementary Figure 5B; Supplementary Table 7) and cluster 9 (data not shown) were enriched in membrane trafficking (KEGG), cell membrane systems (GO), and protein transporter (GO), suggesting that defense against BPH relies on membrane transport systems and cell secretion.In cluster 7, genes were enriched in glycosyltransferases, GTP-binding proteins, phenylpropanoid biosynthesis, and two-component system (Figure 4C; Supplementary Figure 5C; Supplementary Table 7).The phenylpropanoid pathway is involved in the production of various metabolites, including flavonoids, lignin, lignans, and cinnamic acid amide, among others (Dong and Lin, 2021).Genes in cluster 10 were enriched in ribosomes, DNA replication proteins, translation factors, and other transcriptional activity-related factors (Figure 4D; Supplementary Figure 5D; Supplementary Table 7).Overall, the transcriptional system, transmembrane transport system, exosome, and phenylpropanoid pathway are notably stimulated by BPH infestation, particularly infestation with biotype Y BPH, suggesting that they may be crucial for rice resistance to insects.
The expression patterns of several clusters appeared to be desynchronized, indicative of different responses to infestation by either biotype 1 or biotype Y BPH.In clusters 11 to 14, gene expression was increased at the early (6 h) and late (48 h) stages during infestation with biotype 1 BPH, but remained relatively stable during infestation with biotype Y BPH (Figure 3; Supplementary Table 6).These clusters were mainly enriched in zeatin biosynthesis, nitrogen metabolism, cytochrome P450, and plant hormone signal transduction (Figure 4E; Supplementary Figure 5E; Supplementary Table 7).In plants, cytochrome P450 catalyzes various primary and secondary metabolic reactions and is involved in the synthesis and metabolism of terpenoids, alkaloids, sterols, fatty acids, plant hormones, signaling molecules, pigments, flavonoids, and isoflavones, among others (Hansen et al., 2021).Clustering and time-course expression of genes after BPH infestation.The 20 distinct temporal gene expression patterns were computed with Mfuzz.The x axis represents the five treatment groups, and the y axis represents log 2 -transformed normalized intensity ratios for each group.R0, uninfected controls; RT6, YHY15 seedlings infested with biotype 1 BPH for 6 h; RT48, YHY15 seedlings infested with biotype 1 BPH for 48 h; RY6, YHY15 seedlings infested with biotype Y BPH for 6 h; RY48, YHY15 seedlings infested with biotype Y BPH for 48 h.
Such multifunctionality makes cytochrome P450 important in the plant defense against pests, diseases, and abiotic stressors.In addition, zeatin is the primary active component of the phytohormone CK.KEGG analysis revealed that the cytochrome P450-phytohormone signal transduction pathway was activated in YHY15 rice under biotype 1 BPH infestation.In contrast, gene expression in clusters 15 to 17 was increased at the early (6 h) and late (48 h) stages during infestation with biotype Y BPH (Figure 3; Supplementary Table 6).These genes were enriched in proteasome, membrane trafficking, steroid biosynthesis, protein processing in endoplasmic reticulum, exosome, and GTP-binding proteins (Figure 4F; Supplementary Figure 5F; Supplementary Table 7).
These results suggest that cell secretion and extracellular materials are crucial for defense against insects.Finally, clusters 18 to 20 did not exhibit any notable tendencies (Figure 3).

Phytohormonal responses in YHY15 rice to BPH infestation
Both the clustering analysis and KEGG/GO enrichment analysis indicated that phytohormone signaling was differentially affected by infestation with either biotype 1 or Y BPH (Figures 3, 4E).Consequently, we compared the expression of BPH-responsive genes in rice with those induced by exogenous phytohormone application in A. thaliana.According to Hormonometer analysis, a total of 15,603 A. thaliana orthologs of rice genes were selected for comparison (Supplementary Table 8).The results revealed that genes associated with JA-dependent signaling were significantly altered by infestation with biotype 1 BPH, but only moderately changed by infestation with biotype Y BPH (Figure 5).Genes associated with ABA exhibited a similar pattern, except for being negatively correlated with genes whose expression in A. thaliana was elicited 0.5 h after ABA treatment (Figure 5).In contrast, the expression of genes involved in the ET and BR pathways exhibited negative correlations with those responding to phytohormone application in A. thaliana, although with more moderate suppression in biotype Y BPH-infested rice (Figure 5).These results show that genes associated with the JA, ABA, ET, and BR pathways are the most responsive to BPH infestation, particularly in the case of biotype 1 BPH, indicating their crucial roles in the rice defense against insects.
Previous research suggests that JA is involved in BPH resistance (Xu et al., 2021).We found that the JA-dependent pathway was differentially induced by infestation with either biotype 1 or Y BPH, implying that JA plays a crucial role in rice defense (Figure 5).We further analyzed the expression of genes involved in the JA pathway, and found that most of the genes involved in the biosynthesis of JA and jasmonoyl-isoleucine (JA-Ile) (e.g., LOX and JAR) were significantly suppressed by BPH infestation, especially infestation with biotype Y BPH (Figures 6A, B).Accordingly, most JAZs, which inhibit JA signaling, were induced in RT6/48 and RY6/48 (Figures 6A, B).In addition, MYC2 was downregulated in YHY15 rice infested with biotype Y BPH, suggesting that the JA-dependent pathway was inactivated upon infestation with the virulent BPH (Ge et al., 2018).Finally, the expression of many TFs was altered in Comparison of transcriptomic phytohormone signatures between BPH-infested rice and phytohormone-treated Arabidopsis thaliana.Red shading indicates positive correlations and blue shading indicates negative correlations between BPH-infested rice and phytohormone-treated A. thaliana.MJ, methyl jasmonate; ACC, 1-aminocyclopropane-1-caroxylic acid (precursor of ethylene); ABA, abscisic acid; IAA, indole-3-acetic acid; GA3, gibberellic acid 3; ZT, zeatin; BR, brassinosteroid; SA, salicylic acid.R0, uninfected controls; RT6, YHY15 seedlings infested with biotype 1 BPH for 6 h; RT48, YHY15 seedlings infested with biotype 1 BPH for 48 h; RY6, YHY15 seedlings infested with biotype Y BPH for 6 h; RY48, YHY15 seedlings infested with biotype Y BPH for 48 h.Yu et al. 10.3389/fpls.2024.1366515Frontiers in Plant Science frontiersin.orgBPH-infested rice, especially rice infested with biotype Y BPH, including those belonging to the JA-responsive bHLH, ERF, WRKY, and MYB families (Supplementary Figure 6).These distinct response patterns imply that these TFs have functions in BPH defense, although further research is required to verify these functions.The expression profiles of the identified TFs can be found in Supplementary Table 9.

Integrated miRNA and mRNA transcriptional analyses
The functions of the identified miRNAs were initially evaluated by scanning the literature.It has been reported that miR156 negatively regulates BPH resistance in rice by altering the expression of genes related to JA biosynthesis and signaling (Ge et al., 2018).We observed that miR156 expression was significantly decreased in biotype Y BPH-infested rice.Consequently, WRKY53, several MPKs (JA biosynthesis repressors), and several JAZs (negative regulators of JA signaling) were upregulated (Table 2).These transcriptional alterations resulted in low JA expression and improved BPH resistance.According to other research, miR396 suppresses BPH resistance through the miR396b-growth regulating factor 8 (GRF8)-flavanone 3-dioxygenase (F3H)-flavonoid pathway (Dai et al., 2019).Therefore, we examined the expression of miR396, GRF8, and F3H.Mature OsmiR396a and OsmiR396b share the same sequence (Dai et al., 2019), and we observed that miR396a was downregulated by infestation with biotype Y BPH.In contrast, GRF8 and F3H-1 expression increased in RY6 and RY48.Notably, GRF1 and two flavanone 3-dioxygenase 2-like genes exhibited the same trend (Table 2).These results suggest that the miR396b-GRF8-F3H-flavonoid pathway was activated in biotype Y BPH-infested rice.Genetic experiments indicate that miR398b negatively regulates pathogen-associated molecular pattern (PAMP)-induced callose deposition (Li et al., 2010).As miR398 expression was significantly downregulated in RY6/48, we further analyzed the transcription of callose deposition-related genes.Surprisingly, a callose synthase gene was activated following biotype Y BPH infestation, while a series of callose degradation genes were suppressed (Table 2).These results suggest that infestation with biotype Y BPH may lead to callose deposition in rice.Taken together, it appears that the expression of these miRNAs and their regulation of target genes may contribute the ability of YHY15 rice to resist infestation with biotype Y BPH.

Callose deposition is activated by BPH infestation
The combined miRNA-seq and RNA-seq analysis indicated that callose deposition might play a vital role in the rice response to BPH infestation.Callose deposition prevents plant hoppers from ingesting phloem sap (Hao et al., 2008;Du et al., 2009).Specifically, callose deposition ensures that the phloem sieve tubes remain occluded.To determine the role of callose in the rice response to BPH infestation, the expression of genes related to callose synthesis and degradation was examined.Consistent with the RNA-seq

Discussion
The management of rice plant resistance to BPH infestation will be crucial for effective pest management.Understanding the defense strategies employed by resistant rice varieties against both avirulent and virulent BPH populations can provide valuable insights for developing effective pest control strategies.Transcriptomics, which can quantify changes in gene expression and associated regulatory mechanisms, can aid in unraveling these complexities.In plants, miRNAs are involved in various developmental processes and play significant roles in abiotic and biotic stress responses (Zhang et al., 2013;Kumar, 2014;Zhang et al., 2016;Li et al., 2017;Mangrauthia et al., 2017;Natarajan et al., 2018).Specifically, miRNAs regulate targeted gene expression by binding complementary sequences in mRNA molecules, resulting in degradation and/or inhibited translation (Bartel, 2009;Marin et al., 2010;Peng et al., 2014).Integrated mRNA and miRNA transcriptomics analyses have been used to identify miRNA-mRNA networks associated with developmental processes and insect defense (Tan et al., 2020;Li et al., 2021;Mei et al., 2021).However, little attention has been paid to defense strategies employed by resistant rice against differentially virulent BPH populations.In this study, we employed an integrated mRNA and miRNA transcriptomics approach to characterize the defense responses of resistant YHY15 rice (contains Bph15) to infestation with avirulent (biotype 1) and virulent (biotype Y) BPH.Our results revealed that YHY15 rice seedlings exhibited distinct responses under the two infestation scenarios (Figure 1).YHY15 rice is highly resistant to biotype 1 BPH and susceptible to biotype Y BPH, which is in accordance with previous reports (Yang et al., 2004;Jing et al., 2011;Guan et al., 2022).
The BPH-resistance gene Bph15 has been widely applied in rice breeding programs, although the molecular mechanisms underlying Bph15-mediated resistance remain unclear (Lv et al., 2014;Zheng et al., 2021) signaling, and pathogenesis-related proteins (Lv et al., 2014).In addition, 20 upregulated and 3 downregulated miRNAs were identified in resistant rice variety P15 (containing Bph15) compared to susceptible rice variety PC (recipient line) (Wu et al., 2017).Combined with the mRNA transcriptome data, the 67 potential targets of these miRNAs were related to resistance responses to avirulent BPH, including abiotic and biotic stimuli, regulation of plant hormones (GA, SA, ET, and CK), cellulose biosynthesis, amino acid biosynthesis, and protein folding (Cheng et al., 2013;Lv et al., 2014;Wu et al., 2017).Here, comprehensive analysis of miRNA-seq and mRNA-seq data shed light on the underlying mechanisms responsible for the contrasting responses of YHY15 rice to BPH biotypes with different virulence.
Plant defense against BPH is a dynamic and sophisticated process which involves many levels of organizational and functional complexity (Barah and Bones, 2015;Erb and Reymond, 2019).During feeding, the BPH stylet transiently punctures the epidermis and then penetrates the plant cell wall.Subsequently, the insect salivates into the cells and ingests the phloem sap (Hao et al., 2008).According to electronic penetration graph (EPG) waveform recordings, BPH begin feeding on phloem sap within 1-3 hours of settling on rice plants (Hao et al., 2008).A greater number of miRNAs were found to be upregulated in biotype 1 BPH-infested resistant P15 rice (a Bph15 introgression line) than in susceptible PC rice (recipient line) during the early infestation stage (6 h) when the plants had not yet been severely damaged (Wu et al., 2017).In another report, the inducible BPH defense responses (indicated by upregulated DEGs and DEMs) were more robust during the early feeding stages (e.g., 6, 12, and 24 h) in resistant BPH6G rice (BPH6-transgenic rice) than in susceptible Nipponbare rice (wild type, WT) (Tan et al., 2020).Moreover, a miRNA profiling was conducted on resistant IR56 rice (carrying Bph3) under separate infestations of a virulent IR56-BPH and an avirulent TN1-BPH (Nanda et al., 2020).This study revealed that BPH feeding caused significant alterations in miRNA expression profiles of IR56 rice, with a greater number of miRNAs showing downregulation when IR56 rice was infested with TN1-BPH.However, the distinct mechanisms underlying rice plant responses to BPH of varying levels of virulence remains unclear.Here, resistant YHY15 rice plants were exposed to avirulent (biotype 1) and virulent (biotype Y) BPH.Notably, DEMs were only identified in rice plants infested with biotype Y BPH (Figure 2A).In addition, many more DEGs were identified in rice plants infested with biotype Y BPH, regardless of the infestation time (Figure 2E).Together, these results suggest that biotype Y BPH elicit more intense defense responses in YHY15 rice.DaEGs related with cytochrome P450 (KEGG) were more common at the early infestation stage (6 h) in rice infested with both types of BPH (Figure 2F; Supplementary Figures 1A, B).At the later infestation stage (48 h), overlapping DEGs in R0 vs RT48 and R0 vs RY48 were enriched in photosynthesis-related pathways (Figure 2F; Supplementary Figures 3A, B).Specifically, DEGs in R0 vs RY6 were mainly enriched in photosynthetic organisms and DNA replication proteins (Figure 2F; Supplementary Figures 2A, B), indicating that damage to the photosynthetic system occurs earlier in biotype Y BPH-infested rice.Taken together, compared with biotype 1 BPH, infestation with biotype Y BPH results in more serious damage and induces more intense transcriptional responses in YHY15 rice.
In plants, including rice, phytohormone signaling is widely known to be involved in insect defense and resistance (Liu et al., 2016;Jing et al., 2023;Liu Q. et al., 2021).Among phytohormones, JA is critical to the regulation of plant defenses against insect herbivores (Thaler et al., 2012;Howe et al., 2018).Other phytohormones, such as SA, ET, BR, ABA, and CK are also involved in plant responses to herbivory through cross-talk with JA signaling (Bruessow et al., 2010;Schäfer et al., 2015;Pan et al., 2018;Ma et al., 2020;Li et al., 2022;Zhang et al., 2022).In this study, we found that infestation with differentially virulent BPH populations resulted in noticeable effects on the expression of genes associated with phytohormone signaling in YHY15 rice plants.Specifically, infestation with biotype 1 BPH induced JA-and ABA-related signaling pathways, but suppressed ET-and BRrelated signaling (Figure 5).While, infestation with biotype Y BPH induced more moderate expression of JA-and ABA-related pathways.These results suggest that although infestation with virulent biotype Y BPH results in extensive damage to YHY15 rice plants, they exhibit weaker JA-and ABA-mediated defense responses.These findings are in line with previous research suggesting that IR56 BPH can overcome Bph3-mediated resistance by suppressing the transcription of defense-responsive MAPK pathways, phytohormone biosynthesis, and secondary metabolite production (Nanda et al., 2018).Previous studies also suggest that oral secretions, digestive and detoxifying enzymes, and endosymbionts can help insect herbivores adapt to host plants (Simon et al., 2015;Yates and Michel, 2018).We therefore speculate that specific BPH effectors may interact with Bph15 and affect Bph15-mediated immunity in rice.However, further studies will be required to test this hypothesis.
The DEMs and their associated mRNAs identified in this study may potentially play roles in the response of rice to BPH infestation.Among the identified DEMs, miR156, miR5076, miR1856, miR398, miR5072, miR5079, miR408, miR2873, and miR169 exhibited significant variation in amplitude (Supplementary Table 3).In addition, many of them have been implicated in various plant developmental processes and stress responses (Li et al., 2017b) (Sharma et al., 2015;Zhao et al., 2017;Lin et al., 2018;Liu et al., 2019;Liebsch and Palatnik, 2020;Li et al., 2020;Liu X. et al., 2021;Gao et al., 2022;Pachamuthu and Hari Sundar, 2022;Zhao et al., 2022).The altered expression of these miRNAs may enhance BPH resistance in YHY15 rice.Combined with mRNA transcriptome data (Table 2), we found that the miR156-JA, miR396b-GRF8-F3H-flavonoid, and miR398b-callose deposition pathways may contribute to the resistance of YHY15 rice to biotype Y BPH.Physiological tests verified the increased deposition of callose in biotype Y BPH-infested rice (Figure 7).Notably, these miRNAmediated responses only occurred in rice infested with biotype Y BPH, which suggests that they may contribute to the differential resistance of rice against BPH populations with varying levels of virulence.
Finally, research suggests that numerous TFs are involved in the rice response to BPH infestation.For example, OsMYB30 (an R2R3 MYB TF) induces the expression of phenylalanine ammonia-lyase (PAL) enzymes, thereby improving BPH resistance in rice (He et al., 2020).OsMYB22 promotes rice resistance by affecting flavonoid biosynthesis (Sun et al., 2023).The bHLH protein MYC2 is involved in JA-mediated insect resistance (Schweizer et al., 2013;Xu et al., 2021).OsWRKY45, OsWRKY53, OsWRKY70, and OsWRKY89 also mediate herbivore resistance (Chujo et al., 2007;Wang et al., 2007;Hu et al., 2015;Li et al., 2015;Huangfu et al., 2016).A recent study has demonstrated the pivotal involvement of OsWRKY71 in Bph15-mediated resistance (Li et al., 2023).Here, we identified several differentially expressed TFs associated with BPH infestation.Notably, several predominant TF families, including bHLH, MYB, ERF, WRKY, bZIP, NAC, C2H2, TALE, G2-like, HD-ZIP, MYBrelated, HSF, and NF-Y, were differentially responsive to BPH infestation.The expression of most TFs was altered in biotype Y BPH-infested rice, indicating their specific roles in defense against the virulent biotype.Notably, the expression of certain other TFs was disturbed specifically in response to biotype 1 BPH infestation.These findings highlight the important functions of TFs in BPH defense and warrant further research to uncover their specific roles.

Conclusion
In conclusion, our study provides valuable insights into the differential defense strategies employed by resistant YHY15 rice (carrying BPH resistance gene Bph15) against avirulent (biotype 1) and virulent (biotype Y) BPH.The BPH defense response was found to involve the modulation of miRNAs, TFs, phytohormone signaling pathways, and the induction of callose deposition.These responses were most noticeable in biotype Y BPH-infested rice plants.These findings contribute to the elucidation of the molecular intricacies underlying rice-BPH interactions and pave the way for further research into the specific genes, pathways, and regulatory elements involved in plant defense against diverse BPH populations.Studying these defense mechanisms will aid our understanding of the intricate interactions between rice and BPH and allow the development of targeted pest control strategies for improved rice cultivation.It is worth noting that the infestationinduced defense responses of YHY15 rice do not appear to alter the survivability of biotype Y BPH, implying the existence of corresponding adaptive responses in the virulent biotype.This result will be the subject of in-depth exploration in future studies.

7
FIGURE 7 Induction of callose deposition in YHY15 rice plants following BPH infestation.(A) Heat map of the expression of genes associated with callose deposition and degradation.Asterisks indicate statistically significant differences in gene expression at different time points in either biotype 1 or biotype Y BPH-infested rice relative to control (0 h) (*padj < 0.05, **padj < 0.01, via the Benjamini and Hochberg adjustment method).(B) Expression patterns of callose deposition-related genes in BPH-infested rice plants.UBQ was used as a control.(C) Callose deposition (CD) in BPH-infested YHY15 plants.The images were taken 48 h after BPH infestation.Samples were collected from the outermost sheath.Callose deposition values are the means of 20 biological replicates.Asterisks indicate statistically significant differences between BPH-infested and uninfested rice plants (one-way ANOVA, *P < 0.05).(D) Callose deposition in YHY15 plants infested with either biotype 1 or biotype Y BPH.Scale bar = 20 mm.The images were taken 48 h after BPH infestation.Uninfested YHY15 plants in the same condition served as control.R0, uninfected controls; RT6, YHY15 seedlings infested with biotype 1 BPH for 6 h; RT48, YHY15 seedlings infested with biotype 1 BPH for 48 h; RY6, YHY15 seedlings infested with biotype Y BPH for 6 h; RY48, YHY15 seedlings infested with biotype Y BPH for 48 h.

TABLE 2
The miRNA-mRNA interactions related to plant resistance.