Identification of Rice LncRNAs and Their Roles in the Rice Blast Resistance Network Using Transcriptome and Translatome

1. Introduction

Long noncoding RNAs (lncRNAs) are transcripts longer than 200 nucleotides with limited or no protein-coding potential [1]. Emerging evidence indicates that lncRNAs participate in a wide range of biological processes across various organisms, regulating gene expression through diverse mechanisms at the chromatin, transcriptional, posttranscriptional, translational, and posttranslational levels [2,3]. In plants, lncRNAs have been shown to regulate key physiological activities, such as flowering time [4,5], crop yield [6], fruit development [7], photomorphogenesis [8], gene silencing [9], and responses to biotic and abiotic stresses [10,11,12].

Despite these advances, the functional characterization of plant lncRNAs remains relatively limited, with only a few regulatory mechanisms fully elucidated. For example, the lncRNAs COLDWARP, COLDAIR, and COOLAIR mediate vernalization by silencing FLOWERING LOCUS C (FLC) [9,13,14]. The long noncoding RNA LRK Antisense Intergenic RNA (LAIR) interacts with OsMOF and OsWDR5 to promote the expression of leucine-rich repeat receptor kinase (LRK) gene clusters, significantly increasing rice yield [6]. Importantly, the crucial roles of lncRNAs in plant responses to biotic stress and immunity are gradually coming to light. In tomato, lncRNA16397 induces the expression of SlGRX, reducing reactive oxygen species (ROS) accumulation and thereby increasing resistance to Phytophthora infestans [11]. In Arabidopsis, the lncRNA ELENA1 interacts with the mediator complex subunit MED19a to regulate PR1 expression, bolstering immunity against pathogens [10]. In rice, the lncRNA ALEX1 enhances resistance to bacterial blight by modulating the jasmonic acid (JA) pathway [15], whereas the long noncoding RNA SABC1 helps balance plant immunity and growth by regulating salicylic acid (SA) synthesis [16]. Moreover, lncRNA23468 functions as a competing endogenous RNA (ceRNA) that suppresses miR482b accumulation, consequently elevating the expression of NBS-LRR genes and fortifying defense against P. infestans in tomato [17].

Rice (Oryza sativa L.), one of the world’s most important staple crops (particularly in Asia and Africa), faces a major threat from rice blast disease caused by the fungus M. oryzae. Globally, rice blast can lead to yield losses up to 30%, posing a severe challenge to food security and agricultural economies, with estimated annual losses of up to $66 billion [18]. Hence, understanding the molecular mechanisms of rice blast disease is of paramount importance for global food security [19,20].

Earlier research demonstrated that protein-coding genes play significant roles in the defense response of rice against M. oryzae [21,22]. However, the involvement of lncRNAs in this process has not been extensively explored. A recent rice telomere-to-telomere genome annotation revealed an additional 1,373 protein-coding genes, indicating that the genome annotation is still incomplete, especially in intergenic regions [23]. Such gaps might contribute to false positives in lncRNA identification and thus hamper downstream functional analyses.

To more accurately identify and characterize rice lncRNAs expressed during M. oryzae infection, we combined multiple software tools and methods suitable for lncRNA identification in rice and developed a pipeline named RiceLncRNA (https://github.com/njausxl/RiceLncRNA). This pipeline begins by integrating rice translatome data to construct a novel rice protein-coding gene annotation database (the CodingRNA dataset). By merging this dataset with Michigan State University Rice Genome Annotation Project version 7 (MSU v7) [24], we substantially improved annotation completeness (adding approximately 18.9% more gene loci) and accuracy in intergenic regions. With this optimized annotation, we analyzed strand-specific RNA-seq data from IR25 (a monogenic line harboring the blast resistance gene Pikm), LTH (a susceptible variety), and NPB (a conventional variety) at multiple time points postinoculation (0 h, 12 h, 24 h, 72 h). Through systematic analyses of these data, we identified 9,003 high-fidelity lncRNAs, with 605 differentially expressed under inoculation—including 415 that respond directly to stress. Notably, the resistant variety IR25 exhibited 293 specific DELs, whereas LTH showed only 70, suggesting that lncRNAs are substantially involved in the resistance signaling network. Functional enrichment analysis revealed that the target genes of these DELs are significantly enriched in multiple hormones signaling pathways, notably SA, JA, ET, and IAA. By constructing a ceRNA regulatory network and performing WGCNA, we discovered that several core regulatory lncRNAs (e.g., LncRNA.9497.1, LncRNA.9562.1) either directly target resistance (R) genes, receptor kinases, and disease-resistance proteins, or indirectly regulate them by competitively binding specific miRNAs. Through these complementary mechanisms, these lncRNAs synergistically modulate hormone signaling pathways to improve resistance. This study not only refines the rice genome annotation and establishes a reliable lncRNA identification pipeline but also, more importantly, elucidates the regulatory network of lncRNAs in rice–M. oryzae interactions, providing essential theoretical foundations and potential targets for advancing plant immune mechanism research and implementing novel disease-resistant breeding strategies.

2. Results

2.1. Analysis of the CodingRNA Dataset Assembly

On the basis of translatome data, we constructed a rice protein-coding gene annotation database, which we termed “CodingRNA.” Compared with the MSU v7 reference genome alone, this translatome-based annotation covered approximately 57% of the genome, and merging both annotations yielded an additional ~18.9% of gene loci (Table S2). As shown in Figure S1, several protein-coding genes were identified within regions previously annotated as intergenic in MSU v7. For example, in the region from 4,643,010 to 4,659,242 on chromosome 4 (Figure S1A), MSU v7 had originally classified the entire stretch as intergenic; however, our detailed analysis of the corresponding translatome data revealed a novel protein-coding gene with clear ribosome footprint signals (labeled “CD48328”). This discovery confirms that the locus encodes a functional protein rather than a non-coding transcript. In summary, by incorporating such findings into the CodingRNA dataset, we significantly enhance the accuracy and completeness of rice genome annotation by ensuring that protein-coding transcripts are not erroneously classified as lncRNAs.

2.2. Integrated Translatome and Transcriptome-Based Lncrna Identification Pipeline and Analysis

To systematically uncover rice lncRNAs, we analyzed strand-specific RNA-seq data from IR25 (carrying the blast resistance gene Pikm), LTH (susceptible), and NPB (a conventional variety), sampled at different postinoculation time points (0 h, 12 h, 24 h, 48 h, 72 h) under M. oryzae. Overall, 40 RNA-seq libraries that have been validated for strand-specificity were utilized (Table S3), though two (LTH-24 h-2, LTH-24 h-3) were discarded due to poor quality.

After trimming low-quality reads and adapters, we retained approximately 10,070,404,673 high-quality reads with an average GC content of 48.89% (Table S4). The mean alignment rate to MSU v7 was 96.89%, with an 86.23% rate of correctly paired reads.

We eventually assembled 114,267 transcripts (≥200 nt, fragments per kilobase of transcript per million mapped reads (FPKM) ≥0.5). To mitigate false positives, an updated rice lncRNA identification pipeline (Figure 1A) was implemented. Its key innovation lies in employing the refined CodingRNA dataset to avoid misclassifying protein-coding transcripts as lncRNAs. Through GFFcompare, we extracted transcripts labeled with classcodes i, x, o, u, and p (totaling 12,836 candidates; see Figure S2A for an explanation of the classcodes). Subsequent filtering against protein families database (Pfam), RNA families database (Rfam), and NCBI nonredundant protein database (NR) removed additional protein-coding or known noncoding RNAs, giving 11,543 candidates (Figure 1B).

Moreover, we applied Coding Potential Calculator 2 (CPC2) (11,748 transcripts), Predictor of Long non-coding RNAs and mEssenger RNAs based on an improved K-mer scheme (PLEK) (11,133), and Coding-Non-Coding Index (CNCI) (10,098) to evaluate coding potential (Figure 1C). Only transcripts consistently deemed “noncoding” by all three tools were retained (9,350). Taking the intersection of these results, we ultimately identified 9,003 high-fidelity lncRNAs (Figure 1D, Table S5). Classification indicated 2,803 intergenic, 2,076 antisense, 419 bidirectional, 3,056 intronic, and 649 sense lncRNAs (Figure 1E, Table S6).

A more conventional pipeline that did not integrate translatome data identified 10,464 putative lncRNAs, which included the aforementioned 9,003 high-confidence lncRNAs plus an additional 1,461 transcripts (Figure S2B). Further analysis revealed that approximately 96.3% of these 1,461 transcripts exhibited protein-coding potential, mostly located in insufficiently annotated intergenic regions. When we compared these 1,461 transcripts against newly predicted genes in the T2T rice genome (Table S7), about 10.1% (147 transcripts) aligned with those newly predicted coding genes. These findings underscore that incorporating translatome data substantially reduces the false-positive rate in lncRNA identification.

The final set of 9,003 lncRNAs was unevenly distributed across chromosomes (Figure 2A, Figure S3A), with Chr1 bearing the most (1,177) and Chr10 the fewest (549). These lncRNAs featured a mean GC content of 41%, notably lower than the ~53% typical of protein-coding genes (Table S8). Approximately 64% spanned 200–400 bp (median 333 nt) (Figure 2B), much shorter than mRNAs (median 1,008 nt), and displayed lower expression levels (Figure 2C). Notably, ~84.75% contained only one exon (Figure 2D), and they had fewer isoforms than mRNAs (Figure 2E).

Further comparison to publicly available rice lncRNA databases (PlantNATdb, PNRD, RNAcentral, NONCODE, CANTATAdb, GreeNC) revealed that 3,778 (41.96%) of our identified lncRNAs matched known records (E value <1e-10, identity >80%, coverage >50%), of which 2,446 were shared across all four sample sets (Figure S3B).

In summary, the optimized identification pipeline that integrates translatome and transcriptome data not only significantly reduces the false-positive rate and ensures high accuracy of the results, but also provides a convenient workflow and reference standard for future rice lncRNA identification, thereby greatly expanding the rice lncRNA resource.

2.3. DELs and Their Differential Expression (DE) Target Genes Regulatory Networks Reveal the Mechanisms Underlying Rice Response to M. oryzae Infection

To investigate the transcriptional regulation differences in rice varieties under M. oryzae infection and to discover potential resistance- or susceptibility-related lncRNAs, we systematically analyzed the DE of both mRNAs and lncRNAs in the resistant variety IR25 and the susceptible variety LTH at 0 h, 12 h, and 24 h postinoculation. A comprehensive comparison of transcriptomic differences between inoculated and noninoculated conditions and across time points revealed 4,426 differentially expressed mRNAs (DEGs), including 3,821 DEGs and 605 differentially expressed lncRNAs (DELs). Of these, 415 DELs and 3,338 DEGs were stress-responsive (Table S9, Figure 3A).

When comparing conditions before and after inoculation, the resistant genotype IR25 displayed 293 lineage-specific DELs, whereas the susceptible genotype LTH presented only 70 (Figure 3B). Further analysis indicated that IR25 had the highest number of specific DELs (137) at 12 h postinoculation, suggesting a faster transcriptional response during the early stage of M. oryzae stress. By applying cis and trans target gene prediction (20 kb upstream and downstream, |r| > 0.5; RIblast combined with correlation analysis), we identified 61 cis and 341 trans target genes for the 293 IR25-specific DELs (Table S10, Table S11). GO enrichment analysis showed that these target genes are primarily involved in photosynthesis, “small molecule biosynthetic process” (GO: 0006508), and “organic acid metabolic process” (GO: 0016053)—all of which are closely tied to aromatic amino acid metabolism (Figure 3C). KEGG pathway analysis further revealed significant enrichment in the “phenylalanine, tyrosine, and tryptophan biosynthesis” pathway (osa00400) (Figure 3D), with notable enrichment of genes such as LOC_Os01g55870 (chorismate mutase 3, chloroplastic) and LOC_Os09g08130 (indole-3-glycerol phosphate synthase). These genes play critical roles in phenylalanine and tryptophan metabolism, respectively, serving as key precursor nodes for SA and IAA biosynthesis. These results suggest that the DELs specific to IR25 may enhance plant defense against M. oryzae by modulating aromatic amino acid metabolic pathways and expediting the synthesis of resistance-related hormones.

In contrast, the susceptible genotype LTH exhibited only 70 specific DELs, which were associated with 18 cis and 74 trans target genes (corresponding to 15 and 12 DELs, respectively). GO enrichment showed that these target genes were enriched in carbon metabolism-related terms such as “starch metabolic process” (GO: 0019250) and “sugar metabolic process” (GO: 0006006) (Figure 3E). Meanwhile, KEGG analysis indicated that “starch and sucrose metabolism” (osa00500) and “glyoxylate and dicarboxylate metabolism” (osa00630) pathways were significantly enriched (Figure 3F). Thus, under M. oryzae stress, LTH-specific DELs seem more inclined to govern carbon metabolism reprogramming rather than potent defense hormone pathways, suggesting that LTH predominantly undergoes a metabolic adjustment-based stress pattern at early stages, failing to trigger strong hormone-mediated defenses in a timely manner.

Furthermore, we detected 52 DELs shared by IR25 and LTH (Figure 3B), all of which were consistently up- or downregulated in both varieties. GO enrichment analysis on the cis and trans target genes of these 52 DELs uncovered terms related to “rhythmic process” (GO: 0048511) and sugar metabolism, such as “hexose metabolic process” (GO: 0019318) (Figure S4A). Their KEGG enrichment highlighted the “pentose phosphate pathway” (osa00030), “Calvin cycle carbon fixation” (osa00710), and “amino acid biosynthesis” (osa01230) (Figure S4B). These findings point to a conserved role for these commonly responsive genes in energy metabolism and fundamental physiological regulation, helping coordinate basic adaptive responses in both resistant and susceptible genotypes.

We further examined DELs in the NPB variety at 24, 48, and 72 h post-inoculation (Figure S5), observing only a few DELs—just two overlapped with the aforementioned set of 52 DELs—indicating significant differences in how LTH, IR25, and NPB respond at the lncRNA level. This highlights genotype-specific defense mechanisms.

A direct comparison of the susceptible and resistant genotypes (LTH vs. IR25) at 12 h and 24 h postinoculation revealed 21 lncRNAs co-responsive to M. oryzae infection (Figure S4C). Their cis/trans target enrichment suggested potential synergy between JA and ET biosynthesis, as well as IAA modulation. GO enrichment analysis indicated that phospholipase A1 (PLA1) and cysteine synthase genes are closely associated with these 21 lncRNAs, with PLA1 being crucial for JA biosynthesis during early infection stages and cysteine serving directly as a precursor for ET (Figure S4D) [25,26]. Furthermore, within the KEGG-enriched “plant hormone signal transduction” pathway (Figure S4E), two key genes—JAZ (LOC_Os04g32480) and SAUR (LOC_Os02g52990)—were identified (Figure S4F-G) [27,28]. The marked upregulation of JAZ suggests a potential negative feedback mechanism isn JA signaling, whereas the pronounced downregulation of SAUR may indicate suppression of IAA, thereby prioritizing the synthesis of defense-related hormones.

We showed that the resistant variety IR25 exhibits a more robust transcriptomic response to M. oryzae infection than the susceptible LTH, largely through lncRNAs involved in aromatic amino acid metabolism and hormone signaling. Moreover, we identified a subset of lncRNAs that coordinate JA, ET, and IAA pathways, indicating that JA and ET signals play pivotal roles under M. oryzae stress while IAA signaling functions as an auxiliary route, collectively maintaining a dynamic balance between stress tolerance and immune responses.

2.4. Construction of a ceRNA Network Reveals the Role of lncRNAs in Rice Blast Resistance

This study constructed a ceRNA network containing 20 miRNAs, 17 lncRNAs, and 115 mRNAs (Figure 4), systematically analyzing how lncRNAs indirectly regulate resistance genes and hormone signals by “sponging” miRNAs during M. oryzae infection. Four lncRNAs—LncRNA.9497.1, LncRNA.9562.1, LncRNA.13491.1, and LncRNA.33800.3—emerged as major nodes, closely linked to hormone signaling pathways (e.g., JA; IAA; ABA; and gibberellin, GA) as well as the expression of various resistance-related genes.

LncRNA.9497.1 indirectly regulates several hormone-related genes by “sponging” osa-miR395a, osa-miR2864.1, and osa-miR5830. For instance, osa-miR395a targets OsSultr2;2 (a sulfate transporter), potentially bolstering sulfur metabolism to supply more substrates for defense [29]. In addition, osa-miR2864.1 targets a variety of receptor kinases, including OsLRK6 (leucine-rich repeat receptor kinase), OsMRLK16 (wheat germ agglutinin domain kinase), OsRLCK204 (receptor-like cytoplasmic kinase), and SDRLK-40 (receptor-like kinase). Among these, OsLRK6 is vital for immune signal transduction, whereas OsMRLK16, OsRLCK204, and SDRLK-40 facilitate microbe-associated molecular pattern (MAMP) perception and potentially abiotic stress cross-tolerance. Significantly higher expression of these receptor kinases in IR25 at 24 h post-inoculation suggests that certain lncRNAs may facilitate early-stage pathogen recognition through miRNA “sponging”.

Moreover, osa-miR5830 targets OsABA8ox3 and OsMETS2, which regulate ABA degradation and methionine metabolism, respectively. OsABA8ox3 modulates ABA levels, influencing stress adaptability [30]. In IR25 inoculated with M. oryzae, OsABA8ox3 was notably downregulated by 24 h, whereas OsMETS2 was upregulated, thus promoting ethylene biosynthesis and reinforcing disease resistance [31]. Thus, in the resistant IR25 line, LncRNA.9497.1 may indirectly downregulate ABA signaling and upregulate ethylene signaling via miRNA sponging, thereby enhancing the downstream immune response in rice.

LncRNA.33800.3 indirectly regulates OsWRKY70, OsCPS1, and OsGH3-2 by “sponging” osa-miR5075. OsWRKY70 is a key transcription factor in the JA signaling pathway that positively regulates JA biosynthesis while negatively affecting GA synthesis, thus prioritizing defense over growth [32]. OsCPS1 participates in GA biosynthesis [33], and OsGH3-2 is an IAA amino acid synthase implicated in broad-spectrum resistance by limiting auxin levels [34].

In LTH (susceptible) at 12 h postinoculation, OsGH3-2 was significantly downregulated, whereas OsWRKY70 was markedly upregulated, indicating that under early pathogen stress, lncRNAs such as LncRNA.33800.3 may simultaneously reduce auxin and raise JA signals to enhance disease resistance.

Furthermore, LncRNA.9562.1 and LncRNA.13491.1 regulate several resistance genes and hormone signals by targeting osa-miR529a, osa-miR5830, and osa-miR2090.

osa-miR2090 targets OsRLCK42 (receptor-like protein kinase) and Cht5 (chitinase). OsRLCK42 aids in early stress signal transduction, whereas Cht5 degrades fungal cell walls, a process pivotal to JA-mediated immunity [35]. In LTH, OsRLCK42 was strongly downregulated while Cht5 was upregulated, suggesting partial defense activation yet compromised early signaling in a susceptible background.

osa-miR529a targets OsABA8OX2 (ABA metabolic balance) and OsbHLH148 (a transcription factor regulating JA signaling), which contribute to both drought and disease stress responses [36] . In resistant IR25, OsABA8OX2 tended to be downregulated at 24 h postinfection, while OsbHLH148 was upregulated, reflecting a shift from ABA to JA-driven defenses [37,38].

osa-miR5830 targets OsABA8ox3 and interacts with OsMETS2, a methionine synthase supporting ethylene biosynthesis [31,39]. This further underscores intricate hormone crosstalk controlled by hub lncRNAs.

Collectively, these ceRNA interactions indicate that lncRNAs can coordinate hormone pathways (ABA, ET, JA, IAA) and receptor kinases to bolster rice blast resistance.

2.5. WGCNA Uncovers Key Hormone, R genes / Proteins, and Receptor Kinase Networks Under M. oryzae Stress

After removing batch effects, normalizing expression, and filtering outliers, a total of 2,583 lncRNAs and 15,701 genes were retained for WGCNA (Table S12). With a soft threshold power of 18 (scale-free topology index ~0.85; Figure 5A-B), 26 modules were identified (TOMType=Unsigned, deepSplit=2, minModuleSize=30, mergeCutHeight=0.2). Notably, around 85.23% of all genes/lncRNAs gathered into the top 10 modules (Figure 5C).

By correlating module eigengenes (MEs) with phenotypic traits (blast resistance vs. susceptibility) (Figure 5D), three notable correlations emerged. The Darkgreen Module showed strong positive correlation with resistance (r=0.96, p=3×10^-21), while the Lightyellow Module demonstrated negative correlation with resistance (r=–0.87, p=1×10^-12). The Gray60 Module exhibited negative correlation with susceptibility (r=–1.0, p=9×10^-40). Scatter plots of transcript significance (TS) vs. module membership (MM) (Figure 6E, Figure S6A-B) confirmed these relationships, prompting further network analysis of hub lncRNAs and their target genes.

In the darkgreen module (Figure 6A, Table S13), LncRNA.9497.1, LncRNA.21901.1, and LncRNA.35959.1 serve as core nodes, associating with disease/stress resistance genes such as wall-associated kinase 1 (WAK1), putative disease resistance protein (RGA4), OsWRKY125, and mitogen-activated protein kinase 17 (MAPK17). WAK1 mediates cell wall strengthening and defense signal transduction [40]. Meanwhile, RGA4 and RGA3 are NB-LRR-type resistance genes detecting pathogen effectors [41], and OsWRKY125 is a WRKY TF driving the expression of defense genes [42]. MAPK17 sits at a crucial position in MAPK signaling cascades under pathogen attack [43].

Within the lightyellow module (Figure 7C), LncRNA.10688.1 and LncRNA.36066.1 emerged as central regulatory points. LncRNA.10688.1 is directly linked to RGA5 (LOC_Os11g37740, LOC_Os12g37770) and Bph40 (LOC_Os11g39209), both contributing to pathogen-triggered immunity (PTI) and physical defense barriers [44,45]. It also relates to OsLP2 (LOC_Os12g08240), implying cross-disease and environmental stress crosstalk [46]. In LTH inoculated with M. oryzae, RGA5 and OsLP2 were notably upregulated, whereas Bph40 was slightly downregulated, suggesting partial immune enhancement yet weaker mechanical defenses in a susceptible background.

Piks-1 and Piks-2, recognized CC-NBS-LRR pairs defending against M. oryzae [47], and Pb1 (LOC_Os04g06280), a durable blast resistance gene from indica [48], also converged in this module’s coexpression network, implying multiple layers of early PTI and defense synergy. In terms of stress resistance, OsPUB69 (LOC_Os12g33180), an E3 ubiquitin ligase gene associated with LncRNA.21855.1, has been shown to play a crucial role in stress tolerance, particularly in protein degradation and cellular homeostasis regulation [49]. Similarly, OsGLP8-12 (LOC_Os12g28015), associated with LncRNA.36066.1, regulates reactive oxygen species (ROS) levels, which is essential for rice adaptation and disease resistance under oxidative stress conditions [50,51].

The gray60 module features LncRNA.13491.1, LncRNA.9562.1, and LncRNA.9997.2 (Figure 7B), which coordinate multiple receptor kinases (OsRLCK366, OsRLCK5), resistance genes (RGA5-L1), and hormone regulators (OsPP2C19). For example, LncRNA.13491.1 intersects with RGA4, RGA5, RPM1, and Xa1, all tied to effector-triggered immunity (ETI) [52,53]. It also influences OsMADS56 and OsbHLH179, potentially modulating strigolactone (SL) and GA signals [54,55]. Similarly, LncRNA.9562.1 correlates with OsRLCK366 (a receptor kinase) and OsPP2C19 (an ABA regulator) [56,57], facilitating hormone interplay and stress responses [58,59,60]. LncRNA.9997.2 impacts the regulation of GA, ABA, and even cytoskeletal elements via OsRLCK5 and VLN4 [61]. Interestingly, based on the expression heatmap of lncRNA within the module (Figure S6C-E), LncRNA.13491.1 and LncRNA.9562.1 peaked in LTH at 12 h but remained low in IR25, indicating potential genotype-specific temporal regulation of immune-associated lncRNAs.

Collectively, these findings demonstrate that multiple lncRNAs within WGCNA modules coordinate hormone pathways (ABA, GA, SL), receptor kinases (OsRLCK5, OsRLCK366), and R genes/proteins (RGA4, RGA5, Pb1, Piks-1, Piks-2) to bolster both early PTI and downstream ETI, thereby establishing a robust multi-layered defense framework in rice.

2.6. Mechanistic Insights into the Functional Roles of Key lncRNAs in Rice Blast Resistance

To examine the mechanistic roles of the top lncRNAs identified by WGCNA—e.g., LncRNA.9497.1, LncRNA.9562.1, LncRNA.10688.1, and LncRNA.13491.1—we performed gene set enrichment analysis (GSEA) for single key lncRNAs, focusing on pathways implicated in secondary metabolism and hormone signals, and structural diagrams for these four lncRNAs were also generated (Figure S7).

LncRNA.9497.1, annotated as an antisense transcript (Figure 7A), was significantly associated with the phenylpropanoid biosynthesis pathway (osa00940) and broader plant secondary metabolite biosynthesis (osa00999) (Figure 7A). These metabolic processes contribute to lignin and flavonoid production, which reinforce cell walls and produce antimicrobial compounds. Thus, LncRNA.9497.1 likely boosts disease resistance through robust metabolic defense.

Genes linked to LncRNA.9562.1 were mainly enriched in plant hormone signal transduction (osa04075) (Figure 7B, Figure S8), encompassing SA, ET, brassinosteroid (BR), GA, and JA pathways. Critical components like NPR1 and PR1 and ET-related genes appear among its targets, indicating a capacity to orchestrate hormone crosstalk and secondary metabolism, thereby enabling a rapid immune response to M. oryzae.

LncRNA.10688.1 was significantly enriched in the phenylpropanoid biosynthesis pathway (osa00940) and the carotenoid biosynthesis pathway (osa00906) (Figure 7C). Phenylpropanoids foster structural defenses and are precursors to SA, while carotenoids enhance antioxidant capacity and can feed into ABA production, thus regulating plant stress adaptation.

LncRNA.13491.1 exhibited notable enrichment in pigment biosynthetic processes (GO: 0046148) and flavonoid biosynthesis (osa00941) (Figure 7D). Flavonoids help mitigate oxidative stress and may act as signaling molecules during pathogen attacks, suggesting that LncRNA.13491.1 aligns plant defense responses with redox equilibrium.

Altogether, these key lncRNAs enhance plant adaptation and defense under M. oryzae infection by integrating basic metabolic, secondary metabolic, and multiple hormones signaling pathways—particularly the phenylpropanoid, flavonoid, JA, SA, ET, and IAA pathways. Their multifaceted roles highlight potential targets for breeding disease-resistant rice varieties using gene editing or advanced selection techniques.

2.7. qRT–PCR Analysis and LncRNA Cloning Validation

To validate these bioinformatic results (Table S14), four important lncRNAs and seven of their target genes were chosen from key WGCNA modules for quantitative real-time PCR (qRT–PCR) (Figure 8A–C, M, G–I, N-Q). The expression patterns largely matched the RNA-seq data, reinforcing the reliability of differential expression analysis. Another three random DELs and their three target genes were likewise confirmed (Figure 8D–F, J–L), supporting the computational pipeline’s accuracy.

One high-priority lncRNA (lncRNA.9562.1) was subsequently cloned, and Sanger sequencing verified its length and sequence fidelity (Figure S9A-B). The cloned sequence fully matched the predictions, indicating the pipeline’s robust annotation. These validation results pave the way for potential functional analyses, such as overexpression or knockout experiments, to assess the phenotypic impacts of this lncRNA on rice disease resistance.

3. Discussion

3.1. Optimized lncRNA Identification Pipeline and Its Contribution to Rice Genome Annotation

In this study, we developed a more comprehensive rice protein-coding gene database (CodingRNA) by merging translatome data with existing genome annotations, thereby refining the lncRNA identification pipeline (Figure 1A). Compared with conventional annotation approaches, our pipeline demonstrates an enhanced ability to exclude transcripts with protein-coding potential, effectively reducing false-positive lncRNA calls. This improvement aligns with recent telomere-to-telomere genome assembly work [23] reporting the discovery of 1,571 additional coding genes beyond MSU v7, underscoring the importance of updated annotation strategies (Table S7).

We identified 9,003 high-fidelity lncRNAs (Figure 1D, Table S5)—exceeding previous tallies by Wang et al. (4,787) and Priyanka et al. (5,337) [62,63]. Such discrepancies likely stem from differences in rice varieties, stress conditions, sequencing strategies, and lncRNA pipelines. Notably, 40% of our lncRNAs matched those in public databases (Figure S3B), whereas 60% appeared novel, highlighting possible genotype-, tissue-, or stress-specific expression. Our pipeline incorporates translatome data integration for excluding potential coding transcripts, three-pronged coding potential prediction (CPC2, PLEK, CNCI) supplemented by Pfam/Rfam/NR filtering, rice-specific TAD features for lncRNA target gene identification, and RIblast for trans-target detection, replacing older methods.

Despite substantial progress, there remain limitations. The focus on certain genotypes and time points leaves open questions about broader spectrum resistance and lncRNA behavior under other M. oryzae races or environmental contexts. Additionally, although qRT–PCR verified expression changes (Figure 8), functional assays such as overexpression or knockout are still needed to confirm direct phenotypic influences. Lastly, this work primarily addressed leaf tissue, calling for future multi-tissue or single-cell omics analyses to reveal dynamic spatiotemporal patterns of lncRNAs. Nevertheless, by offering an integrated pipeline on GitHub, this study establishes a rigorous foundation for more accurate lncRNA discovery and functional research in rice.

3.2. LncRNAs Participate in Rice Blast Resistance by Regulating JA, ET, and IAA Signaling Pathways in a Coordinated Manner

We observed that during M. oryzae infection, the target genes of DELs — which are themselves DEGs — in IR25 (resistant) and LTH (susceptible) varieties were significantly enriched in the JA and ET signaling pathways. This aligns with prior evidence highlighting the importance of the JA pathway in rice blast resistance. [62]. Yet, our study also unveils the ET pathway’s role. For example, some lncRNA targets in both IR25 and LTH—such as PLA1 (involved in JA biosynthesis) and PCO (affecting ethylene-related transcription factors [25,26])—suggest overlapping ET–JA crosstalk (Figure S4D). This synergy aligns with known JA–ET cooperation in defense against necrotrophic pathogens [64,65].

In IR25-specific responses, many DEL targets clustered within aromatic amino acid metabolism, implying potential enhancement of phenylalanine- and tryptophan-derived SA and IAA. Notably, previous studies have demonstrated that IAA and SA exhibit antagonistic effects, and this balance helps conserve energy during pathogen infection, thereby optimizing plant defense responses, which is consistent with our findings [66]. Meanwhile, JAZ genes (negative regulators of the JA pathway) [27] and SAUR genes (auxin-responsive [28]) also appear among the common DEL targets (Figure S4E-G), implying dynamic hormone crosstalk. Overall, IR25’s greater number of specific DELs (293) versus LTH’s 70 indicates a more robust or earlier activation of these hormone-regulated defenses in the resistant genotype. Conversely, LTH-specific DELs were concentrated in carbon metabolism, illustrating a metabolic shift that might be less effective at early defense induction.

Thus, multiple hormone signals—SA, ET, IAA, JA—likely collaborate under M. oryzae stress, and each genotype’s capacity to coordinate these pathways promptly could decide the strength and timeliness of immune responses.

3.3. lncRNAs Mediate Immune Responses by Regulating RLKs and R Genes/Proteins and by Participating in ceRNA Networks

Our results further indicate that, within the WGCNA modules correlated to blast resistance, RLKs and resistance proteins are markedly enriched. For instance, in the dark green module, LncRNA.13491.1 interacts with OsRLCK366 (a receptor-like kinase), OsRGA5-L1 (a resistance protein), and OsPP2C19 (a signaling regulator). This observation is consistent with findings in mulberry [67], where the lncRNA MuLRR-RLK-AS negatively modulates RLK expression to influence disease resistance.

Moreover, our ceRNA network showed that several lncRNAs—LncRNA.9497.1, LncRNA.9562.1, LncRNA.13491.1, and LncRNA.33800.3—competitively bind specific miRNAs (osa-miR395a, osa-miR2864.1, osa-miR5830), thereby influencing the expression of OsSultr2, OsWRKY70, and OsGH3-2, genes affecting sulfur metabolism, JA signaling, and auxin homeostasis, respectively. These discoveries validate the endogenous target mimic (eTM) hypothesis by Franco-Zorrilla et al. (2007) [68]. Similarly, studies in tomato [69] and sweet sorghum [70] suggest that lncRNAs strengthen resistance by regulating transcription factors and stress response genes via ceRNA-based interactions.

Furthermore, subcellular localization predictions (Table S15) provide further insights into how these lncRNAs may exert their regulatory functions. For instance, LncRNA.9562.1 exhibits high-confidence localization in extracellular, mitochondrial, and nuclear compartments, suggesting multifaceted roles in perceiving pathogen signals or orchestrating nuclear transcriptional events. Conversely, LncRNA.9497.1 is predicted to reside largely in the extracellular region, supporting a putative role in intercellular communication or apoplastic signaling. Such compartment-specific predictions enrich our understanding of how each lncRNA may spatially coordinate defense responses under M. oryzae infection.

In summary, lncRNAs in rice not only directly modulate RLKs and R genes/proteins but also indirectly shape hormone signaling pathways through ceRNA networks, forming an extensive regulatory web to determine the eventual disease resistance phenotype.

4. Materials and Methods

4.1. LncRNA Strand-Specific Library Data Sources

A total of 513 translatome (datasets were collected from various tissues of japonica rice under different experimental conditions [71,72,73,74,75,76], along with 40 strand-specific RNA-seq datasets derived from leaves of multiple japonica rice varieties [77,78] (Table S1, Table S3). These datasets served two major purposes: (i) constructing the rice coding gene dataset for improved annotation and subsequent lncRNA identification, and (ii) conducting differential expression analyses as well as coexpression network inference. Specifically, all 40 strand-specific datasets were used for lncRNA identification and WGCNA, whereas 16 of these datasets were devoted to differential expression analysis. We focused on the monogenic line IR25, harboring the blast resistance gene Pikm, and the susceptible line LTH (Lijiangxintuanhe) [79]. Among the 513 translatome datasets, various experimental replicates spanned stress conditions such as drought, temperature fluctuation, salinity, submergence, and heavy metals, in addition to multiple tissues (roots, stems, leaves, buds, panicles, and flowers). Detailed information on library construction is provided in Supplementary Materials 2.

4.2. LncRNA Identification and Classification

For the 40 strand-specific RNA-seq datasets, we first used FastQC (v0.11.9) to check read quality [80]. Strand specificity was determined via Rseqc (v5.0.1) [81], with any nonstrand-specific libraries excluded (Table S3). Adapter sequences and low-quality reads were removed using TrimGalore (v0.6.7) [82] (quality cutoff q30) and Fastp (v0.20.1) [83]. Bowtie2 (v2.4.2) was then employed to discard ribosomal RNA (rRNA) reads [84], after which HISAT2 (v2.2.1) [85] aligned the cleaned reads to the MSU v7 reference genome with strand-specific parameters (--rna-strandness FR or RF) [24]. Transcriptome assembly was carried out using Stringtie (v2.1.5) under the --rf or --fr options, and all resulting GTF files were merged using Stringtie merge. Transcripts below 200 nt in length or FPKM < 0.5 were discarded [86].

To remove protein-coding transcripts, we used GFFcompare (v0.9.8) to compare the assembled annotation against our integrated rice coding gene fusion annotation (hereafter referred to as the CodingRNA dataset) [87]. Candidate lncRNAs were those with class codes i, x, o, u, and p. Additional filtering steps were taken by comparing these candidates to the protein families database (Pfam) (v37.0) [88] and NCBI nonredundant protein database (NR) databases (e-value <1e-5) to eliminate any potential protein-coding transcripts [89], as well as against the RNA families database (Rfam) (v14.10) database (e-value <1e-5) to remove known small RNAs [90]. Next, three coding potential prediction tools—Coding Potential Calculator 2 (CPC2) (v1.0.1) [91], Predictor of Long non-coding RNAs and mEssenger RNAs based on an improved K-mer scheme (PLEK) (v1.2) [92], and Coding-Non-Coding Index (CNCI) (v2) [93]—were simultaneously applied, and only transcripts consistently predicted to be noncoding by all three tools were retained.

To further classify the final lncRNAs, we compared them once more with MSU v7 via GFFcompare. Class i was defined as an intronic lncRNA, class o as a sense lncRNA, class x as an antisense lncRNA, and classes p and u as general lncRNAs (Figure S2B). Additionally, bedtools closest (v2.30.0) identified bidirectional lncRNAs, while the rest were assigned as intergenic [94].

4.3. Identification of cis- and trans-Targets and Known lncRNAs

We extracted both differentially expressed genes (DEGs) and differentially expressed lncRNAs (DELs) from the expression analyses and used Pearson correlation analysis (the cor function in R) to evaluate DEL–DEG pairs with an absolute correlation >0.5 and p <0.05. For trans-target identification, we employed RIblast (interaction energy < –14 kcal/mol, interaction length ≥ 15 bp), alongside the same correlation threshold, to confirm putative lncRNA–gene associations [95].

Considering the median size of rice topologically associated domains (TADs) is ~35 kb [96], a ±20 kb window was chosen for cis-target searches. Any genes located within 20 kb upstream or downstream of the lncRNA, exhibiting |r| >0.5, were deemed cis-targets. We then downloaded six rice lncRNA databases—PlantNATdb (v1.4) [97], PNRD (v1.0) [98], RNAcentral (v22) [99], NONCODE (v6.0) [100], CANTATAdb (v2.0) [101], and GreeNC (v2.0) [102]—merging them via CD-hit (v4.8.1) [103]. Finally, we utilized Blastn (v2.9.0) [104] (E-value<1e-5, identity>95%) to compare our candidate lncRNAs with these known datasets, thus identifying novel lncRNAs versus those present in public repositories.

4.4. Differential Expression Analysis and Functional Enrichment

We performed differential expression analysis using the DESeq2 (v1.22.1) [105] package in R. Prior to analysis, batch effects were removed, and each group of lncRNAs and mRNAs was analyzed separately. Thresholds of |log2FC| ≥1.5 (for lncRNAs) or |log2FC| ≥2 (for mRNAs) and adjusted p ≤0.05 were used to designate significant differential expression. Principal component analysis (PCA) was then conducted via prcomp package, with visualization through ggplot2 (v3.5.1) [106] and pheatmap (v1.0.12).

Functional enrichment (Gene Ontology, GO; Kyoto Encyclopedia of Genes and Genomes, KEGG; Gene Set Enrichment Analysis, GSEA) was carried out using clusterProfiler (v4.10.1)[107], referencing the org.Osativa.eg.db (v0.01) rice database [108]. We applied a significance cutoff of p-value < 0.05 for all enrichment analyses.

4.5. Transcription Factor Identification and lncrna Localization

Transcription factors (TFs) are crucial in rice’s response to M. oryzae stress [20]. Hence, PlantTFDB (v5.0) [109] was utilized to predict TFs within the DEGs, the ceRNA network, and the WGCNA modules. LncRNA subcellular localization was assessed via RNALocate (v3.0) [110], providing an initial insight into nuclear- vs. cytoplasmic-located lncRNAs.

4.6. Weighted Gene Coexpression Network Analysis (WGCNA)

Weighted gene coexpression network analysis (WGCNA) was undertaken to explore interaction relationships among lncRNAs and mRNAs [111], excluding those already used for differential expression analysis. Following expression normalization (via DESeq2) and batch effect removal, low-expression data and outliers were discarded, retaining genes with higher variance (top 75% by median absolute deviation). A soft threshold power of 18 was used, constructing an adjacency matrix via the adjacency function. A topological overlap measure (TOM) was then calculated to define the similarity matrix of lncRNA–mRNA expression.

Hierarchical clustering was performed, and modules were defined or merged using the dynamic tree cut method (deepSplit=2, minModuleSize=30, mergeCutHeight=0.25). Each module’s eigengene was correlated with specific phenotypes (e.g., rice blast resistance vs. susceptibility), generating correlation matrices. Modules showing an absolute correlation coefficient |r|≥0.8 and p<0.05 were considered significantly associated with the trait. The core hub genes or lncRNAs within those modules were visualized using Cytoscape (v3.10.2) [112].

4.7. Single Key lncRNA Analysis

A normalized, batch-corrected expression dataset of both lncRNAs and mRNAs was used for Spearman correlation analyses. For each chosen lncRNA, all rice genes were sorted by the absolute value of their correlation (descending order), creating an ordered gene list.

Next, gene set enrichment analysis (GSEA) was performed via clusterProfiler (v4.14.4) on this ranked list to determine the biological processes or functional categories most associated with the lncRNA in question [107]. The ggplot, ridgeplot, and gseaplot functions in ggplot2 were utilized to visualize results from GSEA (v1.68.0) [113].

4.8. Competing Endogenous RNA (ceRNA) Network Construction

A total of 713 rice miRNAs were sourced from miRBase, forming the foundation for our ceRNA prediction library. psRNATarget (v2) [114], with parameter settings adapted from Zhang et al. [115], was used to predict lncRNA–miRNA and mRNA–miRNA interactions. The integrated mRNA–miRNA–lncRNA coexpression network was finally visualized in Cytoscape (v3.10.2).

4.9. qRT-PCR Method and lncRNA Cloning

To validate the reliability of the sequencing results, we randomly selected three long noncoding RNAs (lncRNAs), four key lncRNAs, and ten of their target genes for quantitative real-time polymerase chain reaction (qRT-PCR) analysis. Primers for qRT-PCR were designed using Primer3Plus software [116] and verified for specificity using PrimerBlast [117]. All primers were synthesized by Genscript Biotech (primer sequences are provided in Table S16).

qRT-PCR experiments were performed using the Applied Biosystems 7500 Real-Time PCR System, with three biological replicates for each sample. The reaction program was as follows: pre-denaturation at 98°C for 2 minutes, followed by 40 cycles of denaturation at 98°C for 2 seconds, annealing and extension at 59°C for 10 seconds. A melt curve analysis was performed after each run to confirm the specificity of the amplification products. The 18S rRNA gene was used as the reference gene. The relative expression levels of the target genes were calculated using the 2^(-ΔΔCt) method.

Specific primers (Table S16) for amplifying the full-length sequences of target lncRNAs were designed based on the RNA-seq data. Total RNA was reverse-transcribed into cDNA using SuperScript III Reverse Transcriptase (Invitrogen). PCR amplification was performed using PrimeSTAR GXL DNA Polymerase (TaKaRa), with the following reaction conditions: 98°C for 5 minutes; (98°C for 10 seconds, 60°C for 15 seconds, 72°C for 1 minute per kb) × 30 cycles; and 72°C for 10 minutes.

PCR products were separated by agarose gel electrophoresis, purified, ligated into the pNC-Cam1304-35S vector, and then transformed into DH5α competent cells. Positive clones were confirmed by Sanger sequencing (Tsingke). Sequencing results were aligned with the reference genome to verify the accuracy of the cloned sequences.

5. Conclusions

This study provides significant advancements in understanding the roles of lncRNAs in rice’s defense against M. oryzae. By creating an optimized lncRNA identification pipeline—incorporating translatome data and existing annotations—we identified 9,003 high-confidence rice lncRNAs with improved accuracy. In-depth analyses of differential expression, WGCNA, and ceRNA network construction revealed their critical involvement in multiple hormones signaling pathways (JA, SA, ET, IAA), as well as in regulating receptor-like kinases and resistance proteins. Key lncRNAs (such as LncRNA.9497.1 and LncRNA.9562.1) emerged as central regulators enhancing rice adaptability and immune response.

These findings offer a comprehensive framework for delineating lncRNA functions in plant immunity and propose molecular targets for breeding resistant rice varieties. Future efforts should emphasize the functional validation of top candidate lncRNAs through overexpression or knockout lines, exploring their applicability across diverse genotypes and M. oryzae strains.