Genome-wide identification and characterization of long non-coding RNAs involved in flag leaf senescence of rice

Key message This study showed the systematic identification of long non-coding RNAs (lncRNAs) involving in flag leaf senescence of rice, providing the possible lncRNA-mRNA regulatory relationships and lncRNA-miRNA-mRNA ceRNA networks during leaf senescence. Abstract LncRNAs have been reported to play crucial roles in diverse biological processes. However, no systematic identification of lncRNAs associated with leaf senescence in plants has been studied. In this study, a genome-wide high throughput sequencing analysis was performed using rice flag leaves developing from normal to senescence. A total of 3953 lncRNAs and 38757 mRNAs were identified, of which 343 lncRNAs and 9412 mRNAs were differentially expressed. Through weighted gene co-expression network analysis (WGCNA), 22 continuously down-expressed lncRNAs targeting 812 co-expressed mRNAs and 48 continuously up-expressed lncRNAs targeting 1209 co-expressed mRNAs were considered to be significantly associated with flag leaf senescence. Gene Ontology results suggested that the senescence-associated lncRNAs targeted mRNAs involving in many biological processes, including transcription, hormone response, oxidation–reduction process and substance metabolism. Additionally, 43 senescence-associated lncRNAs were predicted to target 111 co-expressed transcription factors. Interestingly, 8 down-expressed lncRNAs and 29 up-expressed lncRNAs were found to separately target 12 and 20 well-studied senescence-associated genes (SAGs). Furthermore, analysis on the competing endogenous RNA (CeRNA) network revealed that 6 down-expressed lncRNAs possibly regulated 51 co-expressed mRNAs through 15 miRNAs, and 14 up-expressed lncRNAs possibly regulated 117 co-expressed mRNAs through 21 miRNAs. Importantly, by expression validation, a conserved miR164-NAC regulatory pathway was found to be possibly involved in leaf senescence, where lncRNA MSTRG.62092.1 may serve as a ceRNA binding with miR164a and miR164e to regulate three transcription factors. And two key lncRNAs MSTRG.31014.21 and MSTRG.31014.36 also could regulate the abscisic-acid biosynthetic gene BGIOSGA025169 (OsNCED4) and BGIOSGA016313 (NAC family) through osa-miR5809. The possible regulation networks of lncRNAs involving in leaf senescence were discussed, and several candidate lncRNAs were recommended for prior transgenic analysis. These findings will extend the understanding on the regulatory roles of lncRNAs in leaf senescence, and lay a foundation for functional research on candidate lncRNAs. Supplementary Information The online version contains supplementary material available at 10.1007/s11103-021-01121-3.


Introduction
Non-coding RNAs are the main products of the eukaryotic transcriptome and possess key regulatory functions (Fabbri and Calin 2010). Long non-coding RNAs (lncRNAs) are a class of endogenous non-coding RNAs longer than 200 nucleotides which lack apparent protein-coding capacity (Zhu and Wang 2012). Based on the location in the genome, lncRNAs can be classified into five groups, including sense, antisense, intronic, bidirectional and long intergenic noncoding RNAs (lincRNAs) (Mattick and Rinn 2015). Compared with protein-coding genes, most lncRNAs are usually expressed at low levels and lack strong sequence conservation between species (Cabili et al. 2011;Necsulea et al. 2014). Increasing evidence has shown that lncRNAs play vital functions on growth, development, disease occurrence and epigenetic regulation in mammals (Sleutels et al. 2002;Rinn et al. 2007;Zhao et al. 2008).
With the rapid advances in high-throughput sequencing, numerous lncRNAs have been identified in plant species. For example, a total of 3857 lncRNAs have been identified during fruit ripening in melon, of which 1601 were differentially expressed between developmental stages and 142 are highly expressed (Tian et al. 2019). A comprehensive view of 833 high-confidence lncRNAs including 652 intergenic and 181 antisense lncRNAs in cassava have been identified under drought stress condition, of which 124 are drought-responsive . Additionally, 2572 lncRNAs related to flower development have been identified in Prunus mume, of which two lncRNAs XR_514690.2 and TCONS_00032517 might contribute to the formation of multiple pistils . In non-heading Chinese cabbage, a total of 4594 putative lncRNAs have been identified with a comprehensive landscape of dynamic lncRNA expression networks under heat stress (Wang et al. 2019a). Moreover, recent studies have also reported the functional studies of lncRNAs in plants. For example, overexpression of lncRNA T5120 in Arabidopsis promotes the response to nitrate, enhances nitrate assimilation, and improves biomass and root development . A total of 567 disease-responsive lncRNAs in rice have been systematically identified, among which, overexpression of lncRNA ALEX1 could activate jasmonate pathway and enhance resistance to bacterial blight (Yu et al. 2020). Taken together, these results all show the essential functions of lncRNAs in various biological processes, but not decipher in leaf senescence.
Emerging evidence has shown that lncRNAs could regulate the expression of protein-coding genes via numerous complex mechanisms to execute their functions (Liu et al. 2015a). For example, a novel antisense long noncoding RNA, TWISTED LEAF, maintains leaf blade flattening by regulating its associated sense gene R2R3-MYB in rice . A tomato lncRNA 16397 could regulate SlGRX gene expression to decrease ROS accumulation and weaken the injury of cell membrane, leading to increased resistance to Phytophthora infestans (Cui et al. 2017). In addition, a rice long non-coding RNA Ef-cd transcribed from the antisense strand of the flowering activator OsSOC1 locus can positively regulate the expression of OsSOC1, and plays roles in balancing grain yield with maturity duration (Fang et al. 2019). Overexpressing lncRNA LAIR transcribed from the antisense strand of neighbouring gene LRK (leucine-rich repeat receptor kinase) could increase grain yield and regulate the expression of LRK gene cluster in rice (Wang et al. 2018a). A rice cis-natural antisense long non-coding RNA cis-NAT PHO1;2 acts as a translational enhancer for its cognate mRNA and contributes to phosphate homeostasis and plant fitness (Jabnoune et al. 2013). Moreover, it has been reported that lncRNAs could act as a class of competing endogenous RNAs (ceRNAs) binding with microRNAs (miRNAs) to condense its repression on target genes (Zhu and Wang 2012). One example of ceRNA was found in Arabidopsis, and the study showed that lncRNA IPS1 could influence the expression level of PHO2 by binding to miR399 (Franco-Zorrilla et al. 2007). In addition, 89 lncRNA-originated endogenous target mimics for 46 miRNAs in tomato have been identified, of which lncRNA 23468 functions as a competing endogenous RNA to modulate NBS-LRR genes by decoying miR482b in the tomato Phytophthora infestans interaction (Jiang et al. 2019). A long non-coding RNA osa-eTM160 could attenuate the repression of osa-miR160 on osa-ARF18 mRNA during early anther developmental stages . Besides, ceRNA analysis also have been clearly illustrated in other species, such as honey bee , cucumber (He et al. 2019) and pepper (Zuo et al. 2019). However, ceRNA network has not been analyzed to study the functions of lncRNAs involved in leaf senescence.
Leaf senescence is a genetically programmed cell death process that constitutes the final stage of leaf development (Leng et al. 2017). Leaf senescence may limit yield in crop plants (Lim et al. 2007). Rice (Oryza sativa L.), a monocotyledonous model organism, is one of the major food crops. Thus, deciphering leaf senescence in rice is of great importance for understanding its molecular regulatory mechanism, and provides means to control leaf senescence for genetic improvement of yield. In the past decades, enough efforts to reveal the molecular mechanisms underlying leaf senescence have been mainly made by transcriptomic (mRNA), proteomic and metabolomic researches . Certainly, the identification of leaf senescence-associated genes (SAGs) and their functional studies were also performed, such as ONAC011 (El Mannai et al. 2017), SPL29  and OsCPK12 (Wang et al. 2019b). However, no systematic identification of lncRNAs involved in leaf senescence in rice has been reported.
In this study, a genome-wide high throughput sequencing approach was applied to investigate the genome-wide transcriptome changes of rice flag leaves during the developing process to senescence. Subsequently, senescence-related lncRNAs were systematically identified and characterized. Integrating the analysis on target gene functions and ceRNA networks, the putative regulatory function of these lncRNAs to mRNAs were predicted and analyzed. These results provided new insights into the regulatory mechanism of lncR-NAs in rice during flag leaf senescence, and laid a foundation for functional research on the candidate lncRNAs.

Plant material and sample collection
The standard laboratory rice cultivar 93-11 (Oryza sativa L.), an excellent maintainer line of indica rice in China, was used as the experimental material in this study. The 93-11 seeds were maintained in our laboratory by strict selfing.
The 93-11 plants were cultivated in paddy field by normal management way. The flag leaves at booting stage (FL1, 9 days before flowering), flowering stage (FL2, 3 days after flowering), early-senescence stage (FL3, 9 days after flowering), mid-senescence stage (FL4, 19 days after flowering) and late-senescence stage (FL5, 29 days after flowering) were collected, respectively (Fig. 1a). Three biological replicates were used for samples in each stage, with each sample pooled with three flag leaves from three independent plants, flash-frozen in liquid nitrogen, and stored at − 80 °C for subsequent analysis.

Detection of chlorophyll content and senescence marker genes
Chlorophyll was extracted from rice flag leaves with ice-cold 80% acetone, and the chlorophyll content per gram of leaf fresh weight (FW) was determined as previously described (Lichtenthaler 1987). Furthermore, the expression levels of senescence marker genes SGR (LOC_Os09g36200), RCCR1 (LOC_Os10g25030), OsNAP (LOC_Os03g21060) and Osl43 (LOC_Os01g24710) (Table S1) were performed by quantitative real-time polymerase chain reaction (qRT-PCR).

RNA extraction, library construction and sequencing
Total RNAs from samples were extracted with TRIzol reagent (Invitrogen, USA). Subsequently, the quality and integrity of total RNAs were detected by using Nanodrop 2000 and Agilent Bioanalyzer 2100 System (Agilent Technologies, USA).
After RNA extraction, a total of 1.5 μg RNA per sample was firstly used to remove ribosomal RNA using a Ribo-Zero rRNA Removal Kit (Epicentre, Madison, USA). Then, the sequencing libraries were constructed using NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina® (NEB, USA) according to the manufacturer's recommendations. Finally, the libraries were sequenced on an Illumina HiSeq 2500 platform, and 150 bp paired-end reads were generated.

Identification of lncRNAs and mRNAs
After RNA-Seq, clean reads were obtained by trimming adaptor sequences and removing low quality reads from raw data. Then, the clean reads were mapped to the rice reference genome based on MSU-v7.0 (http://rice.plant biolo gy.msu.edu/) using HISAT2 software . Next, the mapped reads were assembled using the StringTie software (Pertea et al. 2016). Then, the assembled transcripts were annotated using the gffcompare program. Finally, the protein-coding transcripts were identified and the remaining unknown transcripts were used to screen for putative lncRNAs.
For lncRNAs identification, the screening criteria were as follows: (1) the transcripts with a single exon and less than 200 bp long were removed; (2) the transcripts that were identified as known mRNAs and other small RNAs were removed; (3) the transcripts with protein-coding ability were removed based on the evaluation of coding potential calculator (CPC, CPC score > 0) (Kong et al., 2007), coding-non-coding index (CNCI, CNCI score > 0) (Sun et al., 2013), and coding potential assessing tool (CPAT) ); (4) the transcripts with known protein 1 3 Fig. 1 The senescence feature of rice flag leaves sampled at five developmental stages. a The sampled flag leaves at booting stage (FL1), flowering stage (FL2), early-senescence stage (FL3), midsenescence stage (FL4) and late-senescence stage (FL5), respectively. b The measured chlorophyll content of rice flag leaf. c, d The expression trends of two chlorophyll degradation-associated genes (SGR, LOC_Os09g36200;RCCR1, LOC_Os10g25030). e, f The expression trends of two senescence-associated genes (OsNAP, LOC_ Os03g21060;Osl43, LOC_Os01g24710). The data are expressed as the mean ± SD of three biological replicates. Different letters indicate values are statistically different based on one-way ANOVA analysis domains were also excluded according to protein family database (Pfam) (Finn et al. 2013). Finally, the remaining transcripts were considered as reliable lncRNAs. The different types of lncRNAs, including long intergenic noncoding RNAs (lincRNAs), antisense lncRNAs, intronic lncRNAs and sense lncRNAs were classified based on their location on the genome. Furthermore, the sequences of plant lncRNAs were downloaded from the Green Non-Coding Database (GreeNC) for sequence conservation analysis (Wang et al. 2019a).

Identification of differentially expressed (DE) lncRNAs and DE mRNAs
The fragments per kilobase of transcript per million mapped reads value calculated by StringTie software was used to quantify the expression levels of both lncRNAs and mRNAs in each sample. A pairwise differential expression analysis between any two stages (FL1 vs FL2;FL1 vs FL3;FL1 vs FL4;FL1 vs FL5;FL2 vs FL3;FL2 vs FL4;FL2 vs FL5;FL3 vs FL4;FL3 vs FL5;FL4 vs FL5) was performed using the DESeq R package (1.10.1). LncRNAs and protein-coding mRNAs with a false discovery rate < 0.05 and absolute value of log2 (Fold change) > 1 found by DE Seq were considered to be differentially expressed. Furthermore, hierarchical clustering analysis on DE lncRNAs and DE mRNAs were carried out using heatmap R package.

Target gene prediction of lncRNAs
Potential target genes of the lncRNAs were predicted based on their regulatory patterns, which were classified into cisand trans-acting groups (Stadler 2014). These protein coding genes transcribed within a 100 kb upstream or downstream of lncRNAs were searched as potential cis target genes by using Perl script (Tian et al. 2019). Furthermore, the trans-acting target genes were determined by calculating the expression correlation between the expression level of lncRNAs and mRNAs. And the correlation in expression was evaluated using Pearson's correlation coefficient (|r|> 0.9 and p < 0.01) (Tian et al. 2019).

Co-expression network analysis of DE lncRNAs and DE mRNAs
To further explore the functions of lncRNAs related to leaf senescence, a weighted gene co-expression network analysis (WGCNA) was conducted using the R package WGCNA (Langfelder and Horvath 2008). An unsigned co-expression relationship was built based on the adjacency matrix between DE lncRNAs and DE mRNAs. The one-step network construction and module detection were adopted using the "dynamic hybrid tree cut algorithm" with the power value of 5, minimum module size of 30 and merge cut height of 0.2826. The other parameters were defined as default values. Highly similar modules were subsequently identified by clustering and then merged into new modules on the basis of eigengenes. The correlation of each module was also analyzed and visualized by a heatmap. Finally, the coexpression network was visualized by Cytoscape software.

Function classification for the target mRNAs of interested lncRNAs
Gene Ontology (GO) enrichment analysis of the target mRNAs which co-expressed with interested lncRNAs were implemented using Blast2GO (Conesa et al. 2005), and GO terms with p-value < 0.05 were considered significantly enriched with differential expressed mRNAs. We used R package to test the statistical enrichment of the target mRNAs which co-expressed with interested lncRNAs in Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. The p-value < 0.05 was required for differences to be considered statistically significant. The KEGG pathways were enlisted in order according to the corrected p value, and pathways with a corrected p-value < 0.05 were considered to be most enriched.

Construction of lncRNAs-miRNAs-mRNAs network
To explore the coordinated functions of DE lncRNAs, the ceRNA networks were constructed. Firstly, known mature miRNA sequences for rice were downloaded from miRBase (release 22, October 2018). Subsequently, the identified DE lncRNAs and DE mRNAs in 'ME black' and 'ME blue' module were used as target prediction library. Then, psRNA-Target was used with default parameters to identify miRNA target binding sites (Dai et al. 2018). And regulatory data of miRNAs-DE lncRNAs and miRNAs-DE mRNAs were obtained. In addition, the target mRNAs were compared with the DE mRNAs in specific modules to obtain the crossover mRNAs. Finally, based on interaction relationship among DE lncRNAs-miRNAs-DE crossover mRNAs, ceRNA networks were constructed and visualized by Cytoscape software.

Quantitative real-time polymerase chain reaction (qRT-PCR)
To validate the relative expression level of lncRNAs and mRNAs, qRT-PCR analysis was performed. Firstly, the total RNAs of flag leaf samples at five developmental stages in rice were respectively extracted using TRIzol reagent (Invitrogen, USA). Next, random hexamers were used for cDNA synthesis of lncRNAs and mRNAs. Then, qRT-PCR was performed using RNA-specific primers with the SYBR Green PCR kit (TaKaRa, China) according to the manufacturer's instructions (Table S1). Three genes, UBC (LOC_Os02g42314), ARF (LOC_Os05g41060) and Profilin-2 (LOC_Os06g05880) (Wang et al. 2016), were used as internal reference genes to normalize the qRT-PCR data for mRNAs and lncRNAs. The relative RNA expression levels were calculated by using the 2 −ΔΔCT method (Ren et al. 2018). The reaction was carried out using three biological replicates with three technical replicates.

Analysis of senescence feature during flag leaf development
Chlorophyll degradation is a hallmark feature of leaf senescence. To effectively examine the leaf senescence status, we initially measured the chlorophyll content of rice flag leaves at five developmental stages, including booting stage (FL1, 9 days before flowering), flowering stage (FL2, 3 days after flowering), early-senescence stage (FL3, 9 days after flowering), mid-senescence stage (FL4, 19 days after flowering), and late-senescence stage (FL5, 29 days after flowering) (Fig. 1a). Accompanying with the leaf development from early-senescence stage (FL3) to late senescence stage (FL5), the chlorophyll content significantly decreased from 1392 to 553 μg/g FW (Fig. 1b). To confirm this tendency, the mRNA expression levels of two chlorophyll degradation associated genes SGR and RCCR1 were measured by qRT-PCR and found to increase during the senescence stages (Fig. 1c, d). Meanwhile, the senescence associated genes OsNAP and Osl43, were also higher expressed in the flag leaves of the senescence stage (Fig. 1e, f). These results suggested that there existed apparent leaf developmental differences among analyzed samples before and after senescence, implying that samples at these time points were suitable for subsequent study.

High throughput sequencing
To systematically identify lncRNAs and mRNAs related to leaf senescence, a genome-wide high throughput sequencing was performed for rice flag leaf at five developmental stages with the senescence feature identification described above. Sequencing was done on the Illumina HiSeq 2500 platform and 150 bp paired-end reads were generated. After filtering the adapters and low-quality reads, a total of 288.77 gigabases clean reads were obtained from 15 sequencing libraries (5 samples × 3 replicates) ( Table 1). And clean reads in each library mapped to the Oryza sativa reference genome ranged from 93.75 to 96.36% (Table 1). Furthermore, the Q30 of the clean reads was greater than 94.10% and GC content of the sequencing outputs were 45.67-46.95% (Table 1).

Identification of lncRNAs and mRNAs involved in flag leaf senescence of rice
After clean reads mapping to the Oryza sativa reference genome, the mapped transcripts were assembled and annotated using StringTie software (v1.3.1). And a total of 38757 mRNAs were identified (Table S2). The remaining  (Table S4). Accordingly, the result further indicated the low conservation of lncR-NAs among different organisms. The basic genomic features of these lncRNAs were additionally characterized. Compared with the mRNAs, the lncRNAs were mainly 400-1600 nt in length, and the number of lncRNAs decrease with the increasing transcript length (< 3000 nt) (Fig. 2c). The number of the corresponding exons of lncRNAs was much less than that of mRNAs and mainly below 5 (Fig. 2d), while that of mRNAs ranging from 1 to 30 (Fig. 2d). The length of open reading frames of lncRNAs was mainly below 100 nt ( Fig. 2e), while that of mRNAs was mainly below 900 nt (Fig. 2e). In addition, the expression level of mRNAs was higher than those of lncRNAs (Fig. 2f). These results thus revealed divergent features of lncRNAs compared with those of mRNAs.

Analysis of differentially expressed (DE) lncRNAs and DE mRNAs
To explore the transcriptional changes of lncRNAs and mRNAs at five stages, a total of ten pairwise comparisons groups were analyzed in this study. Based on the screening criteria of log2 (fold change) > 1 and false discovery rate < 0.  (Fig. 3a, c). Eventually, a total of 343 DE lncRNAs and 9412 DE mRNAs were obtained (Tables S5 and S6). In addition, the hierarchical cluster analysis of DE lncRNAs and DE mRNAs suggested that the existing differences of expression patterns at five stages during the process of flag leaf development to senescence in rice (Fig. 3b, d).

Analysis on potential target genes of lncRNAs
LncRNAs have been found to regulate the expression of protein-coding genes through cis-and trans-acting models . The protein-coding genes located within a genomic window of 100 kb upstream and downstream of lncRNAs were searched as potential cis-regulated target genes of lncRNAs (Tian et al. 2019). In this study, a total of 3631 lncRNAs were predicted to have potential cis-regulatory effects on 23273 protein-coding genes (Table S7). Among them, more than 80% lncRNAs targeted ten to fifty protein coding genes, and 22 lncRNAs have more than 40 target genes (Fig. 4a). More than 90% protein-coding genes corresponded to one to ten lncRNAs, and 514 protein-coding genes were predicted to be cis-regulated by more than ten lncRNAs (Fig. 4b). Moreover, lncRNAs may act in trans-, and the prediction of trans-regulation relies on the expression correlation between lncRNAs and mRNAs (Tian et al. 2019). In this study, a total of 2878 lncRNAs and 21307 associated target protein-coding genes were predicted to be trans-regulated (Table S7). Among these lncRNAs, 48% lncRNAs targeted more than 50 protein coding genes (Fig. 4c). Moreover, among the 21307 target protein-coding genes, 3873 (18%) corresponded to only one lncRNAs and up to 515 coding-genes were predicted to be targeted by more than 150 lncRNAs (Fig. 4d). These results may indicate the complex regulation relationships between lncRNAs and mRNAs.

Analysis of lncRNAs co-expressed with mRNAs by WGCNA
To systematically explore the potential regulation functions of lncRNAs associated with flag leaf senescence, WGCNA was performed to analyze the co-expression relationship between 343 DE lncRNA transcripts and 9412 DE mRNA transcripts. A total of seven distinct modules were obtained, in which major tree branches define the modules (labeled with different colors), as shown in the dendrogram (Fig. 5a). Furthermore, the modules closely related to flag leaf senescence were of particular interest to investigate by analyzing the module-trait correlations. Notably, of the seven modules, the 'ME black' module displayed a continuous down-regulation trends accompanying with the flag leaf development from normal to senescence, whereas the 'ME blue' module showed a continuous up-regulation trends during this process of leaf senescence (Fig. 5b). And the heatmaps of the 'ME black' module and the 'ME blue' module showed the comprised transcripts which were the The number of lncRNAs that have potential cis-regulatory effects on protein-coding genes. c The number of trans-target genes regulated by lncRNAs. d The number of lncRNAs that have potential trans-regulatory effects on protein-coding genes most highest expressed separately at booting stage (FL1) and at late senescence stage (FL5) (Fig. 6). Thus, considering the continuous down-and up-regulation trends during the process of flag leaf developing to senescence, the 'ME black' and 'ME blue' module were considered to be highly associated with leaf senescence (Fig. 5b).
The 'ME black' module comprised 22 DE lncRNA transcripts and 812 DE mRNA transcripts that displayed a continuous negative correlation with leaf senescence (Fig. 5b, Table S8). These 22 DE lncRNAs were predicted to target 2695 protein coding genes (Table S9), in which 378 target protein coding genes were DE mRNAs clustered in the 'ME black' module (Table S10). The 'ME blue' module comprised 48 DE lncRNA transcripts and 1209 DE mRNAs transcripts that displayed a continuous positive correlation with leaf senescence (Fig. 5b, Table S8). These 48 DE lncRNAs were predicted to target 3361 protein coding genes (Table S9), in which 941 target protein coding genes were DE mRNAs clustered in the 'ME blue' module (Table S10). The 22 DE lncRNAs from the 'ME black' module and 48 DE lncRNAs from the 'ME blue' module were selected as the most interested lncRNAs involving in flag leaf senescence.
Moreover, statistical enrichment of KEGG was analyzed for the co-expressed mRNAs targeted by interested DE lncRNAs involving in flag leaf senescence. Results showed that the co-expressed mRNAs targeted by down-regulated lncRNAs were significantly enriched in 11 KEGG pathways, in which 'photosynthesis-antenna proteins (ko00196)' and 'starch and sucrose metabolism (ko00500)' were the most enriched (Fig. 8a, Table S12). Meanwhile, the coexpressed mRNAs targeted by up-regulated lncRNAs were significantly enriched in 10 KEGG pathways, in which 'plant hormone signal transduction (ko04075)', 'alpha-Linolenic acid metabolism (ko00592)' and 'galactose metabolism (ko00052)' were the most enriched (Fig. 8b, Table S12).

Identification of interested lncRNAs targeting transcription factors (TFs)
It has been demonstrated that transcription factors (TFs) could play essential roles in senescence process (Podzimska-Sroka et al. 2015). In this study, a total of 378 target DE mRNAs targeted by 22 down-regulated DE lncRNA in the 'ME black' module were submitted to the Plant Transcription Factor Database (http://plant tfdb.cbi.pku.edu.cn/). Totally, 25 DE mRNAs were found to encode TFs, which were targeted by 10 DE lncRNAs in 'ME black' module.  (Table S13). For example, the gene BGIOSGA002217 encoding a C3H transcription factor was the target gene of lncRNA MSTRG.1416.1. Interestingly, one lncRNA could regulate multiple TFs, and one TF could also be regulated by multiple lncRNAs. For example, lncRNA MSTRG.56765.1 was predicted to regulate three TFs, containing BGI-OSGA001872 (MYB family), BGIOSGA027617 (LSD family) and BGIOSGA029710 (C2H2 family) (Table S13). And the NAC transcription factor BGIOSGA036980 was predicted to be regulated by three lncRNAs MSTRG.30264.1, MSTRG.24808.2 and MSTRG.65878.1 (Table S13). Similarly, a total of 941 DE mRNAs targeted by 48 up-regulated DE lncRNAs in the 'ME blue' module were also searched in the Plant Transcription Factor Database. And 86 DE mRNAs were found to encode TFs, which were predicted to be the target genes of 33 DE lncRNAs (Table S13). The The top 25 GO terms for 941 target DE mRNAs regulated by 48 upregulated DE lncRNAs in 'ME blue' module complex regulatory relationship between these DE lncRNAs and their target TFs in 'ME blue' module were also found. For example, lncRNA MSTRG.12677.1 was predicted to regulate 47 TFs, such as BGIOSGA000374 (NAC family), BGIOSGA013831 (C2H2 family) and BGIOSGA023457 (NAC family) (Table S13). And the NAC transcriptional factor BGIOSGA000374 was predicted to be regulated by 14 lncRNAs, such as MSTRG.12677.1, MSTRG.30528.1 and MSTRG.74673.3 (Table S13). These results suggested that senescence associated DE lncRNAs might play roles in leaf senescence through regulating transcriptional factors.

Identification of lncRNAs targeting well-studied senescence-associated genes (SAGs)
So far, a vast number of senescence-associated genes (SAGs) were identified to play roles in leaf senescence. It is interesting to investigate whether these well-studied SAGs were regulated by the identified senescence-associated lncRNAs in this study. Results showed that a total of eight down-regulated DE lncRNAs in the 'ME black' module were predicted to target 12 identified SAGs, of which three SAGs were in the 'ME black' module while the other nine were not (     Zhu et al. (2019) in the 'ME black' module (Table 2). Similarly, a total of 29 up-regulated DE lncRNAs in the 'ME blue' module were found to target 20 identified SAGs, in which 13 identified SAGs were found to locate in the same 'ME blue' module while the other seven SAGs were not (Table 2). For example, DE lncRNA MSTRG.8133.1 and its targeted mRNA BGIOSGA004591 (OsPME1) were both clustered in the 'ME blue' module, while BGIOSGA012281 (PLS2) targeted by the DE lncRNA MSTRG.20126.1, was not located in the 'ME blue' module. Interestingly, it was apparent to find that one identified SAG could be targeted by multiple lncRNAs. For example, the mRNA BGIOSGA023457 (ONAC011) was the target gene of 21 DE lncRNAs (Table 2). Meanwhile, one lncRNA could also be found to target multiple identified SAGs (Table S14). For instance, the DE lncRNA MSTRG.24808.2 in the 'ME black' module was found to target BGIOSGA004283 (SPL28), BGI-OSGA028186 (SPL29), BGIOSGA033278 (OsGDCH) and BGIOSGA014571 (cZOGT2), while MSTRG.77321.3 in the 'ME blue' module could target multiple SAGs such as BGIOSGA023457 (ONAC011), BGIOSGA011859 (OsPAO) and BGIOSGA011953 (OsFBK12) (Table S14). Eventually, the expression levels of down-expressed lncRNA MSTRG.24808.2 and its targeted mRNAs (BGI-OSGA033278 and BGIOSGA014571), and that of upexpressed lncRNA MSTRG.77321.3 and its targeted mRNAs (BGIOSGA023457 and BGIOSGA011859) were verified by qRT-PCR (Fig. 9). These results together implied that the senescence-associated lncRNAs played important roles in leaf senescence by regulating numerous senescenceassociated mRNA transcripts. And these results also supplied evidences for the involvement of identified lncRNAs in flag leaf senescence.

Analysis of lncRNA-associated ceRNA network in senescent leaf of rice
It has been reported that lncRNAs can act as miRNA sponge to form a ceRNA network and thus regulate the expression of target transcripts of the same miRNAs (Cesana et al. 2011). In this study, two ceRNA networks were separately constructed using senescence associated DE lncRNAs and DE mRNAs from the 'ME black' module and 'ME blue' module, combining with miRNAs from miRBase (release 22, October 2018). For down-regulated lncRNAs and their co-expressed mRNA in 'ME black' module, a total of five competitive relationship subnetworks were identified, including 6 lncR-NAs, 15 miRNAs and 51 mRNAs (Fig. 10a-e). Among these 51 mRNAs, none was identified as well-studied SAG. However, it was interesting to found that one key lncRNA MSTRG.44302.1 (purple triangle) may bind with osa-miR5809 to regulate 41 protein coding genes, in which six  Chen et al. (2013) mRNAs are TFs including BGIOSGA028751 (AP2 family), BGIOSGA005075 (SBP family), BGIOSGA022967 (bHLH family), BGIOSGA027807 (NAC family), BGI-OSGA034224 (ZF-HD family) and BGIOSGA023751 (COlike family) (Fig. 10a). Accordingly, a total of nine competitive relationship sub-networks were identified for up-regulated lncRNAs and their co-expressed mRNA in 'ME blue' module, including 14 lncRNAs, 21 miRNAs and 117 mRNAs (Fig. 10f-n). Interestingly, some TFs and well-studied SAGs were found in networks. For example, one key lncRNA MSTRG.62092.1 may combine with osa-miR164a and osa-miR164e to regulate 12 protein coding genes, in which five mRNAs are TFs including BGIOSGA037778 (NAC family), BGIOSGA008492 (NAC family), BGIOSGA016546 (OsNAC2, NAC family), BGIOSGA023457 (ONAC011, NAC family) and BGIOSGA002469 (MYB family) (Fig. 10f) Fig. 10 The lncRNAs-associated ceRNA networks formed by downand up-regulated lncRNAs and their regulated co-expressed mRNAs. a-e The ceRNA subnetwork formed by down-regulated lncRNAs and their regulated co-expressed mRNAs in 'ME black' module. f-n The ceRNA subnetwork formed by up-regulated lncRNAs and their regu-lated co-expressed mRNAs in 'ME blue' module. Triangle, diamond and circle indicated the lncRNAs, miRNAs and mRNAs, respectively. Purple nodes, green nodes and red circle nodes represented the emphasized competitive lncRNAs, miRNAs and known mRNAs (transcription factors and known well-studied SAGs) h BGI-OSGA016313. The data are expressed as the mean ± SD of three biological replicates. Different letters indicate values are statistically different based on one-way ANOVA analysis BGIOSGA016852 (HD-ZIP family), BGIOSGA026766 (bZIP family), and BGIOSGA014332 (NAC family) (Fig. 10g). Moreover, lncRNAs MSTRG.80924.2 and MSTRG.27317.2 may bind osa-miR156k to regulate the chlorophyll-degradation associated gene BGI-OSGA003046 (NYC1) (Fig. 10h).
To validate the regulatory relationship between lncRNAs and mRNAs, qRT-PCR experiments were performed. In the ceRNA network, five DE mRNAs and three DE lncRNAs were analyzed. The results showed that expression levels of lncRNA MSTRG.62092.1 and mRNAs BGIOSGA037778 (NAC family), BGIOSGA008492 (NAC family), and BGI-OSGA002469 (MYB family) increased with the flag leaf developing to senescence (Fig. 11a-d), suggesting that lncRNA MSTRG.62092.1 could bind with osa-miR164a and osa-miR164e to regulate the expression of these three TFs in the leaf senescence process. Similarly, lncRNAs MSTRG.31014.21 and MSTRG.31014.36, and mRNA BGI-OSGA025169 (OsNCED4) and BGIOSGA016313 (NAC family), also showed the increased expression levels during leaf senescence (Fig. 11e-h), implying that lncRNAs MSTRG.31014.21 and MSTRG.31014.36 could combine with miR5809 to regulate protein coding genes in this process. Thus, these key lncRNAs forming the ceRNA networks may play important functions in leaf senescence.

LncRNA is a key regulator during flag leaf senescence in rice
LncRNAs were well demonstrated to play vital roles in various biological processes including gene silencing, abiotic stress response, growth and development (Liu et al. 2015b). To date, a vast number of lncRNAs were identified in plants, such as melon (Tian et al. 2019), cassava , Prunus mume , Arabidopsis , rice (Yu et al. 2020) and cotton (Deng et al. 2018), revealing their important roles in multiple different biological processes. However, no systematic identification of lncRNAs involving in leaf senescence were reported. In this study, a total of 3953 lncRNAs were identified by genome-wide high-throughput sequencing during the developmental process of rice flag leaf from normal to senescence, in which 343 lncRNAs were found to be differentially expressed (Table S5). Among all ten pairwise comparisons of five analyzed stages, group 'FL1 vs FL3' had the highest number of DE lncRNAs (Fig. 3a), implying that the senescence status of flag leaves changed significantly from booting stage (FL1) to early senescence stage (FL3). On the contrast, the senescence status was more uniform in flag leaves from booting stage (FL1) to flowering stage (FL2), and mid-senescence stage (FL4) to late-senescence stage (FL5). Furthermore, the continuous positive and negative modules ('ME black' and 'ME blue') correlated with leaf senescence were found by WGCNA analysis (Fig. 5), in which, DE lncRNAs targeted their coexpressed mRNAs participating in senescence-associated biological processes and pathways (Figs. 7 and 8). And these interested lncRNAs were also predicted to target transcription factors and well-studied SAGs (Tables 2 and S13). We speculated that lncRNAs in these two modules could most possibly play important roles in leaf senescence. In addition, the ceRNA network analysis revealed that some senescence-associated lncRNAs were predicted to regulate their co-expressed mRNAs through the miRNA sponges (Figs. 10 and 11). These results provided a molecular basis for understanding about regulatory functions of lncRNAs on flag leaf senescence, and laid a foundation for further functional studies on candidate lncRNAs.

Involvement of lncRNAs in flag leaf senescence by targeting transcription factors
Transcription factors are important regulatory elements for the transcriptional activation or repression of target genes by recognizing their specific regulatory DNA sequences usually in the promotor region. Studies on the regulatory mechanisms and functions of transcription factors have increasingly been considered as research hotspots for leaf senescence (Podzimska-Sroka et al. 2015). In this study, whether senescence associated lncRNAs might regulate transcription factors were investigated. GO analysis revealed 17 and 69 mRNAs were separately enriched for down-regulated and up-regulated lncRNAs in the term of 'regulation of transcription, DNA-templated (GO:0006355)' (Fig. 7, Table S11). Additional analysis on transcription factors identified 25 and 86 transcription factors, which were separately targeted by 10 down-regulated lncRNAs and 33 up-expressed lncRNAs (Table S13). A complex regulation network might exist among these lncR-NAs and TFs in the process of leaf senescence. Functional analysis of some specific TFs was reported for leaf senescence in rice. OsDOS is a novel nuclear-localized C3Htype transcription factor, and RNAi knockdown of OsDOS showed accelerated leaf senescence, shorter panicle and slightly reduced seed setting whereas its overexpression appeared a delayed leaf senescence and severe sterility (Kong et al. 2006). Gene BGIOSGA002217 was identified to be OsDOS, and it was predicted to be the target gene of lncRNA MSTRG.1416.1 (Table 2). ONAC011 is a NAC-type transcription factor. Transgenic rice overexpressing ONAC011 accelerated heading time and promoted leaf senescence, while blocking the function of ONAC011 Fig. 12 The network models of senescence-associated lncRNAs involved in flag leaf senescence. a A proposed regulatory mechanism involving lncRNAs and their potential and well-studied target genes during flag leaf senescence. The heatmaps represent the expression patterns of lncRNAs and their potential target genes in five analyzed leaf stages. LncRNAs participating in more than four different biolog-ical processes were labeled in different colors, and recommended to be prior candidates for transgenic analysis. b The model of lncRNAs by regulating well-studied mRNAs through ceRNA network to participate in senescence-related biological processes. LncRNAs in the dotted box probably regulating multiple senescence-related mRNAs, were recommended to be prior candidates for transgenic analysis gene via RNA interference could delay both heading time and leaf senescence. Irrespective of early or delayed senescence, transgenic plants showed reduced grain yields (El Mannai et al. 2017). BGIOSGA023457 was identified to be ONAC011, and surprisingly, up to 21 DE lncRNAs were predicted to target ONAC011 (Table 2). As a result, it is speculated that lncRNAs may play key roles in leaf senescence by transcription regulation (Fig. 12a).

Involvement of lncRNAs in flag leaf senescence by targeting genes in hormone pathways
Leaf senescence is a developmental process controlled by complex networks of interacting genes and signaling pathways, and plant hormones play vital roles in regulating this process (Lim et al. 2007). In this study, during the process of investigating the functions of interested lncRNAs, a lot of hormone related GO terms were found, such as 'response to hormone (GO:0009725)', 'hormone-mediated signaling pathway (GO:0009755)', 'hormone metabolic process (GO:0042445)' (Fig. 7, Table S11). It is known to all that jasmonates act as a signaling in numerous developmental processes including senescence, which are formed from alpha-linolenic acid of chloroplast membranes by oxidative processes (Wasternack and Song 2016). Phenylpropanoid was reported to involve in the formation of salicylic acid, which played roles in leaf senescence (Strawn et al. 2007). And tryptophan was involved in the biosynthesis of melatonin, which could delay senescence and protect photosynthetic systems (Arnao and Hernández-Ruiz 2018). KEGG analysis revealed that 'alpha-linolenic acid metabolism (ko00592)', 'phenylpropanoid biosynthesis (ko00940)' and 'tryptophan metabolism (ko00380)' were all enriched in this study, and 'plant hormone signal transduction (ko04075)' was the most significantly enriched (Fig. 8, Table S12).
It is also founded that interested lncRNAs may target identified hormone-related genes involving in leaf senescence. In rice, mutation of COI homologs OsCOI1b, mainly affects leaf senescence by mediating jasmonate signaling, and further reduces grain production by decreasing spikelet fertility and grain weight (Lee et al. 2015). In this study, BGIOSGA019953 (OsCOI1b) was predicted to be the target gene of 7 lncRNAs (Table 2). OsPME1 is an epigenetically regulated gene, affecting the content of jasmonates and regulating leaf senescence. Accelerated leaf senescence and chlorophyll degradation were observed in plants overexpressing OsPME1, while OsPME1 RNAi plants displayed retarded leaf senescence (Fang et al. 2016). The gene BGIOSGA004591 (OsPME1) in 'ME blue' module was predicted to be the target gene of lncRNA MSTRG.8133.1 (Table 2). OsNCED3 is a key gene in regulating abscisic acid biosynthesis, and overexpression of OsNCED3 in transgenic rice could promote leaf senescence and increase abscisic acid content, whereas its knockout using CRISPR/ Cas9 could delay leaf senescence (Huang et al. 2018). Gene BGIOSGA013214 in 'ME blue' module was identified to be OsNCED3, and this gene was predicted to be the target gene of 12 lncRNAs (Table 2). OsRTH1 play roles in the modulation of ethylene responses, and its overexpression substantially prevented ethylene-induced alterations in growth and development, including leaf senescence (Zhang et al. 2012). Gene BGIOSGA004316 in 'ME blue' module was identified to be OsRTH1, and this gene was found to be the target gene of lncRNAs MSTRG.4640.1 and MSTRG.4077.5 (Table 2). Cytokinins was reported to involve in the regulation of various biological processes including senescence (Kim et al. 2006). The putative cis-Zeatin-O-glucosyltransferase gene cZOGT2 plays roles in regulating cytokinin activity of cis-Zeatin, and transgenic rice lines ectopically overexpressing the cZOGT2 gene exhibited delay of leaf senescence (Kudo et al. 2012). Interestingly, gene BGIOSGA014571 was found to be cZOGT2, and this gene was predicted to be the target gene of lncRNAs MSTRG.24808.2 and MSTRG.65878.1 (Table 2). Exogenous application of melatonin to rice plant indicated that melatonin acts as a biostimulator to delay leaf senescence (Liang et al. 2015). Serotonin N-acetyltransferase (SNAT) is a key gene for melatonin biosynthesis, and overexpression of OsSNAT1 in transgenic rice plants confers resistance to cadmium and senescence, and increases grain yield (Lee and Back 2017). BGIOSGA020059 (OsSNAT1) was found to be the target gene of lncRNA MSTRG.56765.1 in this study (Table 2). Based on all above evidences and results, it is proposed that lncRNAs may play important roles in flag leaf senescence by involving in hormone pathways (Fig. 12a).

Involvement of lncRNAs in flag leaf senescence by targeting genes related to redox metabolism
The redox state of cells is important for plant. Reactive oxygen species (ROS), which are continually produced by aerobic metabolism, are thought as cytotoxic molecules when in high levels, while low levels of ROS can act as signaling molecules to regulate plant development, stress responses, programmed cell death, and senescence (Jajic et al. 2015). In present study, interested senescence-associated lncRNAs were found to target a large number of mRNAs enriched in ROS-related GO terms, such as 'oxidation-reduction process (GO:0055114)', 'response to oxygen-containing compound (GO:1901700)' and 'cellular response to oxygencontaining compound (GO:1901701)' (Fig. 7, Table S11). Ascorbate peroxidase (APX) enzyme plays an essential role in the control of intracellular H 2 O 2 levels. And OsAPX4 gene has been reported to play an important role in leaf senescence mediated by ROS signaling (Ribeiro et al. 2017). In this study, BGIOSGA026526 (OsAPX4) was predicted to be the target gene of lncRNA MSTRG.79278.1 (Table 2). NAD plays critical roles in cellular redox reactions and remains at a sufficient level in the cell to prevent cell death. OsNaPRT1, encoding the nicotinate phosphoribosyltransferase in the NAD salvage pathway, protects rice plant from premature leaf senescence by maintaining the redox balance (Wu et al. 2016). BGIOSGA013898 (OsNaPRT1) was predicted to be the target gene of lncRNA MSTRG.44302.1 (Table 2). It is suggested that senescence associated lncRNAs perform functions probably by regulating genes involved in the redox metabolism (Fig. 12a).

Involvement of lncRNAs in flag leaf senescence by targeting lesion mimic genes
Early leaf senescence is usually present as one of the characteristics for lesion mimic mutants. Functional analysis of the genes responsible for lesion mimic mutants have provided a unique tool for dissecting possible molecular mechanism and regulatory pathway during leaf senescence. The rice lesion-mimic gene SPL4 encodes a plant spastin that inhibits ROS accumulation in leaf development and functions in leaf senescence, and its mutant showed low seed setting rate and panicle number, but increased 500-grain weight ). SPL28 appears to be involved in the regulation of vesicular trafficking, and SPL28 dysfunction causes the formation of hypersensitive response-like lesions, leading to the initiation of leaf senescence, small panicles, produced wizened grains, and few viable seeds (Qiao et al. 2010). SPL29 encodes UDP-N-acetylglucosamine pyrophosphorylase 1, and its functional inactivation induces early leaf senescence and defense responses in rice, accompanying with decreased panicle length, grain number, and thousand-grain weight . The OsPLS1 gene encoding vacuolar H + -ATPase subunit A1 (VHAA1), is implicated in leaf senescence through a combination of ROS and SA signals . In present study, these lesion mimic genes BGIOSGA022234 (SPL4), BGI-OSGA004283 (SPL28), BGIOSGA028186 (SPL29), and BGIOSGA020739 (OsPLS1) were predicted to be the target gene of one or two lncRNAs (Table 2). It is proposed to be an important way that lncRNAs participates in leaf senescence though regulating lesion mimic genes (Fig. 12a).

Involvement of lncRNAs in flag leaf senescence by targeting genes in other biological processes
Calcium-dependent protein kinases as a type of calcium sensor, perform multiple biological function in plants including senescence and cell death. OsCPK12 encodes a calciumdependent protein kinase. The mutant of OsCPK12 triggers the premature leaf senescence, while the overexpression of OsCPK12 may delay its growth period and provide the potentially positive effect on productivity in rice (Wang et al. 2019b). In our study, BGIOSGA014554 (OsCPK12) is predicted to be the target gene of lncRNA MSTRG.51442.2 (Table 2). Protein glycosylation may paly critical roles in flag leaf senescence (Huang et al. 2019a). Glycosyltransferases catalyzing the transfer of sugar moieties for glycosylation. PLS2, a putative glycosyltransferase-encoding gene, is involved in leaf senescence (Wang et al. 2018b). In present investigation, BGIOSGA012281 (PLS2) is identified to be the target gene of lncRNA MSTRG.20126.1 (Table 2). Senescence could be induced by abiotic stresses. A salt-induced protein gene salT was found to affect the leaf senescence induced by natural and dark conditions , and BGIOSGA001700 (salT) was predicted to be the target gene of 10 lncRNAs (Table 2). OsGDCH encodes the H-protein subunit of the glycine decarboxylase complex, playing a major role in photorespiration, and knockdown of GDCH gene reveals reactive oxygen species-induced leaf senescence in rice (Zhou et al. 2013). BGIOSGA033278 (OsGDCH) identified in this study was predicted to be the target gene of three lncRNAs MSTRG.44376.1, MSTRG.24808.2 and MSTRG.65878.1 (Table 2). Transgenic rice overexpressing an F-box protein containing a Kelch repeat motif OsFBK12 delayed in leaf senescence and increased seed size, while knockdown lines could promote the senescence program ). BGIOSGA011953 (OsFBK12) identified in this study was predicted to be the target gene of 15 lncRNAs (Table 2). It is seemed that lncRNA can be involved in flag leaf senescence by participating in several different biological processes (Fig. 12a).

Vital roles of key lncRNAs in lncRNAs-miRNAs-mRNAs network
It has been reported that lncRNAs can competitively bind to miRNAs through miRNA response elements and affect the expression of genes (Salmena et al. 2011). In the present study, based on the interested down-regulated and up-regulated lncRNAs, and their co-expressed mRNAs, two ceRNA networks were constructed, from which, some key lncRNAs may regulate a large number of mRNAs through miRNAs to participate in flag leaf senescence (Fig. 10).
The roles of miR164 and its NAC target genes were reported to involve in regulating a numerous of biological processes. It has been reported that miR164 is a conserved microRNA with miR164-targeted NAC genes negatively regulating drought resistance in rice (Fang et al. 2014). In Kiwifruit, a genome-wide analysis of coding and non-coding RNA revealed a conserved miR164-NAC regulatory pathway playing key roles in fruit ripening (Wang et al. 2019c). A conserved role for the NAM/miR164 developmental module revealed a common mechanism underlying carpel margin fusion in Monocarpous and Syncarpous Eurosids (Vialette-Guiraud et al. 2016). In addition, the rice transcription factor ONAC011 could promote leaf senescence and accelerate heading time in rice (El Mannai et al. 2017). The OsNAC2 was found to promote leaf senescence via abscisic acid synthesis (Mao et al. 2017). In ceRNA network constructed by up-regulated lncRNAs and their coexpressed mRNAs, lncRNA MSTRG.62092.1 was predicted to regulate four NAC TFs, BGIOSGA023457 (ONAC011), BGIOSGA016546 (OsNAC2), BGIOSGA037778, BGI-OSGA008492, acting as the target genes of osa-miR164a and osa-miR164e (Fig. 10f). These indicated that a conserved miR164-NAC regulatory pathway may be presented in the process of flag leaf senescence. Phytohormones including jasmonates and abscisic acid, were reported to play key roles in leaf senescence. Epigenetically regulation ofOsPME1 was reported to affect the content of jasmonates and thus play an important role in regulating leaf senescence (Fang et al. 2016). In abscisic acid biosynthetic pathway, 9-cis-epoxycarotenoid dioxygenase (NCED) is the key ratelimiting enzyme. The CRISPR/Cas9 knockout of abscisic acid-biosynthesis genes OsNCED3 was identified to delay leaf senescence in transgenic rice, while its overexpression of OsNCED3 could promote leaf senescence (Huang et al. 2018). And OsNCED5 overexpression increased abscisic acid level, enhanced tolerance to the stresses and accelerated leaf senescence (Huang et al. 2019b). Moreover, the TF ONAC066 can function as a positive regulator to resist drought induced plant death (Yuan et al. 2019). In this study, mRNAs BGIOSGA004591, BGIOSGA025169, and BGIOSGA013656 were identified as OsPME1, an abscisic acid-biosynthesis gene OsNCED4, and ONAC066, respectively. And lncRNAs MSTRG.31014.21, MSTRG.31014.36 and MSTRG.39310.1 were proposed to regulate these three well-studied SAGs and another ten unreported TFs through osa-miR5809 (Fig. 10g).
Taken together, these ceRNA network results provided several key lncRNAs involving in flag leaf senescence (Fig. 12b), and might lay a foundation for investigating the important function of candidate lncRNAs in flag leaf senescence through transgenic technology.

Conclusion
Overall, in this study, a total of 3953 lncRNAs and 38757 mRNAs were identified in flag leaves of rice at five developmental stages from normal to senescence by using the genome-wide high throughput sequencing. And 343 lncR-NAs and 9412 mRNAs were differentially expressed in ten pairwise comparison groups. LncRNAs showing continuous negative/positive correlation with flag leaf senescence were additionally identified by WGCNA. GO enrichment results suggested that the senescence negatively and positively associated DE lncRNAs targeted their co-expressed DE mRNAs involving in many biological processes including transcription, hormone response, oxidation-reduction process and substance metabolism.