Abstract Background Programmed cell death (PCD) can be triggered by biotic stress or abiotic stress in rice, which is a crucial factor throughout the lifecycle. Long non-coding RNAs (lncRNAs) have important roles in a variety of biological processes. However, the regulatory relationship between lncRNAs and rice PCD occurrence is unclear. Methods In this study, messenger RNAs (mRNAs) and lncRNAs associated with PCD regulation were identified in two rice cultivars: japonica rice (Oryza sativa) cultivar zheyoukang639 (zyk639) and its mutant zyk639-lesion mimic (zyk-lm). Besides, these data were used for weighted gene co-expression analysis (WGCNA) and construction of competitive endogenous RNAs (ceRNAs). Results As a result, 5,289 lncRNAs and 39,475 mRNAs were identified, of which 1,873 lncRNAs and 11,930 mRNAs were differentially expressed. In addition, WGCNA showed that the blue module was critical for the development of PCD in rice. The construction of the ceRNA regulatory network associated with the occurrence of PCD in rice identified 101 pairs of relationships, including 9 lncRNAs, 11 miRNAs and 76 mRNAs. Importantly, three major lncRNAs in zyk-lm, including MSTRG.24301.2, MSTRG.17476.1 and MSTRG.24898.1, may be involved in leaf PCD, with MSTRG.24898.1 binding to osa-miR580 to regulate transcription factors. Furthermore, the lncRNAs co-expression mRNAs rpl27-1 and OsRPS27, which were associated with ribosomal genesis, and Os11g0490900, a WRKY family transcription factor, may also regulate PCD genesis. Conclusions Models of protein synthesis regulation associated with zyk-lm PCD were discussed, lncRNAs may serve as key factors in regulating PCD and provide ideas for understanding PCD and genetic messaging. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-025-11844-y. Keywords: PCD, Whole-transcriptome, LncRNA, WGCNA, ce-RNA Introduction Programmed cell death (PCD) serves as the mechanism for cell demise in multicellular organisms triggered by internal developmental cues or external environmental signals at specific temporal and spatial coordinates [[52]1]. PCD is a common phenomenon in the whole process of plant growth and development. When plants begin to develop from embryonic cells, the formation and growth of organs such as roots, stems, leaves, flowers, and fruits are accompanied by cell PCD [[53]2]. In plants, PCD plays a critical role in various developmental events, such as cell death during xylogenesis, aerenchyma formation, several plant reproductive processes, leaf and petal senescence, and endosperm development [[54]3, [55]4]. PCD can be divided into three types: apoptosis, autophagy, and programmed cell necrosis according to the differences in morphology, regulation and execution factors [[56]5]. PCD is typically characterized by morphological changes such as condensation of nuclei and cells, the appearance of apoptotic bodies and genomic DNA fragmentation. Furthermore, the nucleus plays an important role in the process of PCD, the nucleoplasm is condensed, and the chromatin is dispersed at the edge of the nuclear membrane. DNA is subsequently fragmented to form apoptotic bodies under membrane enclosure [[57]6]. The regulatory mechanism of PCD in plants is complex and diverse, involving the interaction of multiple organelles in the cell. There were several signaling molecules involved in the regulation of cell death in the process of PCD in the plant cell system. At present, the known regulatory factors were Ca^2+, ROS, NO, and phytohormones [[58]7–[59]11]. PCD is a gene-regulated process controlled by related genes, proteases and transcription factors. The BAG gene in Arabidopsis has been implicated in the PCD pathway, and members ofBcl-2-related BAG family have been shown to play a central role in regulating autophagy and responding to biotic and abiotic stresses, which have diverse subcellular localizations in the cytoplasm, mitochondria, vesicles, nucleus and endoplasmic reticulum [[60]12, [61]13]. For example, the CaM/AtBAG5/Hsc70 signal transduction complex in mitochondria regulates mitochondrial membrane potential through Ca^2+ conductance and has an important role in regulating plant senescence [[62]14, [63]15].Transcription factors may have a switching role in triggering the development of PCD in plants under stressful environments, transcription factor families such as NAC, WRKY, and MYB have received widespread attention [[64]16, [65]17]. In previous studies, MAPK phosphorylation of the WRKY11 transcription factor induced RBOHB, which acts as a NADPH oxidase leading to a burst of PTI EOS and ETI ROS [[66]17]. Overexpression of OsAP25 and OsAP37, encoding aspartic proteases, leads to the development of PCD, and ETERNAL TAPETUM1, a conserved helix-loop-helix transcription factor (TF), directly regulates OsAP25 and OsAP37 expression in plants [[67]18, [68]19]. There were reports of an animal-like family of cysteoaspartase proteases in plant genomes, and plant metacaspases had strong evidence in Arabidopsis to support association with PCD [[69]20, [70]21]. Recent studies have shown that MicroRNAs (miRNAs) are regulators of many cellular processes and have been shown to modulate apoptosis by recognizing target sequences and interfering with transcriptional, translational, or epigenetic processes [[71]22]. Long non-coding RNAs (LncRNAs) are a class of endogenous non-coding RNAs that are more than 200 nucleotides in length and have been mainly used in epigenetic regulation, cell cycle regulation, and epigenetic processes [[72]23]. In the current study, lncRNAs mainly play important roles in epigenetic regulation, cell cycle regulation, cell differentiation and many other life activities [[73]24].According to the latest research, the interaction among messenger RNA(mRNA), miRNA, and lncRNA affects gene expression, so some scholars have proposed the competing endogenous RNA (ceRNA) hypothesis to explain the mechanism of gene expression regulation [[74]25]. It has been found that lncRNAs can act as a class of ceRNAs that bind to miRNAs to regulate the repression of target genes. For example, the lncRNA osa-eTM160 negatively regulates the repressive effect of osa-miR160 on osa-ARF18mRNA during the early stages of another development [[75]26]. In addition, ceRNAs have been used to reveal the regulatory mechanisms in plant biology, but the regulatory mechanisms in the ceRNAs network regarding the occurrence of PCD in rice are still unclear. Rice is an important food crop in the world, the quality and yield of rice is of paramount importance. According to estimates, the world will need to produce 25% more rice by 2030 to meet the challenges of feeding increasing populations. However, the various abiotic and biotic stresses are becoming serious global threats to sustained rice production [[76]27]. The presence of some disease resistance genes in Lesion mimic mutants (LMMs) formed through natural mutation or artificial variation can enhance rice resistance to pathogens. The investigation showed that a large number of cells died at the site of the lesions, which was a special type of PCD [[77]28, [78]29]. The further study of the mutant will be beneficial to the development of disease-resistant rice seed resources and the study of the molecular regulation mechanism of PCD [[79]30]. In the previous research, Ethyl Methane sulfonate (EMS) was utilized to induce mutations in the japonica rice (Oryza sativa) cultivar zheyoukang639 (zyk639), resulting in the generation of a mutant library. Among the mutants produced, a spontaneous lesion mimic mutant, designated zyk639-lesion mimic (zyk-lm), displayed distinct phenotypic and physiological PCD traits compared to its wild type zyk639. To clarify the relevant regulatory mechanisms of PCD in mutant rice, we performed whole-transcriptome sequencing of zyk639 and zyk-lm and identified differentially expressed RNAs, including miRNAs, mRNAs, and lncRNAs. In the present study, phenotypic analysis and physiological analysis combined with disease resistance analysis were conducted on zyk639 and zyk-lm, which provided a basis for further exploration of the regulatory mechanism of rice PCD occurrence through target gene function analysis and ceRNA network regulation analysis. Materials and methods Experimental materials The rice mutant zyk-lm and its wild cultivar zyk639 are both preserved by our laboratory in Zhejiang Academy of Agricultural Sciences. In mid-June 2024, seed soaking and sprouting treatments were conducted. After seed germination, the seeds were transplanted into seedling cultivation trays and raised in a rice-specific greenhouse. During the seedling stage, the temperature was controlled at 30 °C, the relative humidity was set at 75%, and the water level in the trays was maintained. When the rice seedlings reached the five-leaf stage, they were transplanted to an outdoor experimental field for further cultivation. The soil type was paddy soil. Throughout the rice growth period, conventional field cultivation practices were adopted for water and fertilizer management, with a topdressing of nitrogen fertilizer applied during the tillering stage. It was planted in randomised blocks with 20 cm line spacing and 35 cm column spacing. The first three leaves of two rice materials, wild-type zyk639 and mutant zyk-lm, were subjected to whole transcriptome analysis and named zyk639a, zyk639b, zyk639c, and zyk-lma, zyk-lmb, and zyk-lmc according to the sequential distribution of the flag, inverted one, and inverted two leaves. Three replicates were performed for each sample. Phenotypic, physiological and disease resistance analyses Firstly, the phenotypes of zyk639 and zyk-lm were observed by collecting the leaves at three different positions: the flag leaf, the second leaf and the third leaf, at the adult stage. A total of 10 strains of the mutant zyk-lm and zyk639 with consistent growth were randomly selected to investigate agronomic traits, such as plant height, 1,000-grain weight, and grain number per panicle. In stage-4 rice growth, 10 strains of Xoo(Xanthomonas oryzae pv. Oryzae, xoo) including P2, P3, P4, P5, P6, P7, P8, P9, C4 stored in our laboratory were inoculated into zyk639 and zyk-lm at 1 cm leaf tip by means of leaf clipping. After 30 days of inoculation, the disease situation of the plants was observed and recorded. The lesion length at three leaf positions of the two rice materials was investigated with 3 replicates. The leaves of zyk639 and zyk-lm at the stage of joint initiation, with a length of approximately 5 cm, were excised according to Bowling’s method [[80]31]. Subsequently, they were immersed in a pre-prepared solution of Trypan blue dye and subjected to boiling for a duration of 10 min. After cooling for 6 h, the samples underwent bleaching and rinsing using an acetic acid: anhydrous ethanol mixture (1:1) until decolorization was achieved. Finally, the decolorized samples were observed under an optical microscope. Referring to the staining method described by Wang [[81]32], the zyk639 and zyk-lm leaf samples were placed in 15 ml centrifuge tube with the upper end facing upwards. Subsequently, the pre-configured solution DAB dye solution was added to allow for infiltration of approximately one-fourth of the leaves. The samples were then subjected to light treatment at a controlled temperature of 25℃ for a duration of 10 h. Following the staining process, decolorization was carried out using acetic acid mixed with anhydrous ethanol in a ratio of 1:1, repeated several times until complete removal of excess dye was achieved. Finally, observation under an optical microscope was conducted after decolorization. RNA extraction and strand-specific library construction and sequencing Total RNA was extracted using Trizol reagent kit (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s protocol. RNA quality was assessed on an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) and checked using RNase free agarose gel electrophoresis. After total RNA was extracted, rRNAs were removed to retain mRNAs and ncRNAs. The enriched mRNAs and ncRNAs were fragmented by using fragmentation buffer and reverse transcribed into cDNA with random primers. Second-strand cDNA was synthesized by DNA polymerase I, RNase H, dNTP (dUTP instead of dTTP) and buffer. Next, the cDNA fragments were purified with the QiaQuick PCR extraction kit (Qiagen, Venlo, The Netherlands), end repaired, poly(A) added, and ligated to Illumina sequencing adapters. Then UNG (Uracil-N-Glycosylase) was used to digest the second-strand cDNA. The digested products were size selected by agarose gel electrophoresis, PCR amplified and sequenced using Illumina HiSeqTM4000 by Gene Denovo Biotechnology Co. (Guangzhou, China). To get high quality clean reads, reads were further filtered by fastp(version 0.18.0) [[82]33]. lncRNA-mRNA association analysis To reveal the interaction between antisense lncRNA and mRNA, the software RNAplex was used to predict the complementary correlation of antisense lncRNA and mRNA. The program contains ViennaRNA package, and the prediction of best base pairing was based on the calculation of minimum free energy through thermodynamics structure. The down-stream or 3’UTR region lncRNAs which had been previously annotated as “unknown region” were annotated again. lncRNAs in less than 100 kb up/down stream of a gene were likely to be cis-regulators [[83]34]. The cis target genes were then subjected to enrichment analysis of GO functions and KEGG pathways. We analysed the correlation of expression between lncRNAs and protein-coding genes to identify target genes of lncRNAs. Identification of differentially expressed mRNAs and LncRNAs Differentially expressed genes (DEGs) analysis was performed by DESeq2 software between two different comparison groups [[84]35]. The genes with the parameter of value P < 0.05 and fold change ≥ 1.5 (FC > 1.5) were taken as thresholds to be considered differentially expressed. In addition, GO enrichment analysis and KEGG enrichment analysis were performed for the identified differentially expressed messenger RNA (DEmRNAs). Weighted gene co‑expression network analysis (WGCNA) WGCNA is a systems biology method used to describe the gene association patterns between different samples [[85]36]. The co-expression relationship was established based on the link-connection matrix between lncRNA and mRNA, and lncRNA and mRNA with low expression levels were filtered out (FPKM < 1), the power of 8, TOM Type was unsigned, merged cut height of 0.85, and minimum module size of 50. The other parameters were defined as default values. Co-expression networks were constructed using the WGCNA (v1.47) package in R [[86]37]. The networks were visualized using Cytoscape (v3.7.0). Construction of ceRNA network We intersected mRNAs and lncRNAs in bule module with their DEmRNAs and DElncRNAs to locate the lncRNAs and targeted mRNAs, and then added miRNAs identified in the whole transcriptome to construct ceRNA network. The method of constructing ceRNA network is as follows: (1) Spearman rank correlation coefficient (SCC) was used to evaluate the expression correlation between mRNA-miRNA or miRNA-lncRNA. Selecting SCC > − 0.7 pairs as co-expressed lncRNA-mRNA pairs or mRNA-miRNA pairs, where miRNAs can differentiate mRNA and lncRNA simultaneously. (2) Pearson correlation coefficient (PCC) was used to evaluate the correlation between the expressions of lncRNA-mRNA. The pairs with PCC > 0.9 were selected as the co-expressed lncRNA-mRNA pairs, in which both the mRNA and lncRNA in each pair were targeted and co-expressed with common miRNAs. The networks were visualized using Cytoscape software. Quantitative real-time PCR To verify the differential genes in the RNA-seq data results, the OsActin-(122) gene was used as an internal reference gene for qRT-PCR verification, and primer sequences were retrieved from the NCBI database. Primer sequences were provided in Table S1. Relative expression levels of target genes were calculated using the ΔΔCt method, and primer amplification efficiencies were validated via standard curves to be between 95% and 105%. The primers were designed using primer premier 5.0 software and were synthesized by Beijing Kengke Xinye Biotechnology Co, which were listed in Table S1. The reaction mix consisted of the following: 10-µl final volumes containing 0.2 µl of cDNA, 0.2 µl of each Primer Premier, 5 µl of 2 × U1tra SYBR Mixture, and 4.4 µl RNasefree water. The PCR procedure was as follows: 94℃ for 5 min, followed by cycling for 30 rounds at 94℃ for 10 s, 60℃ for 10 s, and 72℃ for 20 s. The experiment was repeated three times using three replicates. Subcellular localization of the selected gene To determine the functions and expression patterns of the genes potentially involved in PCD of rice, we selected ELF3-2, PAC, PTD, and [87]Q2R432 for subcellular localization. The full-length coding sequences of ELF3-2, PAC, PTD, and [88]Q2R432 were amplified using the cDNA of Nip. The primers utilized for amplification are presented in Table S1. The PCR products were ligated to vector 580 by infusion enzymes and driven for expression by the 35 S: GFP promoter. The 35: GFP vector was chosen as the negative control. The constructed vectors were transfected into the protoplasts of rice stem cells, and the GFP signals of the selected genes were observed using a laser confocal microscope. Determination of reactive oxygen species, antioxidant enzyme activity and osmoregulatory substances The activities of superoxide dismutase (SOD), Peroxidase (POD), and catalase (CAT) were determined by nitrogen blue tetrazole method, guaiacol method, and potassium permanganate titration [[89]38, [90]39]. For SOD activity measurement, reagents included PBS solution (pH 7.8), 0.037 U/L xanthine oxidase, 1.25 mmol/L nitroblue tetrazolium chloride (NBT) solution, and 75 mmol/L xanthine solution. Absorbance was measured at 560 nm. For CAT activity detection, reagents included 0.02 mmol/mL H[2]O[2] standard solution and 50 mmol/L ammonium molybdate tetrahydrate solution. Absorbance was measured at 405 nm. For POD activity measurement, reagents included 100 mmol/L phosphate buffer solution (pH 6.0) and guaiacol. Absorbance of the working solution was measured at 470 nm at 30 s and 90 s. O^2− content and H[2]O[2] content were determined according to the method [[91]40, [92]41]. For O₂⁻ content determination, reagents included hydroxylamine hydrochloride solution, p-aminobenzene sulfonic acid solution, α-naphthylamine solution, and chloroform. 10 µmol/mL sodium nitrite standard solution was diluted into gradient standard solutions, which showed a characteristic absorption peak at 530 nm. The O₂⁻ content in samples was calculated using the standard curve equation. For H₂O₂ content determination, reagents included acetone, 50 mg/mL titanium sulfate, concentrated ammonia, 2 mol/L sulfuric acid, and a 1 µmol/mL H₂O₂ standard solution diluted into gradient standard solutions. Absorbance was measured at 415 nm to establish a standard curve equation for result calculation. Malondialdehyde (MDA) content was determined by thiobarbituric acid (TBA) colorimetry [[93]42]. For MDA content determination, 0.5% thiobarbituric acid was used, and MDA content was calculated using the difference in absorbance at 532 nm and 600 nm. Determination experiments were conducted at the six-leaf stage. The content of proline was determined by extracting sulfosalicylic acid [[94]43]. For proline content determination, reagents included acidic ninhydrin solution, 3% sulfosalicylic acid (extraction solution), glacial acetic acid, and toluene. A 5 mg/mL proline standard solution was diluted into gradient standard solutions. Absorbance was measured at 520 nm, and results were calculated by plotting a standard curve. The soluble sugar content was determined by anthrone colorimetric method and the soluble protein content was determined by Coomassie brilliant blue G-250 staining method [[95]39, [96]44]. For soluble sugar content determination, reagents included 6 mol/L HCl, 6 mol/L NaOH, and DNS reagent. 5 mg/mL glucose standard aqueous solution was diluted into gradient standard solutions. Absorbance was measured at 540 nm, and results were calculated by plotting a standard curve. For soluble protein content determination, Coomassie Brilliant Blue G250 solution was used. 5 mg/mL protein standard solution was diluted into gradient standard solutions, absorbance was measured at 595 nm, and results were calculated by plotting a standard curve. Wuhan ProNets Testing Technology Co. Ltd provided testing services. Results Phenotypic features of zyk-lm and zyk639 To compare the growth of these two rice cultivars, plant morphology analysis and statistics were carried out. The plant height of wild-type zyk639 was significantly higher than that of zyk-lm at the mature period (Fig. [97]1a, d). The 1000-grain weight of zyk-lm was significantly lower than wild-type zyk639(Fig. [98]1c). However, the productive panicles and grains per panicle of the zyk-lm mutant were higher than those of wild-type zyk639 (Fig. [99]1e, f). The leaves of the zyk-lm mutant showed pathogenesis-like characteristics from the seedling stage to the adult stage. There were many reddish-brown lesions on the third leaf of zyk-lm, which can also be observed on the second leaf, but there were almost no lesions on the flag leaves (Fig. [100]1b). Fig. 1. [101]Fig. 1 [102]Open in a new tab Comparison of plant morphology and agronomic traits between the wild type (zyk639) and mutant (zyk-lm). a Plant phenotypes at the adult stage. White scale bar = 20 cm. b Leaf phenotype at the adult stage. Black scale bar = 2 cm. Each set of images represents (from left to right), the three leaves on the left represent the flag leaf, the 2nd leaf (from the top) and the 3rd leaf of zyk639, the three leaves on the right represent the flag leaf, the 2nd leaf, and the 3rd leaf of zyk-lm. Black scale bar = 2 cm. c 1000-grain weight. d Plant height. e Number of productive panicles per plant. f Number of grains per panicle. Values are means ± SD of 10 plants, ** indicated significant differences at p < 0.01 Disease resistance and physiological analysis of zyk-lm Two rice materials, zyk-lm and zyk639, were inoculated with ten strains of xoo to observe and count the length of spots at each leaf position (Fig. [103]2a). We found that zyk639 as a susceptible material had longer spots of the two strains P6 and P10, and the spot lengths showed a tendency to increase in the order of inverted flag leaf, inverted two leaves, and inverted three leaves (Fig. [104]2b). On the contrary, zyk-lm did not show disease phenotypes under infestation with ten leaf blight species. Fig. 2. [105]Fig. 2 [106]Open in a new tab Comparison of disease resistance between the wild type (zyk639) and mutant (zyk-lm). a Bacterial blight lesion length of ten strains infected. b Infection phenotype at three leaf sites. Black scale bar = 2 cm Trypan blue (TB) staining on three leaf sites showed that the number of blue and black spots of zyk-lm increased compared with zyk639. At the same time, small blue spots were also found in the leaf surface area around the plaques, indicating that cell death occurred in the zyk-lm plaques and affected the surrounding cells, and the death phenomenon was aggravated with the aging of the leaves. But there were almost no blue spots in the three leaf positions of zyk639, indicating that there was no cell death in zyk639 leaves (Fig. [107]3a). The above results showed that cell death occurred at all three leaf sites of zyk-lm, which may be the cause of the pathological phenotype of zyk-lm. After rinsing and decolorization by 3,3’ -diaminobiphenyl (DAB) for 10 h, the leaves of zyk639 at all three leaf positions did not show obvious brown markings, while zyk-lm had obvious brown markings, and the inverted three leaves were larger than inverted two and inverted one (Fig. [108]3b). This is consistent with TB staining analysis, which may be due to the ROS accumulation in mutant zyk-lm causing cell PCD and thus the emergence of this phenotype. Fig. 3. [109]Fig. 3 [110]Open in a new tab DAB and Trypan blue (TB) staining on the leaves of zyk639 and zyk-lm sampled at the six-leaf stage. a DAB staining (from left to right) of the flag, second and third leaves of zyk639, zyk-lm. b TB staining (from left to right) of the flag, second and third leaves of zyk639 and zyk-lm Identification of DEmRNAs The leaves of 18 samples were sequenced by Illumina deep sequencing technology, which generated clean reads reached 99.28 −99.82% (Table S2). The findings demonstrated that the quality of transcriptome sequencing meets the requirements for subsequent analysis. Principal component analysis (PCA) also showed a close correlation in expression between replicate samples, confirming the reproducibility of the data and validating the subsequent analysis (Fig S1). We obtained 11,930 DEmRNAs based on the screening criteria P-value less than 0.05 and fold change ≥ 1.5 (Table S3). The total number of differentially expressed genes identified in the nine comparison groups were: zyk639a vs. zyk-lma, zyk639b vs. zyk-lmb, zyk639c vs. zyk-lmc, zyk639a vs. zyk639b, zyk639a vs. zyk639c, zyk639b vs. zyk639c, zyk-lma vs. zyk-lmb, zyk-lma vs. zyk-lmc, zyk-lmb vs. zyk-lmc: 4,478, 7,199, 7,170, 1,566, 2,307, 981, 6,640, 7,416, 1,779. Among them, zyk-lma vs. zyk-lmc had the largest number of DEmRNAs (Fig. [111]4a). Volcano map showed that there was a clear distribution pattern of up-regulated and down-regulated genes in different comparison groups (Fig. [112]4b). Fig. 4. [113]Fig. 4 [114]Open in a new tab Analysis of differentially expressed mRNAs in each comparison group. a The number of significantly regulated genes. b Volcano plots of mRNA expression levels Identification of DElncRNAs A total of 5,289 lncRNAs were obtained, including 95 known lncRNAs and 5,194 new predicted lncRNAs. According to the position of new lncRNAs relative to protein-coding genes on the genome, newly predicted lncRNAs can be divided into five categories (intergenic lnRNAs, bidirectional lncRNAs, intronic lncRNAs, antisense lncRNAs, sense overlapping lncRNAs) (Fig. [115]5a). A total of 1,873 DElncRNAs were obtained (Table S4), among all nine pairwise comparisons zyk639a vs. zyk-lma, zyk639b vs. zyk-lmb, zyk639c vs. zyk-lmc, zyk639a vs. zyk639b, zyk639a vs. zyk639c, zyk639b vs. zyk639c, zyk-lma vs. zyk-lmb, zyk-lma vs. zyk-lmc, zyk-lmb vs. zyk-lmc, group zyk639c vs. zyk-lmc had the highest number of DElncRNAs. Interestingly, except for zyk-lmb vs. zyk-lmc and zyk639a vs. zyk639b, the remaining seven groups had significantly more up-regulated genes than down-regulated genes (Fig. [116]5b). Fig. 5. [117]Fig. 5 [118]Open in a new tab Identification of differentially expressed lncRNAs (DElncRNAs). a Types of lncRNAs. b Number of significantly regulated genes of lncRNAs Correlation analysis between LncRNAs and mRNAs Using RNAplex to predict the complementarity relationships between antisense lncRNAs and mRNAs, we obtained 883 lncRNA-mRNA target gene pairs, which included 874 lncRNAs and 768 mRNAs (Table S5). Only 9 lncRNAs corresponded to 2 mRNAs, while the others maintained a one-to-one relationship. Additionally, lncRNAs with cis-regulatory functions were used to predict mRNA target genes, we identified 12,739 lncRNA-mRNA target gene pairs, encompassing 4,569 lncRNAs and 9,214 mRNAs (Table S6), where over 70% of the predicted lncRNAs target more than 2 mRNAs (Fig. [119]6). Next, we performed trans-regulatory target gene prediction analysis through correlation analysis of expression levels between lncRNAs and protein-coding genes across samples or co-expression analysis methods. From 18 samples of zyk639 and zyk-lm, we obtained 39,213 lncRNA-mRNA target gene pairs, comprising 549 lncRNAs and 4,500 mRNAs (Table S7). Among these lncRNAs, 52% targeted more than 10 mRNAs, and additionally, there were 976 mRNAs corresponding to only a lncRNA. Furthermore, we found that Os11g0490900, a WRKY family transcription factor, targets 53 lncRNAs, while MSTRG.3292.1 targets 2 R2R3-MYB family transcription factors, suggesting that lncRNAs associated with rice PCD may influence PCD occurrence through the regulation of transcription factors. Fig. 6. [120]Fig. 6 [121]Open in a new tab Analysis of cis-and trans-target genes of lncRNAs. a The predicted number of lncRNAs with cis relationship with mRNA. b The predicted number of mRNAs with cis relationship with lncRNA. c The number of lncRNAs with trans relationship with mRNA. d The number of mRNA with trans relationship with lncRNA Identification of conserved DEGs via weighted gene co-expression network analysis (WGCNA) To fully explore the potential regulatory relationship between lncRNA and PCD in rice leaves, we filtered the gene expression matrices of lncRNA transcripts and mRNA transcripts from wild-type zyk639 and mutant zyk-lm, obtaining a total of 16,485 genes for WGCNA (Table S8). It was divided into 18 modules via constructing a gene clustering relationship, and the gene expression profiles in each module were correlated with all samples (Fig. [122]7a). A heat map of that module sample matrix was generated (Fig. [123]7b). The gene expression profiles in the blue module were significantly different between samples with and without PCD (Fig. [124]7c). There were 5,870 genes in the blue module, which were mainly related to protein synthesis (Table S9). Furthermore, GO and KEGG enrichment analysis for co-expressed mRNAs targeted by lncRNAs in the blue module were performed. Various GO terms have been discovered for lncRNA-targeted co-expressed mRNAs, including “Cellular macromolecule metabolic process (GO0044260)” “nitrogen compound metabolic process (GO0006807)” “binding (GO0005488)” “cell (GO0005739)” “intracellular (GO0005622)” “organelle (GO0043226)” (Fig. [125]8a). In the KEGG pathway analysis, the main pathways involved in the blue module were “Ribosome (ko03010)” “Spliceosome (ko03040)” “RNA transport (ko03013)” “Ribosome biogenesis in eukaryotes(ko03008)” “DNA replication (ko03030)” (Fig. [126]8b). Fig. 7. [127]Fig. 7 [128]Open in a new tab WGCNA of genes identified in zyk639 and zyk-lm. a Tree graphs of DEGs by hierarchical clustering of topological overlapping dissimilarities. b Heat map showing the sample expression patterns in the modules. c The expression profile of all the co-expressed genes in module blue. The color scale represents the Z-score. The bar graph shows the consensus expression pattern of the corresponding co-expressed genes in this module Fig. 8. [129]Fig. 8 [130]Open in a new tab GO enrichment classification and KEGG pathway-enriched analysis for co-expressed mRNA targeted by lncRNA in the blue module. a mRNAs were significantly enriched in GO terms (P-value < 0.05). b KEGG pathways involving the top 20 terms. Q-value ranges from 0 to 1. The closer to zero the Q-value is, the more significant the enrichment is. A larger RichFactor value indicates a higher degree of enrichment To further understand and explore the relationship between these DEGs, we constructed gene co-expression network. We selected genes that were enriched in the blue module related to ribosome synthesis and DNA replication and drew the relationship between the differential genes with the weight value of the first 100 (Fig. [131]9). In this network, the key genes such as rpl27-1 (Os08g0404300) and OsRPS27,60s ribosomal protein L27 and 40 s ribosomal protein S27 was significantly enriched in ribosome pathway. Fig. 9. [132]Fig. 9 [133]Open in a new tab Co-expression regulatory network analysis of the bule module Analysis of lncRNA-related CeRNA networks in rice PCD occurrence In this study, the mRNAs, lncRNAs, and the selected DElncRNAs, DEmRNAs from the blue module were taken as the intersection and combined with miRNA to construct a ceRNA network using Cytoscape software. We identified101 relationships, including 9 lncRNAs, 11 miRNAs, and 76 mRNAs, and constructed 8 complex competitive relationship subnetworks (Fig. [134]10, Table S10). We identified three key lncRNAs: MSTRG.24301.2, MSTRG.17476.1 and MSTRG.24898.1, which can bind to miR8175-z, miR11033-z and osa-miR5809, respectively. MSTRG.21792.1 regulated Os02g0557000, Os03g0850400, Os06g0157000, Os08g0512400 in combination with novel-m0237-5p. MSTRG.5369.2 regulated Os01g0308300, Os02g0220500, Os03g0773300, Os06g0609700 in combination with miR477-x. MSTRG.31242.3 regulated with Os05g0393400, Os07g0666300 in combination with miR5225-z, miR3627-z. MSTRG.16622.1 regulated with Os03g0749800 in combination with miR5257-z. It is worth noting that the gene of Os01g0566100 can be co-regulated by MSTRG.7176.1 and MSTRG.28624.1 (Fig. [135]10). We found that Os09g0334500 (WRKY family TF) and Os09g0483600 (JMJ family TF) regulated target genes after binding MSTRG.24898.1 and osa-miR5809 (Fig. [136]10). Fig. 10. [137]Fig. 10 [138]Open in a new tab A CeRNA network constructed with mRNAs, lncRNAs, and the selected DElncRNAs, DEmRNAs from the blue module and miRNAs Validation of RNA-Seq data by qRT-PCR analysis A total of 18 genes were randomly selected for validation through qRT-PCR. The transcriptional level of the tested genes based on qRT-PCR correlated well with their respective abundance estimated by RNA-Seq, indicating that the RNA-Seq dates obtained in this study are reliable and can support the transcriptomic analysis presented above (Fig. [139]11). Fig. 11. [140]Fig. 11 [141]Open in a new tab The expression profiles of qRT-PCR were compared with the results of RNA-Seq analysis Subcellular localization analysis of selected genes By combining the predicted substitute genes with GFP in the protoplasts of rice 9311, we were able to clearly display the target genes’ location within the cell using a laser confocal microscope. We constructed four subcellular vectors for the DEGs PAC (Os05g0393400), PTD (Os05g0588200), ELF3-2 (Os01g0566100), and [142]Q2R432 (Os11g0490900) in rice protoplasts. Under the laser confocal microscope, we found that PTD was in the cytoplasm, ELF3-2were located in the nucleus (Fig. [143]12). These results suggested that multiple cellular organelles may be involved in triggering and regulating PCD. Fig. 12. Fig. 12 [144]Open in a new tab Subcellular localization analysis of selected genes Analysis of antioxidant enzyme system and osmoregulatory substances Many reports indicated that excessive accumulation of ROS and H[2]O[2] can cause PCD in plants [[145]45, [146]46]. Therefore, we inferred the corresponding antioxidant system in rice PCD occurrence by detecting the peroxide content, antioxidant enzyme activity, and osmotic regulation substance content in zyk639 and zyk-lm. As shown in Tables [147]1 and [148]2, there are significant differences in antioxidant enzyme activity and osmotic regulation substance content between different leaf positions in zyk639 and zyk-lm (p < 0.05). It is obvious that the SOD, POD, CAT, and MAD enzyme activity of zyk-lm was greater than that of zyk639, and the antioxidant enzyme activity of zyk-lmb was greater than that of zyk-lma and zyk-lmc. The content of O[2]^− and H[2]O[2] was basically consistent with the trend of antioxidant enzyme activity. Proline content of zyk-lma was 35.97 µg/g, which was significantly higher than other samples and was 222% of zyk-lmb and 213% of zyk639c. Interestingly, the soluble protein of zyk639 was greater than zyk-lm, and the soluble protein content of zyk639a was 144% of zyk-lmc, while the soluble sugar content of zyk-lmb was the lowest at 17.92 mg/g (Table [149]2). Table 1. The effect of zyk-lm on antioxidant activity Variety SOD (U/g) CAT (mmol/h/g) POD (U/g) H[2]O[2] (µmol/g) O[2]^− (nmol/g) MDA (nmol/g) zyk639 a 409.21 ± 15.48a 138.55 ± 3.71b 1618.31 ± 37.54c 5.75 ± 0.16b 69.18 ± 1.65c 15.15 ± 0.33c zyk639 b 296.13 ± 6.06d 134.58 ± 2.85b 2500.15 ± 4.94b 4.08 ± 0.09d 60.61 ± 0.26d 16.61 ± 0.30b zyk639 c 242.28 ± 10.82e 108.38 ± 4.39c 3749.96 ± 145.01a 2.67 ± 0.11e 76.56 ± 3.18b 15.51 ± 0.72c zyk-lm a 221.02 ± 9.37e 158.78 ± 0.11a 3770.45 ± 147.95a 7.73 ± 0.11a 65.94 ± 1.23c 15.92 ± 0.73bc zyk-lm b 377.45 ± 17.62b 151.59 ± 4.30a 3821.39 ± 72.77a 5.74 ± 0.07b 81.48 ± 3.13a 19.73 ± 0.15a zyk-lm c 349.79 ± 16.72c 140.55 ± 5,23b 3733.08 ± 24.29a 5.38 ± 0.11c 76.17 ± 1.00b 18.96 ± 0.74a [150]Open in a new tab Different lowercase letters indicate significant difference (P < 0.05) Table 2. Effect of zyk-lm on the content of osmoregulatory substances Variety Proline content (µg/g) Soluble protein content (mg/g) Soluble sugar content (mg/g) zyk639 a 18.51 ± 0.48c 16.52 ± 0.28a 28.28 ± 1.17a zyk639 b 17.80 ± 0.46 cd 13.28 ± 0.49bc 28.35 ± 0.83a zyk639 c 16.88 ± 0.52d 12.81 ± 0.35c 22.41 ± 1.07b zyk-lm a 35.97 ± 1.59a 12.07 ± 0.40d 30.08 ± 1.42a zyk-lm b 16.19 ± 0.53d 13.67 ± 0.40b 17.92 ± 0.51c zyk-lm c 21.45 ± 0.65b 11.50 ± 0.40d 21.96 ± 0.89b [151]Open in a new tab Different lowercase letters indicate significant difference (P < 0.05) Discussion We found that the leaves of the mutant rice zyk-lm showed a similar phenotype of lesions, and similar phenotypes occured in other plant species because of the occurrence of a PCD response. Therefore, the mutant zyk-lm will provide important material for further study on the regulation mechanism of rice PCD. At the adult stage, we got the agronomic traits of the wild type and the mutant. The plant height of zyk-lm was significantly lower than that of wild-type zyk639 at the mature period (Fig. [152]1a, d). The productive panicles and grains per panicle of the zyk-lm mutant were higher than those of wild-type zyk63, which indicated that the mutant had excellent cultivation advantages in agronomic traits. It is of great significance to explore its molecular regulation mechanism in PCD to reduce lodging and improve yield. Furthermore, the analysis of resistance to disease using zyk-lm and zyk639 showed that the length of the disease spots was significantly different, which further indicated that there were significantly different resistant genes in zyk-lm material compared to zyk639. The TB and DAB staining showed that the class pathological phenotype in zyk-lm might be caused by ROS accumulation in the plant, leading to cell PCD. The production of ROS is directly related to cellular nucleic acids, and the level of ROS is not only closely related to the degradation of rRNA and the destruction of ribosomes but also related to the oxidative stress response and the apoptosis pathway. This was consistent with the results of TB and DAB staining of zyk-lm in this study. In qRT-PCR analysis, BIP101 and BIP102, which encode proteins related to ribosome formation, were significantly up-regulated in the mutant zyk-lm. Thus, increased levels of these genes may lead to PCD of zyk-lm, and we speculate that genes related to ribosome formation affect the degradation of related nucleic acids, the translation of proteins, and the normal metabolism of cells in rice. Reactive oxygen species are produced in cells through aerobic metabolism, and when ROS exceed normal levels, they act in cells. Biological macromolecules, including DNA, are damaged when ROS accumulates continuously and reaches a certain death threshold. Apoptosis, autophagy and other programmed death pathways are activated, leading to cell death. In the analysis of antioxidant systems, our study showed that the antioxidant oxidase activity and peroxide content in the leaves of zyk-lm were significantly greater than that of zyk639, especially zyk-lm b. Our findings suggest that zyk-lm may offer a novel model to investigate regulatory mechanisms underlying rice PCD. LncRNAs, as the hotspot in biogenetics, have important roles in epigenetic regulation, cell cycle regulation, gene expression and many other life activities [[153]47]. However, the functional studies of the vast majority of lncRNAs are still unclear, and there is no systematic study to identify lncRNAs associated with PCD in rice. In the current study, by target gene prediction of newly predicted lncRNAs in the transcript species, there was a complex regulatory relationship between lncRNAs and mRNAs. Interestingly, transcription factors play a critical role in the occurrence of PCD in rice, one transcription factor can regulate multiple lncRNAs, and one lncRNA can also modulate multiple transcription factors. In the present study, we found that Os11g0490900, a WRKY family transcription factor, targets 53 lncRNAs, while MSTRG.3292.1 targets 2 R2R3-MYB family transcription factors, suggesting that lncRNAs associated with rice PCD may influence PCD occurrence through the regulation of transcription factors. WRKY transcription factors are a diverse group of transcriptional regulators that integrate plant responses to environmental stress and regulate development [[154]48]. Several WRKY TFs are involved in the regulation of cell death during biotic stress. In tobacco, WRKY1 was firstly identified as a positive regulator of HR PCD, following its phosphorylation and activation by the salicylic acid (SA) induced kinase SIPK [[155]49]. In addition, by WGCNA analysis of modules associated with PCD occurrence in rice, lncRNAs targeting their co-expressed mRNAs were involved in PCD-related biological processes and pathways. The results showed that some lncRNAs associated with rice PCD occurrence were regulated by miRNA binding to their co-expressed mRNAs influencing each other. We identified a total of 101 relationships, including 9 lncRNAs, 11 miRNAs, and 76 mRNAs, and constructed 6 complex competitive relationship subnetworks (Fig. [156]10a-e, Table S10). We identified three key lncRNAs: MSTRG.24301.2, MSTRG.17476.1 and MSTRG.24898.1, which can bind to miR8175-z, miR11033-z and osa-miR5809, respectively. Therefore, these results provided a molecular basis for the study of lncRNA regulation of PCD in rice. As a highly complex and evolutionarily conserved organelle, the ribosome plays an important role in the cell [[157]50]. A number of studies have shown that the nucleolus is not only involved in the maturation, assembly and export of ribonucleoprotein (RPN); it is also involved in the regulation of the cell cycle, cell senescence and cellular stress response [[158]51]. Under normal physiological conditions, ribosomes assemble proteins in an orderly manner in the nucleolus, and when stress occurs, the ribosomes become abnormal, cell cycle arrest and even apoptosis, a process called nucleolar stress [[159]52]. In the present study, many GO entries related to ribosome generation were found such as “ribosome (GO:0005840)”, “ribosome biogenesis (GO:0042254)”, “ribonucleoprotein complex (GO:1990904)”, “ribosomal subunit (GO: 0044391) “. Ribosomal stress is closely related to TFs under ribosome stress conditions in plant cells, RPs are released from the nucleolus into the nucleoplasm or RBFs mutate, there may be a p53-like transcription factor in plants that is like that in mammals. For example, the transcription factor NAC family SOG1 (suppressor of gamma response 1) is activated and phosphorylated by the PI3K-related protein kinases ATM and ATR. At this time, plant ribosome stress and DNA damage, as well as the NAC family transcription factor ANAC082, cause abnormal growth and development in plants [[160]53]. In the correlation analysis of lncRNA and mRNA, it was found that there was a complex regulatory relationship between lncRNA and TFs. DNA replication represents a crucial process in transmitting genetic information from parental generation to offspring. Errors that accumulate during DNA replication and repair can cause genomic instability and disrupt the cell cycle [[161]54].Ultimately, this can lead to PCD. Furthermore, researchers discovered that sustaining a state of DNA hyperexcitation can significantly decrease the rate of DNA replication and cell proliferation. This strategy may alleviate metabolic stress and inhibit the onset of PCD [[162]55]. Increasing evidence shows that the stability of plant DNA structure and transcription levels correlates with epigenetic modifications [[163]56]. Epigenetic modification can regulate the function of organisms through multiple pathways. However, the inhibition of the expression of an expression-modifying enzyme results in an important effect on the growth of the plant, PCD is even stimulated to cause the death of plants [[164]55, [165]57]. Plant histones undergo rapid and reversible changes in chromatin when subjected to stress. Currently, histone acetylation changes the conformation of nucleosomes and increases the binding capacity of DNA templates to transcriptional regulators, it can successfully activate transcription, resulting in the condensation of rDNA chromatin, this finding was confirmed in maize aleurone layer 45 s rDNA and Arabidopsis thaliana HDA6 mutants [[166]57, [167]58]. However, the molecular mechanism of epigenetic regulation of PCD remains to be further studied [[168]59]. Interestingly, a sharp change in chromatin conformation was observed upon treatment with the transcription inhibitor actinomycin D, and the degree of 5’-terminal deficiency of transcription accumulates, accompanied by DNA methylation, and the level of histone H3 decreases. Histone acetylation levels are elevated, suggesting that epigenetic alterations are associated with the failure of 45 S rDNA condensation [[169]60]. In this study, it was shown that most of the genes involved in encoding the epigenetic modification enzymes, enriched to “N-methyltransferase activity(GO:0008170)”, “iysine N-methyltransferase activity(GO:0016278)” in the GO, for example, qRT-PCR results showed that the FIP genes (LOC4331675, LOC4343938) of the m6a methyltransferase family were expressed in the mutant zyk-lm. We predicted that changes in these kinases may be involved in the epigenetic modification that activates zyk-lm. Consistent with our previous KEGG analysis, DEGs are involved in lipid metabolic pathways, and lipoxygenase-encoding genes respond. Lipoxygenase promotes the degradation of rice cells, and one of the major products of this process is ROS, which leads to the accumulation of ROS. Thus, further facilitating PCD. Our analysis showed that there was significant differential gene enrichment in spliceosome as the key passway in zyk-lm leaves. Certain splitters engage directly with introns based on their sequence, one or more interacting proteins alter the folding of the intron RNA, which enhances cleavage. Alternatively, they may change the conformation of the intron RNA through indirect interactions with the RNA tether. For instance, RNC1 and WTF1 have been demonstrated to collaboratively facilitate the splicing of most group II introns within chloroplasts. Subsequently, ZmWHY was discovered to engage with CRS1 through an indirect RNA-mediated mechanism to enhance the splicing of atpE [[170]61, [171]62]. The complex protein composition of SR proteins and RNA-dependent ATPases, or helicases, as spliceosomes has received much attention. SR proteins are characterized by the presence of one or more RNA structural domains (RBDs), or RNA recognition sequences (RRMs), and highly enriched serine/arginine alternating repeats, they have been shown to be involved in promoting pre-mRNA participation in splicing [[172]63]. RNA-dependent ATPases and/or helicases play a key role in regulating spliceosomal RNA-RNA interactions and spliceosomal conformational transitions [[173]64]. The selection of the splicing site of the spliceosome is the key to the splicing of the plant. Pre-mRNA, the selection of 3′ cleavage sites and the enrichment of continuous A and U sequences are essential for the function of the pre-messenger spliceosome. It also determines the accuracy and efficiency of pre-mRNA cleavage [[174]65, [175]66]. In our analysis, GO analysis showed that lncRNAs associated with rice PCD were involved in “spliceosomal complex (GO: 0005681)”, “helicase activity (GO:0004386)”, “spliceosomal tri-snRNP” complex (GO: 0097526). KEGG showed that PCD-associated lncRNA targeting involves mRNA with “ATP-dependent RNA helicase UAP56/SUB2(Ko12812)”. In conclusion, the functional expression of the spliceosome ensures the conformational change of mRNA precursors and the nuclear-cytoplasmic transport of RNA. It is obvious that lncRNA may regulate PCD generation by regulating the key role of spliceosome in the synthesis of new proteins. Based on the above results, we propose a hypothetical working model for triggering PCD in zyk-lm (Fig. [176]13). In this model, the upregulation of ribosome-related gene expression affects nucleic acid degradation and protein translation in rice. Which, in turn, leads to PCD. In addition, DNA replication, mRNA splicing, and RNA transport may play important roles in the rice mutant zyk-lm. First, errors in sequence accumulation during DNA replication can be regulated by epigenetic modifications. The lipid metabolism pathway is closely related to epigenetic modification. Lipoxygenase can trigger the degradation of cells, leading to the accumulation of ROS and PCD. In addition, under the action of the splicing body, the accuracy of the splicing site ensures the expression of eukaryotic genes. It contributes to RNA processing and transport, plays a key role in cell fate decision-making, and is closely related to PCD. Finally, RNA transport depends on a series of carrier proteins, and abnormal RNA transport may lead to dysregulation of apoptosis-related gene expression. Thereby affecting cell apoptosis. Taken together, DNA replication, pre-mRNA splicing, and RNA transport form a coherent process in biology. DNA replication ensures the stable transmission of genetic information, while pre-mRNA splicing realizes the diversity and complexity of genetic information. RNA transport ensures that this information is transmitted accurately. These processes together constitute the basis of the flow and expression of genetic information in organisms. It is not only of great significance for maintaining the diversity of organisms but also determines the growth and development processes of organisms. Fig. 13. [177]Fig. 13 [178]Open in a new tab The putative working model for triggering PCD Conclusions In this study, there were 1,873 DElncRNAs and 11,930 DEmRNAs in 9 pairwise groupings. In addition, WGCNA identified lncRNAs associated with the occurrence and presentation of PCD in rice. GO enrichment results in the bule module showed that lncRNAs targeted mRNAs co-expressed by lncRNAs involved biological processes such as ribosome generation, DNA replication, RNA transport, and splicing, we propose a hypothetical working model for triggering PCD in zyk-lm (Fig. [179]13). Finally, 6 ceRNAs network of related DElncRNAs involved in rice PCD was constructed by DEmiRNAs regulating DEmRNAs expression, Therefore, we speculate that MSTRG.24301.2, MSTRG.17476.1, MSTRG.24,898 1 and MSTRG.24898.1 and Os11g0490900, a WRKY family transcription factor, may be the crucial factors involved in the regulation of PCD occurrence in rice. Supplementary Information [180]12864_2025_11844_MOESM1_ESM.docx^ (22.6KB, docx) Supplementary Material 1: Table S1. List of primer sequences of qRT-PCR. [181]12864_2025_11844_MOESM2_ESM.xlsx^ (197.5KB, xlsx) Supplementary Material 2: Table S2. The leaves of 18 samples'reads stat. Table S3. The obtained 11930 DEmRNAs. Table S4.The obtained 1873 DElncRNAs. [182]12864_2025_11844_MOESM3_ESM.xlsx^ (3.6MB, xlsx) Supplementary Material 3: Table S5. antisense analysis between lncRNA and mRNA. Table S6. cis-regulatory analysis between lncRNA and mRNA. Table S7. Trans- analysis between lncRNA and mRNA [183]12864_2025_11844_MOESM4_ESM.xlsx^ (389.2KB, xlsx) Supplementary Material 4: Table S8. lncRNA transcripts and mRNA transcripts from wild-type zyk639 and mutant zyk-lm for WGCNA. Table S9. All expressed genes in bule module. Table S10. Identified a total of 101 relationships for CeRNA network. [184]Supplementary Material 5.^ (381.7KB, tif) Acknowledgements