Integrated transcriptomics and miRNAomics provide insights into the complex multi-tiered regulatory networks associated with coleoptile senescence in rice

Coleoptile is the small conical, short-lived, sheath-like organ that safeguards the first leaf and shoot apex in cereals. It is also the first leaf-like organ to senesce that provides nutrition to the developing shoot and is, therefore, believed to play a crucial role in seedling establishment in rice and other grasses. Though histochemical studies have helped in understanding the pattern of cell death in senescing rice coleoptiles, genome-wide expression changes during coleoptile senescence have not yet been explored. With an aim to investigate the gene regulation underlying the coleoptile senescence (CS), we performed a combinatorial whole genome expression analysis by sequencing transcriptome and miRNAome of senescing coleoptiles. Transcriptome analysis revealed extensive reprogramming of 3439 genes belonging to several categories, the most prominent of which encoded for transporters, transcription factors (TFs), signaling components, cell wall organization enzymes, redox homeostasis, stress response and hormone metabolism. Small RNA sequencing identified 41 known and 21 novel miRNAs that were differentially expressed during CS. Comparison of gene expression and miRNA profiles generated for CS with publicly available leaf senescence (LS) datasets revealed that the two aging programs are remarkably distinct at molecular level in rice. Integration of expression data of transcriptome and miRNAome identified high confidence 140 miRNA-mRNA pairs forming 42 modules, thereby demonstrating multi-tiered regulation of CS. The present study has generated a comprehensive resource of the molecular networks that enrich our understanding of the fundamental pathways regulating coleoptile senescence in rice.


Introduction
Plant senescence, the terminal phase of plant development, is a regimented degradation process which ensures efficient nutrient recycling, fitness and reproductive success. It is a natural phenomenon primarily governed by developmental age and has a multi-layered regulation in which phytohormones play a pivotal role in integrating developmental signals with extrinsic factors, thereby fine-tuning the senescence program by an intricate network of signaling pathways (Lim et al., 2007). Abscisic acid (ABA), jasmonic acid (JA), salicylic acid (SA), ethylene, brassinosteroids (BRs) and strigolactones are wellacknowledged positive regulators of plant senescence, while cytokinins, gibberellins and auxins are generally known to negatively regulate plant senescence. Eventually all the cells and tissues undergo slow death in the plant organs that are destined for senescence. Among different plant organs, extensive studies have been carried out for understanding the molecular mechanism underlying the regulation of leaf senescence in plants. The senescence-associated deterioration process is tightly controlled at transcriptional, post-transcriptional and epigenetic level. The transcriptional regulation of senescence refers to the changes in the expression pattern of genes which are known as senescence-associated genes (SAGs). Several transcriptome studies have revealed massive reprogramming of a large repertoire of genes predominantly belonging to families such as transcription factors (TFs), proteases, receptor-like kinases, and nucleases. Members of TF families such as NAC, WRKY, bHLH, Zinc finger, MYB, TCP and AP2 participate in the senescence programme and many of these represent nodal points for mediating crosstalk between external and internal factors to initiate the aging process. Like many other processes of plant development, the last stage is also regulated at posttranscriptional level by small non-coding RNAs (Thatcher et al., 2015). High-throughput sequencing has fuelled the discovery of a large number of miRNAs associated with leaf senescence in multiple plant species Huo et al., 2015;Wu et al., 2016). Three miRNA regulatory modules and their cognate TF targets i.e., miR164-ORE1, miR319-TCP4 and miR390-TAS3-ARF2 play critical role in regulating leaf senescence in Arabidopsis (Ellis et al., 2005;Kim et al., 2009;Li et al., 2013). Senescence is also regulated epigenetically by DNA methylation, histone modifications and chromatin alterations (Ay et al., 2014). Overall, the complex regulatory interactive network comprising SAGs (mainly TFs), miRNAs, hormones and chromatin modifications modulate leaf senescence in plants.
Among different plant organs, coleoptile is the first organ to undergo senescence post-germination in cereal plants. It is a protective sheath that emerges two days after sowing and covers emerging shoots in monocotyledonous plants. In addition to providing physical protection to first true leaves, the coleoptile also provides nutrition, mainly carbohydrates, to the growing parts of developing seedlings (Fröhlich and Kutschera, 1995). During seed germination under submerged conditions, coleoptile rapidly elongates and acts like a 'snorkel' to uptake oxygen from the water surface and transfer it to the underwater parts for their optimal growth (Narsai et al., 2015). It is evidently clear that the coleoptile, along with the first leaf, play a crucial role in seedling establishment, which is a key determinant of plant fitness and productivity of the cereal crops. Interestingly, the coleoptile continues to elongate under anaerobic conditions, however, on transfer to aerobic conditions it rapidly attains maturation followed by growth cessation and initiation of senescence (Inada et al., 2002). Under optimum growth condition, the coleoptile turns green from white within 3 days of seedling transfer to air. This is followed by emergence of the first leaf along with splitting of coleoptile in the next 24 h. Senescence of coleoptile is initiated rapidly afterwards and is evident by the formation of lysigenous aerenchyma, chlorophyll degradation and DNA fragmentation (Inada et al., 2002). Other than coleoptile, leaf senescence also occurs naturally in rice plants. Though deterioration events such as chloroplast breakdown followed by nuclear condensation are conserved in coleoptile and leaf senescence, aging within an individual leaf is highly asynchronous, thereby making it difficult to accurately characterize its senescence program and associated regulatory pathways (Inada et al., 1998). On the contrary, the senescence in coleoptile is relatively synchronous because of its simple cellular architecture consisting of only two vascular bundles and few layers of mesophyll cells. Owing to the simple cellular architecture, short life-span and quick and uniform pattern of aging, coleoptile is often proposed as a model organ for studying the mechanism of senescence in cereals.
Comprehensive understanding of the senescence process in coleoptile has been gained with respect to the pattern of cell death events and major changes that occur at tissue and cellular level, however, till date there is no genomic-scale study on discovering the gene expression changes during coleoptile senescence in any cereal plant. With an aim to gain insights into the components of regulatory pathway(s) and possible interaction among these components/pathways that contribute to coleoptile senescence, we employed an integrated omics approach and performed high-throughput sequencing of coding as well as non-coding RNAs to identify the differentially expressed genes and their regulation by miRNAs. A comparative analysis of the coleoptile transcriptome dataset with the publicly available leaf senescence dataset revealed the unique and common components between the two senescing organs. Efforts were made to validate some of the transcripts and miRNAs that exhibited differential expression during senescence. Regulatory networks involving differentially expressed genes and their probable miRNA regulators were integrated with protein-protein interaction (PPI) datasets, which vividly supported a multi-tiered regulation of senescence program in rice coleoptiles.

Plant material
The surface-sterilized rice seeds (var. Pusa Basmati 1) were submerged in autoclaved water for seven days at 28 0 C in test tubes (8-10 cm below the water surface) placed in a growth room. After seven days of growth under submergence, the seeds with elongated coleoptiles were transferred to cotton beds soaked with rice growth media (Yoshida et al., 1976). Immediately after transferring to aerobic conditions, the unsenesced elongated coleoptile tissues (control-Day 0) were harvested. Senescence of coleoptiles was monitored on Day 1 and Day 2 under aerobic conditions, and tissues were harvested at each time point. Tissues from three biological replicates were harvested, frozen in liquid nitrogen and stored in -80°C for future use. Out of the three replicates, two replicates of control and Day 2 tissues were employed for high-throughput sequencing of RNA and small RNA. For quantitative PCR experiments three biological replicates comprising all three timepoints (control, Day 1 and Day 2) were utilized. mRNA and small RNA library preparation and high-throughput sequencing Total RNA was isolated from the harvested tissues using a modified GITC method (Chomczynski and Sacchi, 1987). The quantity and quality of the isolated RNA were evaluated using a spectrophotometer (Bio-Rad, USA) and MOPS-formaldehydeagarose denaturing gel electrophoresis. Four RNA-Seq libraries were constructed from control (Control_BR1 and Control_BR2) and Day 2 (Day 2_BR1 and Day 2_BR2) samples followed by their 101 bp paired-end sequencing using Illumina's HiSeq4000 sequencing platform (Illumina, Inc., USA). For mRNA library preparation, RNA integrity of the control and Day 2 samples were assessed using Agilent Bioanalyzer 2100 system (Agilent Technologies, USA). TruSeq RNA Sample Preparation Kit v2 (Illumina, Inc., USA) was used according to the manufacturer's instructions to construct the RNA sequencing libraries. The transcriptome library preparation and sequencing were outsourced to M/s Bionivid Technology Pvt. Ltd., India.
For small RNA library preparation, total RNA was resolved on TBE-urea-PAGE gel and the 15-30 nt portion corresponding to the low molecular weight (LMW) RNA fraction was excised and extracted for LMW RNA. The extracted LMW RNA was ligated with RNA adapters (5' and 3') using T4 RNA ligase and subsequently used for first-strand cDNA synthesis using Superscript III reverse transcriptase (Invitrogen, USA). A 3' adapter specific reverse primer was used for the first-strand cDNA synthesis (Supplementary Table 11). The cDNA library was enriched and purified before quality check using Bioanalyzer (Agilent Technologies, USA). The small RNA-cDNA libraries were sequenced using Illumina HiSeq4000 sequencing platform by M/s Bionivid Technology Pvt. Ltd., India. The transcriptome and small RNA sequencing datasets were submitted to the SRA database as accession nos. GSE199139 and GSE199238, respectively.

Differential gene expression analysis, gene ontology and functional enrichment analysis
The sequencing data obtained by RNA sequencing was analyzed for quality assessment. In brief, filtering of high quality (HQ) reads was performed with 'NGS QC Toolkit (Patel and Jain, 2012). The HQ reads with an average Phred score ≥20 and ≥70% high-quality bases were included in the downstream analysis. Oryza sativa-Nipponbare-Reference-IRGSP-1.0 was used as the reference genome for mapping HQ adapter-trimmed reads. TopHat v2.0.11 was employed for mapping with default parameters (Trapnell et al., 2009). Subsequently, Cufflink (Cufflink_2.2.1) and Cuffmerge (Trapnell et al., 2009) tools were used for reference-based assembly of the reads. Fragments per transcript kilobase per million fragments mapped (FPKM) for each transcript was calculated from the number of mapped reads of the transcripts. Transcripts with normalized FPKM value ≥1 were considered for differential expression analysis with the Cuffdiff module of Cufflink. A comparative analysis of FPKM values of all the transcripts from Day 2 with that of control tissue was performed for identifying differentially expressed genes (DEGs). Identification of DEGs in coleoptile senescence was performed by calculating log (base 2) of FPKM of senescence sample divided by FPKM of control condition sample. A fold change of ±2 in the log2 scale was set as a cutoff for the categorization of DEGs. For the transcripts which were either specifically expressed in control or Day 2, fold change could not be calculated owing to FPKM value 0 in one of the samples. Since there were two biological replicates p-value was calculated to identify the transcripts which exhibited significant differential expression during senescence in the two biological replicates and transcripts with p-value <0.05 were shortlisted as a sub-dataset for downstream analyses.

Small RNA analysis and target prediction
The UEA Small RNA Workbench (version 2.4) was used to analyze purity-filtered raw reads (Stocks et al., 2012). The adapter sequences were removed from the raw reads and the resulting reads were converted into FASTA format by using the Adapter removal tool. Further, the reads were filtered according to the length and complexity. The invalid sequences, tRNAs, rRNAs and other non-coding RNA sequences were removed by filtering tools. For predicting mature miRNAs (novel and conserved), the miRCat tool was employed. Further, the Reference genome (Oryza sativa -Nipponbare-Reference-IRGSP-1.0) was used for mapping. miRBase v21 was used to identify the conserved miRNAs wherein two mismatches in the miRNA sequence were allowed. After filtering conserved miRNAs, novel miRNA predictions were carried out in the remaining small RNA population. Furthermore, the putative miRNAs were also screened manually according to criteria proposed by Meyers et al. (2008) and Kozomara and Griffiths-Jones (2014). To identify differentially expressed miRNAs, read counts were normalized and fold change was calculated as described by Xin et al. (2010). For target prediction, Plant small RNA target server (psRNATarget; https://www.zhaolab.org/psRNATarget/) was employed with default parameters and the expectation value was set to 4. MSU rice genome annotation, version 7 was used as the reference genome dataset. To generate a regulatory network of miRNA and target genes, we identified the predicted target IDs in the coleoptile transcriptome dataset generated using the VLOOKUP function in Microsoft excel.

Gene expression profiling and validation of miRNAs
For the validation of transcriptome data, SYBR green based-qPCR method was employed, and the rice eEF1a was employed as an internal control for normalization. For validation and evaluation of expression of predicted rice miRNAs, poly A tailing method of qPCR method was employed with SYBR chemistry. 2 µg total RNA was poly-adenylated by using poly A polymerase enzyme (Ambion, USA) followed by a purification step using the RNAeasy MinElute Cleanup Kit (Qiagen, GmBH). The RTQ primer (Supplementary Table 11) along with Superscript III reverse transcriptase (Invitrogen, USA) was used to reverse transcribe poly A-tailed RNA population. The universal reverse primer and miRNA-specific forward primer were used for Real-time PCR amplification in Mastercycler Realplex 2 (Eppendorf, Germany) with KAPA PROBE FAST Universal qPCR kit (KAPA Biosystems, USA). For the normalization, rice 5S ribosomal RNA was used as an endogenous control. Relative expression levels were calculated after normalization against the control sample. Fold change was calculated using 2 -DDct method (Livak and Schmittgen, 2001) and an average of all the replicates was plotted along with the calculated standard error (SE).

Protein-protein interaction network creation and visualization
The protein-protein interaction of the high confidence putative miRNA target genes which exhibited inverse expression patterns in the transcriptome data with respect to the corresponding miRNA profile were analyzed using the STRING 11 database. The output obtained was integrated with the miRNAs and their predicted targets to generate the miRNAtarget-PPI network using Cytoscape_v3.7.2 software.

Characterization of coleoptile senescence at morphological and molecular level
With an aim to obtain elongated coleoptiles with synchronized growth and senescence, rice seeds were grown under submerged conditions for seven days. To induce rapid maturation and initiate senescence in coleoptiles, submerged seedlings were transferred to aerobic conditions (taken as 'control'), and progression of senescence was subsequently monitored on day 1, day 2 and day 3 as depicted in Figure 1A. Closer examination of morphological alterations during coleoptile senescence revealed that wilting, curling and yellowing/browning of the coleoptile, especially at tips, was discernible at day 2 stage. On day 3, while the first leaf emerged, there was extensive wilting and browning of coleoptiles indicating late stage of senescence. Physiomorphological changes during progression of coleoptile senescence were confirmed at the molecular level by determining relative gene expression levels of senescenceassociated marker genes (SAGs): RED CHLOROPHYLL CATABOLITE REDUCTASE (RCC reductase) and WRKY88. Increased expression of RCC reductase gene encoding for RCC reductase enzyme (which catabolizes chlorophyll) is a hallmark of senescence in leaves and other plant organs (Tang et al., 2011). Similarly, increased expression of WRKY88 is reported during flag leaf senescence in rice (Lee et al., 2017). Quantitative PCR analysis demonstrated significant upregulation in steady state mRNA levels of both the genes during the progression of coleoptile senescence. The level of expression dramatically increased on day 1 and day 2. Notably, RCC reductase exhibited higher levels of induction at all the stages as compared to that of WRKY88 ( Figure 1B). These findings ratify the staging of coleoptile senescence in the present study.
In view of these observations, day 2 was considered as the middle stage of coleoptile senescence which could be utilized for studying the changes occurring at molecular level by transcriptome and small RNAome analysis.  Table 1). Scatter plots plotted for comparing gene expression between the two independent biological replicates of control and day 2 samples showed good correlation (Supplementary Figures 1A, B). As evident from the volcano plot, the number of significantly upregulated transcripts (1970) were higher compared to the downregulated ones (1195) (Figures 2A, B). Overall, 3166 transcripts were differentially expressed (DEGs) in senesced coleoptile samples when compared to the control, while 44 and 273 genes exhibited specific expression in control and senesced coleoptiles, respectively ( Figure 2B; Supplementary Table 3). These observations indicate that a large number of genes are regulated, during progression of senescence in rice coleoptiles. The genes exclusively regulated during senescence at Day 2 were also included under DEGs for further downstream analyses. To predict the function of genes exhibiting altered expression, GO analysis was carried out and GO terms could be assigned to 1792 out of 3439 DEGs. Maximum number of genes were associated with molecular function (77%) followed by cellular components (71%) and biological processes (58%) ( Figure 2C, Supplementary Table 3). The top terms enriched under molecular function, cellular component and biological process GO categories were "integral component of ATP binding", "integral component of membrane", and "transcription" respectively ( Figure 2D).

Functional category enrichment analysis of differentially expressed genes during coleoptile senescence
To predict the functional categories associated with the DEGs, enrichment analysis was performed using the MapMan annotation tool (Thimm et al., 2004). Out of the 3439 DEGs, 438 and 823 genes mapped to the metabolism and regulation pathways, respectively ( Figures 3A, B). Genes belonging to different functional categories such as transcriptional regulators (RNA category), transporters, hormone metabolism, redox related, stress response and signalling pathways are already known to play a role in the senescence of leaves (Guo et al., 2004;Gregersen and Holm, 2007;Zhang et al., 2014;Lin et al., 2015) and our analyses showed that 176, 37, 183, 134, 143, 94 and 105 transcripts belonging to signaling, redox, transcription factors, transporters, stress, cell wall and hormone signalling were differentially expressed during coleoptile senescence, respectively (  Table 4). Detailed analysis of the selected categories revealed that out of the 176 signalingrelated genes only 3 genes were specifically expressed in senescing coleoptiles. A greater number of signaling genes displayed downregulation (108) as compared to those that exhibited upregulation (65)  Coleoptile senescence in rice. (A) Different stages of coleoptile senescence in rice seedlings after transfer from anaerobic (submergence) to aerobic conditions. (B) Expression profile of senescence associated genes (SAGs) at different stages of coleoptile senescence using quantitative PCR. For normalization eEF-1a was used as an internal control. Three biological replicates and two technical replicates were included in the study. White arrows highlight the progression of senescence in the coleoptile. Asterisks indicate a significant difference between the control and senescence tissues (Student's t-test, *P < 0.05, **P < 0.01). RCC reductase: red chlorophyll catabolite reductase.
Out of the 143 stress-associated genes, 56 belonged to abiotic stresses category (37 upregulated, 15 downregulated and 4 senescence-specific), while 87 were related to biotic stresses (44 upregulated, 34 downregulated and 9 senescence-specific) ( Figure 3E; Supplementary Table 4). One of the limitations of the MapMan tool is that many annotated genes are not associated with any specific pathway and are listed as 'unspecified'. Consequently, 56 and 19 genes under biotic and abiotic stress categories, respectively, were not assigned any sub-category. However, from the ones that were assigned a sub-category, PR proteins (13 upregulated, 12 downregulated and 3 senescencespecific) constituted a large number of differentially expressed biotic stress genes. In case of abiotic stresses, many heat stress (13 upregulated, 4 downregulated) and drought/salt stress (5 Functional annotation of DEGs using MapMan tool. The functional annotation of DEGs were analysed using MapMan tool and overviews of metabolism (A) and regulation (B) are depicted. In the overviews, red boxes represent up-regulated/senescence specific transcripts and the blue boxes represent down-regulated transcripts. The major functional categories of genes: Transporters (C), Transcription factors (D), Stress (E) and Cell wall-related (F) are represented as stacked bar graphs. In the stacked bar graphs red bar indicates upregulated genes, green bar shows down regulated genes and orange bar represents senescence-specific genes (day 2-specific). upregulated, 6 downregulated) related transcripts showed altered expression during coleoptile senescence. Because of extensive changes that occur in the cell wall during senescence (Serafini-Fracassini et al., 2002), we also analyzed genes associated with the term 'cell wall' and found that 94 genes associated with the cell wall exhibited differential expression. These included 42 downregulated, 40 upregulated, and 12 senescence-specific genes. Overall, genes associated with the degradation and cellulose synthesis were largely up-regulated, whereas genes associated with modification and cell wall proteins were down-regulated during senescence ( Figure 3F; Supplementary Table 4).

Comparison of transcriptome datasets of coleoptile and leaf senescence in rice
To delineate the core and specific components involved in the two physiologically different senescence programs in plants, we compared the DEGs of coleoptile transcriptome (3439 DEGs) with the publicly available combined transcriptome data of second and flag leaf senescence, henceforth referred to as LS (3778 DEGs; Lee et al., 2017). While the majority of genes were specifically up-and down-regulated in LS and CS (1603, 1691 up-regulated and 1638, 913 down-regulated in LS, CS, respectively), a smaller, albeit significant, number of genes exhibited similar expression pattern in both LS and CS (194 showed upregulation and 85 were downregulated). These genes probably constitute the core components of senescence in rice. Based on these observations it is evident that CS is remarkably distinct from the LS program in rice. This is further supported by an inverse expression pattern of a substantial number of genes (258) in LS and CS ( Figure 4A; Supplementary Table 5). To obtain insights into the biological function of core components of the senescence process, GO analysis of the commonly upregulated genes was carried out. The analysis revealed that maximum number of upregulated genes were associated with cellular components (enrichment of subcategories such as "integral components of membrane", "nucleus" and "membrane") followed by molecular function (subcategories such as "DNA binding" and "TF activity") ( Figure 4B; Supplementary Table 5). Functional enrichment analysis identified several transporters (amino acid transporters, aquaporins, MATE efflux protein, ion transporters), TFs (MYB, NAC, homeobox domain, Zinc finger proteins, and WRKY family proteins), cytochrome P450, heat shock transcription factors or HSFs (HSFA7a, HSFA3a and HSFC1b), hormone signaling pathway genes (ABA, auxin, ethylene and gibberellin) and chlorophyll degradation genes (RCC reductase and senescence-inducible chloroplastic stay green 1 protein) among the commonly upregulated genes ( Suppleme ntary Table 5). Genes encoding for 1aminocyclopropane-1-carboxylate oxidase proteins, TFs such as MYB, bZIP, WRKY118, FBXs (FBX27, FBX 28), receptor like protein kinases, lectin-like receptor kinases and certain enzymes such as GDSL-like lipase/acid hydrolase and glycosyl hydrolase showed down regulation in both CS and LS. Among the genes that exhibited inverse correlation in expression during LS and CS, a large number of genes encoded for proteins belonging to subcategory "integral component of membrane" followed by 'plastid'. Besides, during CS several genes associated with "secondary metabolite synthesis", "photosynthesis", "metal ion binding", "oxidoreductases" and "response to stress" categories were significantly upregulated whereas they exhibited downward trend in expression in LS. On the other hand, several genes coding for proteins of subcategories "plasma membrane", "DNA binding" and "transcription" showed significant induction in LS while their expression declined in CS (Supplementary Table 5). It is noteworthy that several membrane proteins, especially transporters (amino acid transporter, aquaporin, malic acid transporter, phosphate transporter and potassium transporter), were commonly upregulated, while some membrane proteins (citrate transporter, amino acid transporter) exhibited inverse correlation in their expression in CS and LS.
Validation of rice coleoptile senescence transcriptome data RNA-seq derived expression data was reconfirmed with the help of qRT-PCR by examining the expression of 8 DEGs. Relative expression of genes belonging to functional categories such as transporters (amino acid, ABC transporter), signalling (receptor kinase), hormone signalling and transduction (hydrolase, cytochrome P450 and 12-oxophytodienoate), cell wall (pectinesterase), and photosynthesis (ribulose bisphosphate carboxylase) was determined in control and at day 1 and day 2 of senescing coleoptiles. Both transporters and hydrolase genes showed consistent upregulation in day 1 and day 2 of senescence ( Figures 5A, B, D). Expression of receptor kinase, cytochrome P450 and 12-oxophytodienoate increased massively on day 1, followed by a slight decline (although upregulated) on day 2 (Figures 5C, E, F). Both pectinesterase and ribulose bisphosphate carboxylase showed significant downregulation during progression of senescence. However, pectinesterase exhibited a more severe decline in expression levels on day 2 as compared to ribulose bisphosphate carboxylase ( Figures 5G, H). The expression patterns of the genes by qPCR analysis were consistent with that of the RNA-sequencing results.
Identification, digital expression profiling and validation of conserved and novel miRNAs during coleoptile senescence Small RNAs are known regulators of gene expression in diverse biological processes in plants such as growth, development and response to biotic and abiotic stresses . The role of miRNAs in regulating induced and natural senescence in leaves has been demonstrated by several studies (Huo et al., 2015;Kim et al., 2018), however, no miRNAs have been identified which have a potential role in controlling the progression of senescence in coleoptile, which is the first organ to senesce in cereals. Therefore, small RNA libraries of Comparison of DEGs identified in coleoptile senescence with those reported in natural leaf senescence in rice. Leaf senescence data (both second and flag leaf) was retrieved (Lee et al., 2017) and was compared with coleoptile senescence dataset obtained in the present study. DEGs identified in the 36 days post heading (36 DPH) tissues of LS dataset were used for the analysis. (A) Venn diagram of DEGs shows overlap between CS and LS. (B) GO enrichment analysis of DEGs which were commonly upregulated during leaf as well as coleoptile senescence in rice.
Similar to the conserved miRNAs, novel miRNAs were also regulated during coleoptile senescence. A total of 21 novel miRNAs were differentially expressed, of which 9 miRNAs were upregulated and 12 miRNAs exhibited downregulation ( Figure 6B; Supplementary Table 7). Among the upregulated novel miRNAs, miRNA_16, miRNA_24, miRNA_26, miRNA_79, and miRNA_82 showed more than 4-fold upregulation. Further, greater than 3-fold downregulation was observed in the expression levels of miRNA_18, miRNA_22, miRNA_25, miRNA_45, miRNA_47, miRNA_49, miRNA_52 and miRNA_75 ( Figure 6B; Supplementary Table 7). Expression Quantitative PCR-based validation of DEGs associated with coleoptile senescence in rice. Relative expression of genes belonging to functional categories such as transporters: ABC transporter (A), amino acid transporter (B); signaling: Receptor kinase (C); hormone signaling and transduction: Hydrolase (D), cytochrome P450 (E) and 12-oxophytodienoate (F); cell wall: Pectinesterase (G); and photosynthesis: Ribulose bisphosphate carboxylase (H) was determined in control and at day 1 and day 2 of senescing coleoptiles. For normalization eEF-1a was used as an internal control. Three biological replicates and two technical replicates were included in the study. Asterisks indicate a significant difference between the control and the senescence sample (Student's t-test, **P < 0.01).
of 5 conserved and 5 novel miRNAs was validated by poly(T) adaptor RT-PCR method. It was observed that expression of osa_miR1428e-3p and osa_miR171h elevated on day 1 and although their levels declined slightly on day 2, they remained upregulated ( Figure 6C). On the other hand, levels of osa_miR396c-5p were lower on day 1, but increased rapidly on day 2 of senescence ( Figure 6C). Among the downregulated conserved miRNAs, osa_miR398b exhibited significantly higher decline than osa_miR319a-3p.2-3p during senescence ( Figure 6C). The expression of novel miRNAs, miRNA_24 and miRNA_27 increased moderately, while miRNA_47 and miRNA_49 downregulated significantly during coleoptile senescence ( Figure 6D). Notably, the expression pattern of conserved and novel miRNAs as observed by qRT-PCR was similar with their digital expression profiles.

Comparison of coleoptile and rice flag leaf senescence miRNAs
To gain insights into the conservation or differences in the miRNA-dependent regulation of senescence in coleoptile and Expression profiling of miRNAs that are differentially expressed during coleoptile senescence in rice. Heat map of differentially expressed miRNAs belonging to 'conserved miRNAs' (A) and 'novel miRNAs' (B) categories during coleoptile senescence in rice. Quantitative PCR-based validation of few conserved (C) and novel (D) miRNAs. For normalization, 5S ribosomal RNA was used as internal control. Three biological replicates and two technical replicates were included in the study. Asterisks indicate significant differences between the control and the senescence sample (Student's t-test, *P < 0.05, **P < 0.01).
flag leaf, expression of 66 conserved miRNAs was compared with the 38 miRNAs previously identified by our group in flag leaf senescence (FLS) (Sasi et al., 2019). The analysis revealed that 17 conserved miRNAs were detected in both CS and FLS (Supplementary Table 8). While 3 conserved miRNAs (osa-miR164a, osa-miR171b, osa-miR535-3p) showed similar expression patterns, 6 miRNAs (osa-miR1428e-3p, osa-miR160e-5p, osa-miR168a-5p, osa-miR171h, osa-miR1432-5p, and osa-miR1861e) displayed an inverse expression pattern in CS and FLS datasets. The remaining 8 miRNAs did not display any remarkable comparative pattern during the two senescence programs. We further compared the expression trends of the 41 true novel miRNAs in CS with 203 true novel miRNAs in FLS and found that 3 miRNAs [(miRNA_2, 24 and 30 which were named as N_osa_619, N_osa_595 and N_osa_176, respectively, in Sasi et al. (2019)] were present in both the datasets. Out of the 3 commonly present novel miRNAs, miR_30 was slightly downregulated in both CS and FLS, while miR_23 was slightly downregulated in FLS, but was upregulated in CS. miR_2 exhibited no significant change in expression level in FLS, however, it was downregulated in CS. The comparative study led to the identification of few miRNAs which probably comprise the core regulatory pathway of senescence programs in rice. It would be worthwhile to determine the expression of their putative targets in the senescence datasets to identify the functionally conserved miRNA-mRNA modules regulating senescence in different tissues of rice.

Computational prediction of targets of miRNAs, their GO analysis and expression profiling-based validation of selected targets
To elucidate the biological function of miRNAs and other small RNAs it is pertinent to identify the genes targeted by them. The information generated is utilized to unravel the complex regulatory network of miRNA-target interaction and their possible role(s) in regulating different biological processes. We predicted 301 targets for 62 miRNAs (41 conserved and 21 true novel) that were differentially expressed during coleoptile senescence in our dataset (Supplementary Table 9). Targets could not be predicted for four conserved (osa-miR169b, osa-miR390-3p, osa-miR394 and osa-miR5788) and one novel (miRNA_47) miRNAs. Except osa-miR1428e-3p, osa-miR169f.1, osa-miR169h, osa-miR319a-3p.2-3p, osa-miR396a-3p, osa-miR5510 and miRNA_23 which were predicted to hit only a single target gene, all other miRNAs were predicted to target multiple genes. The maximum number of targets were 44 for osa-miR531a, followed by 14 for osa-miR156j-3p. Gene ontology clustering analysis of the predicted targets revealed that the maximum number of targets were associated with molecular function (207) followed by biological process (145) (Supplementary Table 9). The top sub-categories in each gene function cluster belonged to ATP-binding, transferase activity, membrane proteins, protein phosphorylation, oxidationreduction process and regulation of transcription. All these observations indicate that these miRNAs might regulate diverse processes in coleoptile senescence program.
To validate the predicted regulatory relationship between the senescence-associated miRNAs and their target, we performed qRT-PCR of osa-miR164d, osa-miR169r-3p, miRNA_45 and miRNA_75 and their corresponding putative target genes. The expression patterns of all miRNAs and their targets exhibited consistency with the digital expression data obtained in RNAseq and small RNA-seq experiments (Figure 7). osa-miR164d was predicted to target OsNAC2 or OsORE1 (encoded by LOC_Os04g38720) which was significantly upregulated during progression of coleoptile senescence. Similarly, osa-miR169r-3p exhibited downregulation and its target ABC transporter (encoded by LOC_Os11g07600) accumulated to high levels. One of the targets of novel miRNA_45 was MYB TF encoded by LOC_Os05g07010 whose expression levels were upregulated. The expression profiling analyses of miRNA-mRNA interactions indicated the high confidence level of the prediction of the targets and the requirement of testing more such pairs to gain an insight into the regulation of coleoptile senescence in rice.

Regulatory network of miRNAs and their target genes
To identify the key regulatory components and their interactors in coleoptile senescence, a combinatorial approach was employed to integrate the regulatory dataset of RNA-seq and small RNA-seq with publicly available protein-protein interaction (PPI) dataset (Szklarczyk et al., 2021). Expression of putative miRNA-targets was searched in the transcriptome data and the target genes whose expression correlated inversely with that of corresponding miRNAs were interrogated for their protein interactions (Supplementary Tables 9, 10). The PPI analysis output was used to generate miRNA-target regulatory networks using the Cytoscape software (Figure 8). It was noted that 105 known miRNA-target pairs and 35 novel miRNA-target pairs showed an inverse correlation in the expression profile ( Figure 8). The biological regulatory network model consisted of 14 and 28 modules composed of upregulated and downregulated miRNAs targeting the genes that exhibited significant decline and elevation in expression level during coleoptile senescence, respectively (Figure 8). An important module was controlled by an upregulated miR531a whose all the 15 target genes exhibited decline in their expression levels during senescence. Among these 15 targets, two genes encoded for LRR transmembrane protein kinase and an expressed protein and they are known to interact which has been validated experimentally (Szklarczyk et al., 2021). It would be worthwhile to validate the targets and identify which of these contribute to coleoptile senescence. Similarly, other downregulated conserved miRNAs such as miR397a and miR156j could target 9 genes each, while miR164a was found to target 7 genes, expression of all of which was suppressed during senescence. Two novel miRNAs, miRNA_45 and miRNA_75 were predicted to target 6 genes each and these genes also exhibited an inverse correlation in expression. Majority of the modules in the regulatory network consisted of senescence-repressed miRNAs and their corresponding senescence-upregulated targets and therefore protein-protein interactions were largely predicted within these modules.
At the center of the predominant hub were two members of the MIR164 family, miR164a and miR164d. The sequences of miR164a and miR164d are similar except that miR164a had an extra nucleotide at its 3' end. However, this did not affect the ability of these miRNAs to target the same genes ( Supplementary  Tables 9, 10). Interestingly, proteins coded by two of their targets i.e. LOC_Os01g62490 (coding for a laccase precursor) and LOC_Os04g38720 (a NAC TF), appeared to be interacting with each other. These proteins further served as connectors with four additional modules; a) NAC protein interacted with OsSPL2 protein encoded by LOC_01g69830 (a target of miRNA_41 and osa-miR529b; b) laccase precursor interacted with a plastocyanin-domain containing protein encoded by LOC_06g11490 (a target of osa-miR528-3p); c) laccase precursor is also a target of osa-miR397a, thereby connecting the osa-miR164 and osa-miR397a nodes; d) laccase precursor interacted with an auxin-responsive protein encoded by LOC_Os05g48270 (a target of osa-miR169n). Additional connections between different modules were seen such as the Molybdenum cofactor biosynthesis protein 1 encoded by LOC_Os02g15870 (target of osa-miR529b) interacted with CTP synthase protein encoded by LOC_Os01g46570 (target of osa-miR156j-3p). Intermodular connection was also observed between the module regulated by miRNA_75 and module regulated by miRNA_41 and osa-miR529b, wherein 1,3-betaglucan synthase component domain containing protein (encoded by LOC_Os03g03610) was a target of both osa-miR529b and miRNA_75. Intramodular connections were seen in the module regulated by miRNA_41 and osa-miR529, both of which shared two targets i.e., LOC_Os01g69830 and LOC_01g22954. Only few interactions among different proteins encoded by miRNA-targeted genes in the dataset were observed. However, most of the PPIs identified in the network are predicted on the basis of their co-expression or database or text mining. It is, therefore, necessary to validate these interactions experimentally and then characterize their role in coleoptile senescence in rice. It is expected that with the enrichment of PPI data in the public domain, more connectivity between the modules (both inter-modular and intra-modular) will be observed which would highlight the functional and dynamic complexity of proteins and miRNA regulation of their target genes.

FIGURE 7
Validation of expression profile of miRNAs and their predicted target genes during coleoptile senescence in rice. Expression analysis of selected miRNAs (A) and one of their predicted targets genes (B) using qPCR. For normalization of real-time PCR results, eEF-1a was used as internal control. Three biological replicates and two technical replicates were included in the study. Asterisks indicate significant differences between the control and the senescence sample (Student's t-test, *P < 0.05, **P < 0.01). Sasi et al. 10.3389/fpls.2022.985402 Frontiers in Plant Science frontiersin.org

Discussion
Senescence is a precise and highly orchestrated degradation process that is essential for reallocation of resources to developing organs and overall fitness of the plants. The genetic complexity and molecular mechanisms underlying natural leaf senescence have been well-studied, however, owing to slow progression of flag leaf senescence and long experimental duration with flag leaf being the last leaf to emerge, studying flag leaf senescence is challenging. On the other hand, coleoptile the first organ to senesce after seed germination in monocotyledonous plants undergoes a rapid transition from growth to death phase and due to which there is more uniformity in the coleoptile senescence as compared to that of leaves (Kawai and Uchimiya, 2000). Nevertheless, the key changes at tissue and cellular levels associated with senescence in coleoptiles and leaves are conserved (Inada et al., 2002). The unique features of coleoptile senescence thus advocates the possibility of its utilization as a convenient and rapid system for studying the cellular and molecular mechanisms of senescence in cereals. With an aim to unravel the functional biological networks of coleoptile senescence in rice, an integrative multi-omics approach was designed which combined high-throughput RNA-seq and small RNA-seq data of the senescing coleoptile generated in the present study with publicly available protein-protein interaction (PPI) datasets.
Our study reports identification of 3166 differentially expressed and 273 senescence-specific genes (jointly referred to as DEGs) in rice coleoptile senescence. Several of these DEGs encodes for transporter proteins, more interestingly, proportion of upregulated transporters (112) was nearly 5 times higher than those exhibiting downregulation (22). Transport of ions and other macromolecules across cellular boundaries is a key event that ensures remobilization of nutrients from source to sink (Guo et al., 2021). Previous studies on elucidating the molecular changes occurring during leaf senescence reported an increased Integration of regulatory miRNA-mRNA network with PPI network during coleoptile senescence in rice. Protein-protein interaction network for the predicted target genes of conserved and novel miRNAs was generated using Cytoscape visualization software. The expression of miRNAs and target genes is presented based on the RNA-seq and small RNAome data obtained for coleoptile senescence in rice, respectively. Upregulated and downregulated miRNAs are represented by red and green boxes, respectively. Upregulated and downregulated predicted target genes are represented by yellow and white boxes, respectively. The predominant hub identified in the present study is highlighted. expression of genes encoding for several transporters, especially ABC and MATE type, amino acid and ion transporters (Buchanan-Wollaston et al., 2005;Graaff et al., 2006;Lee et al., 2017). ABC, peptides, amino acids, metabolites, sugar, and metal transporters were abundantly expressed during CS, which is in agreement with earlier reports (Guo et al., 2004;Zhang et al., 2014). Among the upregulated transporters, 11 genes each for MATE (Multidrug And Toxic compound Extrusion) and aquaporins (AQPs) were identified. MATE protein family in rice has 46 members and are involved in transport of secondary metabolites, ions and phytohormones and are therefore associated with regulation of several agronomic traits by m e d i a t i n g g r o w t h , d e v e l o p m e n t a n d r e s p o n s e t o environmental stimuli (Upadhyay et al., 2019;Du et al., 2021). Overexpression of ELS1 (a MATE transporter) in Arabidopsis accelerated the senescence (Wang et al., 2016). Further, ABS3 (abnormal shoot 3), a MATE transporter protein, interacts with ATG8 (autophagy-related protein 8) and promotes senescence under natural and carbon-deprivation conditions (Jia et al., 2019). Functional characterization of the differentially expressed MATE transporters will help in identifying the transported molecules during CS. Another transporter family which attracted our attention was aquaporins (AQPs), which are expressed under a myriad of physiological conditions and regulate plant development and stress in plants (Maurel et al., 2008). These channel proteins are actively involved in water transport and trafficking of small solutes such as silicon, urea, reactive oxygen species (ROS) and even gasses such as CO 2 and ammonia (Maurel et al., 2008). Interestingly, CS involved upregulation and downregulation of 11 and 1 aquaporins, respectively. Functionally, AQPs are involved in plant senescence as the suppression of TIP1;1, a member of tonoplast AQPs, compromised carbohydrate transport, had excessive water loss and accelerated leaf senescence (Ma et al., 2004), while overexpression of the cotton TIP1;1-like protein in Arabidopsis promoted premature bolting and delayed senescence of rosette leaves (Cheng et al., 2022). AQPs promote diffusion of H 2 O 2 across the membranes during senescence which may act as a signal for triggering senescence in the adjoining cells (Zentgraf et al., 2022). During senescence, generation of ROS, especially H 2 O 2 is one of the early responses and it acts as a potent signaling reactive oxygen species. H 2 O 2 further reprograms expression of AQPs and other SAGs. In addition to regulating the expression, H 2 O 2 is also believed to modulate the AQPs transportation efficiency, structure and localization (Luu and Maurel, 2005). Owing to the multitasking roles of different transporters in plants, it would be worth exploring the function of these transporters in plant senescence.
Members of TF families such as NAC, MYB, WRKY, TCP, AP2/EREBP, bZIP and bHLH are associated with senescence (Guo et al., 2004;Gregersen and Holm, 2007;Breeze et al., 2011;Zhang et al., 2014;Lin et al., 2015;Moschen et al., 2016;Hinckley and Brusslan, 2020;Dong et al., 2021). The number of upregulated members of TF families AP2/EREBP, Zinc finger, homeobox, MYB, bHLH and WRKY TFs was greater than the down regulated members during CS, thus conforming to the previous observations (Zhang et al., 2014;Lin et al., 2015). Heat shock transcription factors (HSFs) have been indirectly associated with senescence and an early senescence phenotype was reported in the knock-out mutant of AtHSFB1 in Arabidopsis (Breeze et al., 2008). Wu et al. (2012) demonstrated that JUNGBRUNNEN1 (JUB1), a NAC TF, functions as a negative regulator of senescence through its action on HSFA3 which elevates accumulation of HSPs and lowers H 2 O 2 levels (Wu et al., 2012). Our results indicate that HSF family members play a more crucial role in regulating senescence as 5 members of HSF family (OsHSFA2b, OsHSFA3a, OsHSFA7a, OsHSFB2b and OsHSFC1b) were significantly upregulated during CS. Comprehensive expression profiling of HSFs has also revealed differential regulation of several HSFs in natural and induced flag leaf senescence in rice (Sasi et al., unpublished). Further, we found that the putative promoter regions of approximately 5% of SAGs in rice possess canonical and non-canonical heat shock elements (HSEs), thus indicating that HSFs might regulate the expression of SAGs by binding to their promoters (Sasi et al., unpublished).
Phytohormone homeostasis genes, together with an active involvement of regulatory factors such as TFs and small RNAs, fine tune senescence program in plants by integrating the developmental and environmental signals (Staden, 1995;Schippers et al., 2007;Khan et al., 2014). We observed exclusive upregulation for several SA and JA pathway genes. Additionally, the proportion of upregulated genes associated with cytokinin metabolism were largely downregulated, but those related with the auxin metabolism were upregulated. The levels of cytokinins decline during leaf senescence and exogenous application of cytokinins delay senescence (Guo et al., 2021). Although, upregulation of a few IAA (indole acetic acid) biosynthesis genes during leaf senescence has also been previously observed (Quirino et al., 1999), auxin also appears to play a negative role by suppressing expression of WRK57 and its downstream target SAG12, thereby delaying senescence in Arabidopsis (Noh and Amasino, 1999). A comparative transcriptome study revealed that SA, JA and ethylene are predominantly associated with developmental and dark-induced leaf and cell-suspension senescence, respectively (Buchanan-Wollaston et al., 2005). Transcript profiling of hormone pathway genes in senescing leaves in Arabidopsis and cotton demonstrated the involvement of several hormones such as salicylic acid (SA), ethylene, auxin, jasmonic acid and cytokinin, ABA and BR (Graaff et al., 2006;Lin et al., 2015). BRs accelerate senescence in detached cotyledons and leaves in wheat and pea (Saglam-Çag, 2007;Fedina et al., 2016) and acts synergistically with auxin to induce senescence in soybean cotyledons (Baris and Saglam-Cag, 2016).  (Pruzinskáet al., 2005;Kusaba et al., 2007;Park et al., 2007;Schelbert et al., 2009;Kim et al., 2011;Yang et al., 2011;Liang et al., 2014;Takasaki et al., 2015;Mao et al., 2017).
Unfavorable growth conditions such as salinity, extreme temperatures, dehydration, oxidative and nutrient deprivation induce premature senescence in plants. In fact, gene expression analyses have shown that several stress-responsive genes, including regulatory genes, are induced during progression of senescence (Breeze et al., 2011;Dong et al., 2021;Li et al., 2021). It is, therefore, believed that the molecular pathways of leaf senescence cross-talk with stress responsive pathways in plants (Guo and Gan, 2012). Our study revealed that several PR genes, conventionally associated with biotic stress, were differentially expressed in senescing coleoptiles, which is consistent with the previous studies where PR protein encoding genes were categorized as SAGs (Barth et al., 2004;Sillanpää et al., 2005). Accumulation of PR proteins was also observed in the intercellular spaces of barley leaves during natural senescence (Tamaś et al., 1998). Protein profiling of apoplastic proteins of naturally, but not dark-induced, senescing leaves of Arabidopsis revealed that the most abundant SAG proteins were PR2 and PR5 (Borniego et al., 2020). Pathogens interfere with the host development by modifying the signalling pathways and regulating progression of plant senescence to meet their nutritional demands, whereas plants counteract by inducing a hypersensitive response to kill its own cells. When plants are challenged by a pathogen, induction of PR proteins along with several TFs is necessary for triggering the defense-related responses (Häffner et al., 2015). It would be worth comparing the induction patterns of PR genes in genotypes with contrasting senescence timing as well as correlating the progression of senescence with steady state levels of PR genes in pathogenresistant versus susceptible genotypes. It is, therefore, pertinent to understand how the complex interaction of host and pathogen triggers senescence so as to pinpoint the cross-talk, signaling branches and regulatory nodes that play crucial roles in both pathogen response and senescence. Under abiotic stress category, heat stress-related genes such as DnaJ-related heat shock proteins (HSPs), several HSP70s, small HSPs (HSP18.2, HSP17.4, mitochondrial HSP23) exhibited significant upregulation during senescence. Although the potency of heat shock response (HSR) declines at the terminal stages of development, several components of HSR accumulate specifically during aging (Taylor and Dillin, 2011). One class of such proteins is molecular chaperones, which prevent protein aggregation observed during senescence and stress (Calderwood et al., 2009;Trivedi and Jurivich, 2020). Induction of HSFs and HSPs during stress and senescence thus indicate the existence of cross-talk between the senescence and abiotic stress pathways.
Although senescence of coleoptile and flag leaf supplies nutrition to growing seedlings and developing seeds respectively, the overall physiology of senescence in these tissues is vastly different. Whereas coleoptile emerges from the seed and exhibits rapid senescence, second leaf and flag leaf emerge from photosynthesizing plants and have a slow progression of senescence. To explore the extent of similarity and highlight the core components involved in senescence of two spatially and temporally different organs in rice (coleoptile and leaves) DEGs of CS and LS were compared. We observed that approximately 92% of the DEGs had different expression profiles in CS and LS, which is a reflection of their divergent physiologies at the level of gene expression. Among the 8% genes that were similarly regulated, 3 AQP-coding genes were commonly induced during both leaf and coleoptile senescence, however, overrepresentation of upregulated AQPs in CS (11) as compared to LS (4) also highlight specific roles of the aquaporins during coleoptile senescence. It would be worth performing additional functional studies to understand the convergent and divergent roles that the common and specific AQPs play during senescence in rice. Among the TFs, few members such as WRKYs, MYBs, NACs and Zinc finger proteins along with 3 HSFs (OsHSFA3a, OsHSFA7a and OsHSFC1b) were commonly upregulated in CS and LS. The roles of transporters and TFs during plant senescence have been discussed previously in this article. Since AQPs, HSFs and other TFs are also involved in stress-responsive pathways, stress and aging pathways appear to be closely related.
Several genes encoding ELIPs (early light-induced proteins) were also upregulated in both CS and LS. ELIPs bind to free chlorophyll released from pigment-protein complexes and provide protection against oxidative stress during senescence (Binyamin et al., 2001;Hutin et al., 2003). Similarly RCCR (red chlorophyll catabolite reductase) induced in both types of senescence is involved in chlorophyll degradation, an integral component of plant senescence (Pruzinskáet al., 2005). Another important gene associated with the initial step of the chlorophyll degradation pathway, OsSGR (stay-green), was also expressed at high levels during senescence in coleoptiles and flag leaves in rice. OsSGR encodes a vital enzyme, magnesium (Mg) dechelatase that removes Mg from chlorophyll-a resulting in its conversion to pheophytin A (Phein A), an intermediate in the universal chlorophyll breakdown pathway (Shimoda et al., 2016). Furthermore, SGR binds to LHCII and accelerates the breakdown of chlorophyll-protein complexes (Sakuraba et al., 2012). Several CYPs (cytochrome P450 monooxygenases) exhibited upregulation in both CS and LS. CYPs are known as 'versatile biocatalysts' because of their involvement with biosynthesis of antioxidants, secondary metabolites and phytohormones (Pandian et al., 2020), whose levels change with progression of senescence. The role of cytochrome P450 members in brassinosteroid signaling is well-documented in the studies on rice dwarf mutant, dwarf11 (Tanabe, 2005). Cell wall breakdown and a decline in photosynthesis are often associated with the senescence process (Mohapatra et al., 2010) and we also observed that several genes belonging to the cell wall such as pectinesterase and photosynthesis-related categories such as Ribulose bisphosphate carboxylase small chain precursor, displayed significant decline in expression levels. Pectinesterase is involved in cellular adhesion and stem elongation  and its level decreases in senescing strawberries (Castillejo et al., 2004). Ribulose bisphosphate carboxylase small chain precursor is the core component of photosynthesis (Bannenberg et al., 2009) and systematic degradation of chlorophyll during senescence could possibly be related to the decreased expression of Ribulose bisphosphate carboxylase small chain precursor. One and two genes encoding for 1aminocyclopropane-1-carboxylate oxidase (ACO) family were upregulated and downregulated, respectively, during both CS and LS. ACO members are differentially expressed developmentally and in tissue-specific manner (Houben and Van de Poel, 2019) and regulate various developmental processes and stress responses in plants (Houben and Van de Poel, 2019). Since ACOs exert a precise control on spatial and temporal levels of ethylene it will be interesting to correlate change in expression of individual ACO members with ethylene levels and senescence programs.
In addition to the genes that were co-expressed, multiple genes showed an inverse expression profile during the two senescence programs. GO enrichment revealed that the top category of genes belonged to the 'integral component of membrane' followed by DNA binding/transcription. Additionally, several TFs belonging to the MYB family, CYPs, GSTs (glutathione-S-transferase), LTPLs (lipid transfer proteins), glycosyl hydrolases, transporters and stressresponsive proteins exhibited differential expression patterns during progression of senescence in two tissues. We believe that finding the downstream constituents of the inversely correlated genes will help in finding the basis of the different physiologies seen in these senescence programs. An important gene involved in ABA synthesis, i.e., NCED1 encoding for 9-cisepoxycarotenoid dioxygenase accumulated in CS, while its levels declined in LS. Modulating expression of NCED genes can affect the timing of senescence e.g., overexpression of OsNCED3, and OsNCED5 genes that are strongly induced by dark, accelerate leaf senescence in rice (Huang et al., 2018;Huang et al., 2019). NCED1 is a rate-limiting enzyme in the ABA biosynthesis pathway and since ABA is a positive regulator of senescence in plants, upregulation of NCED1 explains the positive correlation between the NCED1 expression and CS. However, the reason for its downregulation in LS is not known. The downregulation of NCED1 in LS was observed by comparing the expression levels at 36 days post heading (DPH) with 4 DPH stage and it is possible that NCED1 levels elevate at a later stage of senescence in flag leaves. The effect of modulating the expression levels of other members of NCED gene family on natural and darkinduced senescence in rice is not known and worth exploring. Since ABA is synthesized from b-carotene, we analyzed the expression pattern of carotenoid biosynthesis genes in CS as well. At least 16 genes involved in carotenoid production have been identified in rice (Chaudhary et al., 2010), of which two genes, Phytoene synthase (LOC Os09g38320) and Lycopene єcyclase (LOC Os01g39960), were found to be up-regulated in our dataset. During a-carotene biosynthesis, lycopene є-cyclase catalyzes the cyclization of linear lycopene to produce cyclic lycopene (Bai et al., 2009). Phytoene synthase is the major ratelimiting enzyme in the carotenoid biosynthesis pathway (Simpson et al., 2018). It is well documented that ABA can induce Phytoene synthase to activate the production of its precursors (Simpson et al., 2018). This suggests that phytoene synthase and carotenes may participate in an auto-activating ABA biosynthesis pathway during CS in rice.
Small RNAs play an important role in regulating gene expression during growth and development and response to different stresses in plants . A number of miRNAs associated with regulation of leaf senescence have been identified in multiple plant species Huo et al., 2015;Thatcher et al., 2015;Qin et al., 2016;Wu et al., 2016;Swida-Barteczka and Szweykowska-Kulinska, 2019;Sasi et al., 2019). However, the contribution of miRNAs in controlling coleoptile senescence has not been elucidated so far. Therefore, in the present study miRNAs that are differentially expressed during progression of senescence in coleoptiles and their target genes were predicted. Further miRNAs and their corresponding target genes which exhibited an inverse expression pattern during coleoptile senescence were identified. Overall, 41 conserved and 21 novel miRNAs exhibiting differential expression patterns in senescing coleoptiles were identified. To find out the commonality in the senescence programs of different tissues of rice, we compared the miRNAs identified in the CS and flag leaf senescence (Sasi et al., 2019) and found that 3 miRNAs were commonly downregulated in both datasets. Interestingly, when the expression pattern of the predicted target genes of these miRNAs was determined in CS (present study) and FLS (Lee et al., 2017) RNA-seq datasets, we found that the targets of 2 c omm only downregulat ed miRNAs (osa-m iR164a: LOC_Os04g38720 which encodes for NAC2 and osa-miR535-3p: LOC_Os05g42250 which encodes for cyclic nucleotide-gated ion channel 2 or OsCNGC16) exhibited remarkable upregulation which indicates that these miRNA-mRNA pairs could play an important role in plant senescence. OsCNGC16 is a plasma-membrane localized protein which triggers an increase in calcium levels when plants are exposed to extreme temperatures (heat or chilling) and imparts thermotolerance in rice (Cui et al., 2020). Calcium also plays an important role in triggering senescence in plants as exogenous application of calcium delays senescence in detached leaves (Poovaiah and Leopold, 1973). The importance of calcium in senescence is further evident by the action of calcium-dependent protein kinase 1 (CPK1) in phosphorylating ORE, the master regulator of senescence in Arabidopsis (Durian et al., 2020). It would be interesting to validate these miRNAs and their respective targets in different senescence tissues and perform their functional characterization for possible involvement in regulation of senescence in plants.
Integration of expression datasets of RNA-seq and small RNA-seq of coleoptile senescence unraveled the important miRNA-mRNA modules and provided a thorough understanding of the associated multiple regulatory networks. Only high confidence miRNA-target pairs which exhibited inverse expression patterns in miRNAome and transcriptome datasets generated by us, respectively, were utilized for the construction of regulatory network. Further, publicly available protein-protein interaction (PPI) data was also integrated with the miRNA-mRNA network and visualization of regulatory network unravelled 140 miRNA-mRNA modules, while only 6 protein interactions were identified. Based on the network generated we concluded that: a. A predominant hub consists of 9 upregulated miRNAs, including two members of MIR164 family at the center, targeting 42 genes during coleoptile senescence. Previous studies have reported few of these reliable miRNA-target pairs in rice and other plant species: osa-miR164a/d-NAC2, osa-miR397a-laccase, osa-miR528-plastocyanin, osa-miR529b-SPL2 (Fang et al., 2014;Pan et al., 2017;Tsuzuki et al., 2019;Zhu et al., 2020). Some of these miRNAs such as osa-miR164, osa-miR528 and their targets genes NAC2 and COPPER ION BINDING PROTEIN1 have been implicated in regulation of leaf senescence (Yuan et al., 2015;Mao et al., 2017). SPLs are also associated with regulation of flowering time, phase transition and developmental aging (Jung et al., 2016). b. The coregulatory miRNA-target pairs are connected by the characteristic nature of miRNAs targeting more than one gene or one gene targeted by several miRNAs. For example: While osa-miR164a/d target NAC2 and laccase, laccase is also a target of osa-miR397a. SPL2 and serine carboxypeptidase are targets of both osa-miR529b and miRNA_41. 1,3 beta glucan synthase, being the common target of osa-miR529b and miRNA_75, links the two miRNA modules. The autolysis of cellular components during PCD has been linked to serine carboxypeptidases (Schaller, 2004). However, no reports were found that directly implicate 1,3 beat glucan synthase and laccases in plant senescence.
c. Several targets of miRNAs are TFs which regulate diverse biological processes at different stages of plant development and few of these have been implicated in regulation of senescence such as WRKY, NAC, SPL, MYB, TCP, and F-box protein as discussed earlier (Woo et al., 2001;Schommer et al., 2008;Zhang et al., 2011;Jiang et al., 2014;Liang et al., 2014;Jung et al., 2016). d. miR169n is predicted to target two genes, LOC_Os05g48270 (DOMON domain-containing protein-cytochromeb561/ ferric reductase transmembrane protein/auxin responsive protein) and LOC_Os06g37300 (cytochrome P450 or CYP). In rice, CYPs belong to a large multigene family, the members of which are associated with metabolism of fatty acids, phenylpropanoids and steroids and more importantly, with hormone homeostasis, thereby regulating plant development and stress responses (Wei and Chen, 2018). A recent report demonstrated that a mutant of ELL1 which encodes for CYP monooxygenase displayed high accumulation of ROS and cell death in rice (Cui et al., 2021). DOMON gene was identified as a strong candidate in the QTLs for functional stay-green (FSG) trait in rice (Lim and Paek, 2015). However, it is important to investigate the exact function of LOC_Os05g48270 and how it regulates coleoptile senescence in rice. e. Two copper-containing proteins, laccase (multi-copper) and plastocyanin (mono-copper), are part of the hub wherein one member (LOC_Os01g62490) of laccase family is targeted by three miRNAs, osa-miR164a/d and osa-miR397a, while the two other members (LOC_Os01g61160, LOC_Os01g62480) are targeted by osa-miR397a only. The gene (LOC_Os06g11490) encoding endomembrane plastocyanin is targeted by osa-miR528-3p. In plants Cu-miRNAs such as miR397, miR408, miR528 and miR857 target coppercontaining proteins: laccases, plastocyanin, polyphenol oxidase (PPO), ascorbate oxidase (AAO), amino oxidase (AO), Cu/Zn superoxide dismutases (CSDs; Zhu et al., 2020). These Cu-containing proteins are actively involved in diverse processes such as ROS metabolism, Cu homeostasis and stress tolerance . In Arabidopsis, a pair of interacting proteins belonging to blue copper protein family: plantacyanin (PCY) and SAG14 are induced during dark-induced senescence which is crucial for chloroplast copper efflux and this module is regulated by miR408 and PIF3/4/5 (Hao et al., 2022). Plant laccases are Cu-oxidases which catalyze lignin and anthocyanin biosynthesis and are involved in diverse plant processes: growth and development, defense against biotic and abiotic stresses (Ranocha et al., 2002;Wan et al., 2022). However, as per our understanding no studies have directly linked laccase with regulation of senescence in plants. In the present study we show that three laccase-encoding genes are upregulated and 3 miRNAs that are possibly regulating their expression are clearly downregulated. Our results highlight the possible involvement of laccases in regulating coleoptile senescence in rice, however, the underlying mechanism needs to be explored. It is hypothesized that the copper-proteins might contribute to regulation of senescence by modulating Cu and redox homeostasis in plants.

Conclusions
We employed an integrated omics approach by combining transcriptome and small RNAome sequencing of coleoptile senescence in rice. Analysis of the RNA-seq dataset identified 3439 DEGs many of which were TFs, transporters, and involved in hormone metabolism, redox and stress responses. Comparison of CS dataset with publicly available leaf senescence (LS) dataset discovered a set of core senescence associated genes. Small RNA sequencing identified many differentially expressed miRNAs, targets for some of which displayed an inverse expression profile in transcriptome dataset. The regulatory network generated by integrating the high confidence miRNA-RNA target pairs with PPI dataset was comprised of 42 modules. A distinguished hub of miR164a/d targeting copper-containing proteins, laccases and plastocyanin, indicated that copper and ROS homeostasis regulate coleoptile senescence in rice. Overall, our study reports a comprehensive atlas on the expression dynamics of genes, miRNA regulators and their interactions and provides insights into the fundamental pathways regulating coleoptile senescence in rice.

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.