Programming of Plant Leaf Senescence with Temporal and Inter-Organellar Coordination of Transcriptome in Arabidopsis1[OPEN]

RNA-seq analysis of total and small RNAs throughout the lifespan of Arabidopsis leaves revealed that leaf senescence proceeds with tight temporal and distinctive inter-organellar coordination of transcriptomes. Plant leaves, harvesting light energy and fixing CO2, are a major source of foods on the earth. Leaves undergo developmental and physiological shifts during their lifespan, ending with senescence and death. We characterized the key regulatory features of the leaf transcriptome during aging by analyzing total- and small-RNA transcriptomes throughout the lifespan of Arabidopsis (Arabidopsis thaliana) leaves at multidimensions, including age, RNA-type, and organelle. Intriguingly, senescing leaves showed more coordinated temporal changes in transcriptomes than growing leaves, with sophisticated regulatory networks comprising transcription factors and diverse small regulatory RNAs. The chloroplast transcriptome, but not the mitochondrial transcriptome, showed major changes during leaf aging, with a strongly shared expression pattern of nuclear transcripts encoding chloroplast-targeted proteins. Thus, unlike animal aging, leaf senescence proceeds with tight temporal and distinct interorganellar coordination of various transcriptomes that would be critical for the highly regulated degeneration and nutrient recycling contributing to plant fitness and productivity.

Most organisms undergo age-dependent developmental changes during their lifespans. The timely decision of developmental changes during the lifespan is a critical evolutionary characteristic that maximizes fitness in a given ecological setting (Leopold, 1961;Fenner, 1998;Samach and Coupland, 2000). Plants use unique developmental strategies throughout their lifespans as opposed to animals. In plants, most organs are formed postnatally from sets of stem cells in the seed. In addition, plants are sessile and cope with encountering environments physiologically, rather than behaviorally. Thus, they have developed highly plastic and interactive developmental programs to incorporate environmental changes into their developmental decisions (Pigliucci, 1998;Sultan, 2000).
The leaf is an organ that characterizes the fundamental aspects of plants. Leaves harvest light energy, fix CO 2 to produce carbohydrates, and, as primary producers in our ecosystem, serve as a major food source on the earth. Leaves undergo a series of developmental and physiological shifts during their lifespans. A leaf is initially formed as a leaf primordium derived from the stem cells at the shoot apical meristem and develops into a photosynthetic organ through biogenesis processes involving cell division, differentiation, and expansion (Tsukaya, 2013). In the later stages of their lifespans, leaves undergo organ-level senescence and eventually death. Organlevel senescence in plants involves postmitotic senescence and is a term used similarly as "aging" in animals. During the senescence stage, leaf cells undergo dramatic shifts in physiology from biogenesis to the sequential degeneration of macromolecular and cellular structures, including the breakdown of the photosynthetic machinery, the major source of nitrogen and carbon (Noodén, 1988;Nam, 1997). Degenerative activities occur concomitantly with a massive remobilization of the hydrolyzed molecules to the growing plant parts such as young leaves, developing seeds, and fruits. Thus, the entire lifespan of leaves requires a proper balance between biogenesis and degenerative processes for plants' fitness (Hensel et al., 1993;Smart, 1994;Buchanan-Wollaston, 1997;Lim et al., 2007).
Leaf development involving complex, but highly regulated molecular processes would be better understood as a continual transition of functional and regulatory networks throughout the entire lifespan. Recent molecular genetic studies have identified key gene regulatory networks that control cell proliferation, expansion, and differentiation during the early leaf development, including miR319, CINCINNATA-TCPs, miR396, and GROWTH-REGULATING FACTORs/ GRF-INTERACTING FACTORs (Nath et al., 2003;Palatnik et al., 2003;Horiguchi et al., 2005;Rodriguez et al., 2010;Wang et al., 2011a;Debernardi et al., 2012). Combining molecular genetic analysis with systems biology approaches has also revealed complex gene regulatory networks governing leaf senescence, such as a trifurcate feedforward loop including ETHYLENE-INSENSITIVE 2, miR164, and ORESARA1/ANAC092 . EIN3 and additional NAC transcription factors (TFs) were further found to be involved in the leaf senescence network (Balazadeh et al., 2011;Hickman et al., 2013;Li et al., 2013;Kim et al., 2014a;Qiu et al., 2015). On the other hand, comprehensive profiling of molecular programs has been revealed by time-course transcriptome analysis during early and late leaf development (Beemster et al., 2005;Efroni et al., 2008;Breeze et al., 2011;Andriankaja et al., 2012;Thatcher et al., 2015). For example, transcriptional changes in the genes involved in chloroplast differentiation, such as photosynthesis and chloroplast retrograde signaling, were found to be tightly coordinated with the onset of cell expansion (Andriankaja et al., 2012). The mRNA expression profiling obtained from high-resolution time-course microarray data in Arabidopsis (Arabidopsis thaliana) leaves revealed a detailed chronology of transcriptional changes during leaf senescence (Breeze et al., 2011). Although these previous studies are highly valuable for investing key regulatory genes and inferring gene regulatory networks, the analyses were mostly limited to illustrate the functional dynamics representing the chronology of biological processes during the early or late developmental stages.
Developmental shifts during lifespan involve temporal changes in expression of multiple types of RNA, including noncoding RNAs (ncRNAs) as well as coding mRNAs (Graveley et al., 2011;Brown et al., 2014). ncRNAs found in plants include rRNAs, tRNAs, small nuclear RNAs, microRNAs (miRNAs), and long ncRNAs such as antisense transcripts. In addition, plant cells, unlike animal cells, harbor the following three genomes: nuclear, mitochondrial, and chloroplast. Chloroplasts contain over 4,000 proteins used for its photosynthetic and other physiological functions, yet the chloroplast genome encodes only 88 proteins, and the rest are encoded in the nuclear genome and translocated into the chloroplast. It is an intriguing question: How are the various RNA species transcribed from the three different genomes and regulated during the whole lifespan of leaves that undergo a fundamental physiological shift from biogenesis to degeneration? However, holistic analysis of multiple types of RNA, for the sake of gaining an understanding of the coordinated transcriptional programs of multidimensional transcriptomes, has not yet been systematically undertaken.
Previous transcriptome analyses have been concerned with how biogenesis or degeneration is achieved and mostly illustrated the functional dynamics representing the chronology of biological processes during the early or late developmental stages (Beemster et al., 2005;Efroni et al., 2008;Breeze et al., 2011;Andriankaja et al., 2012;Thatcher et al., 2015). In this study, we explored how leaf transcriptome is differentially regulated during two distinctive developmental stages, growth versus senescence, by performing multidimensional analyses of Arabidopsis transcriptomes over the entire lifespan of the leaf organ using strand-specific total RNA and small RNA sequencings (total RNA-seq and smRNA-seq, respectively). Unlike animal aging, the senescence stage in plants showed tighter temporal coordination of the various transcriptomes. Furthermore, the chloroplast transcriptome turned out to be a key constituent in the leaf lifespan program and showed strong interorganellar coordination with the nuclear transcriptome during the senescence process. We also observed preferential production of antisense transcripts from genes for chloroplast-targeted proteins. Our transcriptome data would provide a cornerstone to understand fundamental transcriptional programs underlying age-dependent developmental decisions and functional and regulatory shifts during the lifespan of Arabidopsis leaves.

Multidimensional Analyses of Arabidopsis Leaf Transcriptomes
We generated multidimensional leaf transcriptome data by performing strand-specific total RNA-seq and smRNA-seq analyses of Arabidopsis leaves at various ages (age dimension) to obtain transcriptome information for various protein-coding mRNAs as well as long and small ncRNAs (RNA-type dimension) transcribed from nuclear, chloroplast, and mitochondrial genomes (organelle dimension; Fig. 1, A and B). We sampled the fourth true leaves of Arabidopsis from the leaf ages of 4 to 30 d, encompassing the entire lifespan of a leaf. Sequenced reads from time-course sequencing experiments were aligned to the Arabidopsis nuclear, chloroplast, and mitochondrial genomes (Supplemental Table S1).
In the nuclear genome, 60.7% of all annotated transcripts (24,596 of 40,506 nuclear transcripts in the TAIR 10 database) were expressed during the leaf lifespan (Fig.  1C; Supplemental Table S2), and 46.3% of the expressed nuclear transcripts showed differential expression during the leaf lifespan ( Fig. 1C; Supplemental Table S3). For the chloroplast genome, 87 of the 88 annotated transcripts were expressed in the leaf, and 69 of the 87 were differentially expressed during the lifespan ( Fig. 1D; Supplemental Table S3). For the mitochondrial genome, 103 of 122 annotated transcripts were expressed in the leaf ( Fig. 1E; Supplemental Table S3). A lower portion of the expressed mitochondrial transcripts (26 of 103) than that of the nuclear or chloroplast transcripts was differentially expressed (Fig. 1, C-E). These multidimensional transcriptomes provide a comprehensive and integrated resource for understanding transcriptional programs underlying age-dependent use of diverse types of RNAs from the three organellar genomes.

Functional Transition during the Entire Leaf Lifespan
Exploring the transcriptional programs during the entire leaf lifespan can provide insights into the functional and regulatory dynamics of the two contrasting stages of the lifespan involving early biogenesis and late degeneration processes. For this purpose, we divided the leaf lifespan into the following two stages (Fig. 1A): the growth-to-maturation stage (G-to-M; 4-18 d) and the maturation-to-senescence stage (M-to-S; 16-30 d), when biogenesis and degeneration activities, respectively, are prominent (Oh et al., 1996;Woo et al., 2010Woo et al., , 2013Breeze et al., 2011). We then examined the profiles of transcripts showing age-dependent expression changes at the G-to-M and M-to-S stages. Among the 11,548 differentially expressed transcripts (DETs) during the whole lifespan, 4,827 and 7,723 transcripts were identified as DETs at the G-to-M and M-to-S stages, respectively (Supplemental Table S3). The larger number of DETs at the M-to-S stage than at the G-to-M stage indicates that biochemical and/or physiological transitions are more dynamic during the senescence process. The DETs at each stage were categorized into six major clusters, with three down-regulated (D1-D3) and three up-regulated (U1-U3) groups based on their kinetic expression patterns over time (Fig. 2, A and B), which accounted for approximately 97.7% of the DETs.
Gene ontology biological process (GOBP) enrichment analysis of the major clusters in Figure 2, A and B revealed that processes involved in biogenesis and assembly were predominantly down-regulated during the lifespan (Fig. 2C) as previously reported (Breeze et al., 2011;Brown and Hudson, 2015;Lin et al., 2015). For example, cell proliferation-associated processes (cell cycle, chromatin organization, and translation) were down-regulated at the G-to-M stage. Diverse anabolic processes, including photosynthesis, protein complex assembly, and cell wall organization, were continuously down-regulated from G-to-M to M-to-S, whereas carbohydrate biosynthetic processes were upregulated during G-to-M, but down-regulated later at M-to-S. In contrast, processes related to disassembly were mostly up-regulated throughout the leaf lifespan (Fig.  2C). Catabolic processes of chlorophylls and proteins were up-regulated during G-to-M, whereas lipid catabolic processes were up-regulated at M-to-S. Responses to senescence-associated hormones were up-regulated throughout the lifespan (Fig. 2C). Collectively, our leaf transcriptome data suggest that leaves undergo substantial functional transitions from early biogenesis to late degeneration processes during aging.

Tighter Temporal Coordination among Various Biological Processes at the M-to-S Stage
We then asked whether age-associated processes are temporally coordinated and how the processes are differentially coordinated in the early and late stages of the leaf lifespan. To assess transcriptional coordination among the processes enriched in the two distinctive developmental stages, we estimated the degree of coexpression between all pairs of GOBPs in Figure 2C (Supplemental Table S4). This analysis showed characteristics of transcriptional coordination among key processes at each developmental stage (Fig. 2D). At the G-to-M stage, biogenesis-associated biological processes, such as cell cycle, chromatin organization, and translation, showed higher degrees of coexpression than at M-to-S, suggesting strong transcriptional coordination among these processes during early leaf development (Fig. 2D, green box). Metabolic processes associated with energy production, including photosynthesis, cell redox homeostasis, and generation of precursor metabolites and energy, were tightly coordinated at the M-to-S stage (Fig. 2D, red box). Overall, a much higher degree of coordination in gene expression among biological processes was observed at the M-to-S stage; 8.7% and 36.3% of all GOBP combinations were significantly coexpressed at G-to-M and M-to-S, respectively (Supplemental Table S4). The high degrees of coexpression among multiple biological processes at the M-to-S stage suggest that plants employ regulatory programs enforcing tighter temporal coordination of the transcriptome during aging, which would be necessary for highly ordered senescence processes including macromolecule degradation and remobilization to maximize fitness of plants.

TF-Mediated Transcriptional Regulation Is More Actively Utilized at the M-to-S Stage
To explore how expression of thousands of genes is regulated and coordinated during the leaf lifespan, we investigated the importance of transcriptional regulation by TFs. To this end, we examined the dynamics of TF genes at the two developmental stages during the leaf lifespan. Among 2,403 TF transcripts annotated in Arabidopsis, 484 were up-regulated during the leaf lifespan (197 and 287 at G-to-M and M-to-S, respectively), whereas 293 were down-regulated (106 and 187 at G-to-M and M-to-S, respectively; Fig. 2, A and B). The larger number of differentially expressed TFs at M-to-S than at G-to-M suggests that leaves utilize the changes in TFs more actively at the senescence stage. TF families enriched in major clusters during early and late stages of the lifespan were identified (Fig. 2E). C3H, CCAAT-HAP2, and homeobox TF families were enriched in an up-regulated cluster (U2) at G-to-M, consistent with previous reports of the functions of these TFs in the regulation of early leaf development (Lincoln et al., 1994;Wang et al., 2003;Kim et al., 2014b). The ARF TF family, which regulates auxin signaling (Guilfoyle and Hagen, 2007), was enriched in D3 at M-to-S, supporting the involvement of auxin in senescence-associated processes . Two TF families, WRKY and NAC, were up-regulated at both stages, consistent with their known roles as key regulators in leaf senescence and stress responses (Miao et al., 2004;Guo and Gan, 2006;Kim et al., 2009;Besseau et al., 2012).
TFs differentially expressed (DE-TFs) during the leaf lifespan can contribute to the regulation of ageassociated biological processes. To test this notion, potential target genes of TFs were first obtained from interactions between TFs and target genes previously reported (Yilmaz et al., 2011). Based on these TF-target interactions, transcriptional regulatory circular maps were generated to represent regulatory edges from DE-TFs to their tentative direct targets included in ageassociated GOBPs (Fig. 2C) at the G-to-M and M-to-S stages ( Fig. 2F; Supplemental Table S5). One of the most notable features was that the target genes of DE-TFs were found across most GOBPs, but more significantly enriched in represented GOBPs (red/purple/blue) than nonrepresented GOBPs (gray) in each stage ( Fig. 2F; P-value of 1.5 3 10 26 versus 0.49 at G-to-M and P-value of ,10 27 versus 1.7 3 10 23 at M-to-S). Interestingly, more TF-target interactions were observed at M-to-S than G-to-M (P-values of 0.22 and , 10 27 at G-to-M and M-to-S, respectively), highlighting that TF-mediated transcriptional regulatory programs were more actively utilized for regulating age-associated gene expression in aging leaves. Together, these results suggest that the dynamic expression of TFs contributes to transcriptional regulatory programs for transitions of biological processes during the leaf lifespan.

Genes with Antisense Transcripts Preferentially Encode Chloroplast-Targeted Proteins
Our strand-specific total RNA-seq of Arabidopsis leaves provided information on changes in antisense transcriptomes during lifespan. Of the annotated nuclear genes, 6.9% (2,226/32,440) were transcribed in the antisense direction, whereas chloroplast (94.3%, 83/88) and mitochondrial (65.6%, 80/122) genomes showed higher proportions of genes expressing antisense transcripts (Fig. 3, A and C; Supplemental Table S6). Notably, gene ontology cellular compartment (GOCC) enrichment analysis revealed that genes encoding proteins targeted to chloroplasts were mainly enriched in the 2,226 nuclear genes generating antisense transcripts ( Fig. 3D; Supplemental Fig. S1). A, Representative wild-type Arabidopsis fourth rosette leaves at 2-d intervals from 4 to 30 d after emergence. The leaf lifespan was divided into two stages, G-to-M and M-to-S. B, Multidimensional transcriptome analysis involving three dimensions: age, RNA-type, and organelle. C to E, Proportions of expressed transcripts in the total leaf transcriptome (darkcolored bar) and DETs among the expressed transcripts during the lifespan (light-colored bar) in nuclear (C), chloroplast (D), and mitochondrial (E) transcripts. The numbers of expressed transcripts and the DETs, which were categorized into RNA types, are indicated next to each bar.
Among the total 2,389 genes generating antisense transcripts, we identified 97 (50 up-regulated and 47 down-regulated) and 292 (161 up-regulated and 131 down-regulated) genes with differentially expressed antisense transcripts (DE-ATs) at G-to-M and M-to-S, respectively (Fig. 3E, top panels). We then performed GOBP and GOCC enrichment analyses of genes with DE-ATs. No GOBPs and GOCCs were significantly enriched by the 97 genes with DE-ATs at G-to-M. At M-to-S, photosynthesis and defense response were enriched by genes with down-regulated antisense transcripts, whereas lipid catabolic process and response to senescence-associated hormones, such as ethylene and abscisic acid, were enriched by genes with upregulated antisense transcripts (Fig. 3E, bottom). Interestingly, genes encoding proteins categorized to the photosystem were strongly enriched among genes with down-regulated antisense transcripts at M-to-S ( Fig. 3E, bottom). These findings suggest that the plant leaf utilizes antisense transcripts as a regulatory program for controlling biological processes during the leaf lifespan, particularly for genes involved in chloroplast functions.

Dynamic Changes in Small ncRNA Transcriptomes
During plant growth and development, small ncRNAs play important regulatory roles in diverse cellular processes (Vaucheret, 2006). We identified Arabidopsis smRNAs with sizes of 18 to 48 nt during the lifespan. smRNAs mapped to rRNAs were the most abundant, followed by those mapped to intergenic regions, miRNAs, tRNAs, and mRNAs ( Fig. 4A). Of the smRNA species, 29.7% showed age-dependent abundance changes (Fig. 4B). Notably, miRNA was ranked as the RNA type with the highest proportion among the smRNA species showing age-dependent abundance changes. We examined changes in the composition of smRNA species during the leaf lifespan (Fig. 4C). The abundance of smRNAs mapped to rRNAs decreased significantly during the leaf lifespan, whereas that of smRNAs mapped to tRNAs increased drastically (Fig. 4C, left). The overall abundances of miRNAs and smRNAs mapped to mRNAs decreased and increased, respectively, during the leaf lifespan ( Fig. 4C, right).

The Differential Utilization of Small ncRNA-Based Regulatory Networks during the Leaf Lifespan
We identified Arabidopsis smRNAs with sizes of 18 to 48 nt during the lifespan and examined the dynamics and functions of small ncRNA-based regulatory networks. miRNAs and transacting small interfering RNAs (tasiRNAs) are well-known regulators of various physiological processes in plants (Vaucheret, 2006). It is also believed that 21-nt smRNAs act as regulatory smRNAs via the RNA interference process (Elbashir et al., 2001;Wang et al., 2011b). Thus, we examined agedependent abundance changes in miRNAs, tasiRNAs, and 21-nt smRNAs at G-to-M (4-16 d) and M-to-S (16-28 d). smRNAs with age-dependent changes in their abundance (DE-smRNAs) at M-to-S outnumbered those at G-to-M: 1,196 versus 821 ( Fig. 5A; Supplemental Table S7). Similar numbers of DE-smRNAs showed up-regulation and down-regulation at G-to-M, but a majority of the DE-smRNAs were down-regulated at M-to-S (Fig. 5B). Interestingly, DE-smRNAs at M-to-S showed higher degrees of coexpression than those at G-to-M (Fig. 5C). These results suggest that smRNA regulatory programs are highly coordinated at the M-to-S stage, as shown in ageassociated GOBPs (Fig. 2D).
We further investigated age-dependent regulatory networks involving the 55 miRNAs with age-dependent abundance changes (DE-miRNAs). Using three target prediction tools (Allen et al., 2005;Bonnet et al., 2010;Dai and Zhao, 2011), we identified 265 predicted mRNA targets that showed strong anticorrelation (Pearson correlation # 20.7) with DE-miRNAs in the expression patterns (123 and 200 miRNA-target pairs at G-to-M and M-to-S, respectively; Supplemental Table S8). We then constructed miRNA regulatory networks using the miRNA-target pairs at G-to-M and M-to-S (Fig. 5D). The potential targets at the G-to-M stage were associated with a broad spectrum of biological processes, including responses to abiotic stimuli or plant hormones and chlorophyll or lignin biosynthesis. By contrast, at the M-to-S stage, miRNAs appeared to be involved in biological processes such as response to abscisic acid and disaccharide metabolism. This result indicates that miRNAs contribute distinctly to the two stages by controlling genes involved in different biological processes.
Next, we investigated another well-known class of regulatory smRNAs, tasiRNAs. Using the same method used for miRNAs, 18 and 117 tasiRNA-target pairs were found at G-to-M and M-to-S, respectively (Supplemental Table S8). GOBP analysis of the tasiRNA targets at M-to-S suggested that the decrease in tasiRNAs released the suppression of protein phosphorylation and hexose metabolic processes (Supplemental Fig.  S3). As another potential regulatory RNA, we examined 21-nt smRNAs. We identified 153 and 235 pairs of 21-nt smRNAs and targets at G-to-M and M-to-S, respectively, by adapting the same method as described above and analyzing Arabidopsis degradome sequences (Addo-Quaye et al., 2008;German et al., 2008;Supplemental Fig. S4; Supplemental Table S8). The GOBPs potentially regulated by 21-nt smRNAs were also distinct in the two stages: at G-to-M, predominant GOBPs included responses to abiotic stimulus, salt stress, and auxin stimulus, but at M-to-S, predominant GOBPs included lipid biosynthesis and oxidation/ reduction. Interestingly, the comparison showed that the GOBPs for predicted targets of 21-nt smRNAs were quite different from those for predicted targets of miRNAs. This result suggests that plants differentially utilize miRNAs and 21-nt smRNAs for regulating  Figure 2C. E, Enrichment of TF families in 12 clusters at G-to-M and M-to-S. Scale bar, P-values of enrichment for TF families calculated from the hypergeometric test for DETs in each cluster. F, distinct processes to coordinate transcriptional programs during the leaf lifespan.
TFs were significantly enriched in predicted miRNA targets at the two stages (regulation of transcription in Fig. 5D; P-values of 7.9 x 10 23 and 7.5 x 10 24 at G-to-M and M-to-S, respectively). The miRNA targets included 42 TFs (20 and 29 TFs at G-to-M and M-to-S, respectively; Fig. 5E; Supplemental Table S9). The targets of tasiRNAs and 21-nt smRNAs also included 11 and 19 TFs, respectively (Fig. 5E). Notably, the predicted targets at M-to-S included 10 TFs belonging to the NAC TF family, which is known to play a major role in regulating leaf senescence. This finding indicates that smRNA regulatory networks involving TFs are important in coordinating various developmental processes, including senescence, during the leaf lifespan.

The Chloroplast Transcriptome Is a Key Constituent in Leaf Lifespan Programs
Our data demonstrated that a large proportion of the leaf transcriptome was derived from the chloroplast genome: 36.1% of the whole leaf transcriptome originated from chloroplast-encoded genes, whereas 63.6% and 0.4% were nuclear and mitochondrial transcripts, respectively (Fig. 6A). This finding is quite intriguing, considering the numbers of genes encoded in the nuclear and chloroplast genomes. Moreover, the chloroplast transcriptome exhibited the most pronounced changes during the lifespan: 78.4% of the annotated chloroplast transcripts were differentially expressed, whereas 28.3% and 21.3% of annotated nuclear and mitochondrial transcripts, respectively, were differentially expressed (Fig. 6B). Furthermore, 46.5% of the annotated nuclear transcripts for chloroplast-targeted proteins were differentially expressed during the leaf lifespan, whereas only 22.0% of the annotated transcripts for nonchloroplast-targeted proteins were differentially expressed (Fig. 6C).
One of the unique characteristics of the Arabidopsis transcriptome is the declined amount of chloroplast transcripts during the leaf lifespan. The levels of chloroplast-encoded DETs decreased substantially during senescence (Fig. 6D). To determine how the dynamic changes in the chloroplast transcriptome are coordinately regulated throughout the lifespan, we examined diverse regulators involved in controlling the expression of chloroplast-encoded genes at multiple levels and protein localization in chloroplast. Strikingly, a majority of these regulators showed patterns of decreased expression during the leaf lifespan (Fig. 6E).
These results indicate that overall, expression of regulatory units for chloroplast-encoded genes declined during the leaf lifespan, although not all of the genes encoding chloroplast-targeted proteins showed decreased expression (Supplemental Fig. S5). Our results further indicate that coordinated decline in the chloroplast regulatory units might be responsible for the changes of chloroplast transcriptome along aging.

Tight Coordination between Nuclear and Chloroplast Transcriptomes in Senescing Leaves
Tight communication between the nucleus and organelles is essential for the maintenance of appropriate functions of individual organelles (Susek et al., 1993;Pesaresi et al., 2007). We accordingly examined the degree of coexpression of genes encoding ribosomal proteins (RPs). The majority of RPs are encoded by the nuclear genome but are utilized in the cytoplasm, chloroplast, or mitochondria (Sormani et al., 2011). By contrast, chloroplast and mitochondrial genes encoding RPs are transcribed and utilized only in the corresponding organelles. We compared the correlation coefficients for transcripts encoding RPs ultimately localized in the cytoplasm, chloroplasts, or mitochondria (Fig. 6F). Notably, nuclear and chloroplast genes encoding RPs localized to chloroplasts displayed the highest degree of coexpression compared to those encoding mitochondrial (P-value of , 0.001) or cytoplasmic (P-value of , 0.001) RPs. This result suggests that strong transcriptional coordination between the nucleus and chloroplasts is required for proper chloroplast function.
Photosynthetic complexes include PSI and PSII, cytochrome b 6 f, ATP synthase, and Rubisco complexes. These complexes contain 35 nucleus-encoded and 24 chloroplast-encoded proteins. Strikingly, most of the genes (79 of 88) encoding the components in these five complexes showed shared expression patterns during the leaf lifespan ( Fig. 6G; Supplemental Fig. S6). Comparison of the degree of coexpression between all genes in Figure 6G at G-to-M and M-to-S further revealed higher transcriptional coordination among genes encoding the photosynthetic complexes toward senescence and death (P-value of , 0.001 among N-encoded genes, among C-encoded genes, and between N-and C-encoded genes) when compared with correlation coefficients at G-to-M and M-to-S (Fig. 6H). Our results indicate enhanced transcriptional coordination of genes encoding the components of the photosynthetic complexes during the senescence stage is necessary to ensure optimal transition from energy generation to nutrient remobilization in chloroplasts in the senescence processes. Collectively, these findings suggest that the plant leaf has evolved a system specialized for transcriptional coordination between the nucleus and chloroplasts.

DISCUSSION
The plant leaf undergoes a series of developmental transitions throughout its lifespan in which tight coordination of molecular programs for biogenesis and degeneration is required. In this study, we explored how leaf transcriptome is differentially regulated during two distinctive developmental stages by performing transcriptome analysis of total and small RNAs throughout the entire lifespan of Arabidopsis leaves. Interrogation of our Arabidopsis leaf transcriptome carta provided a comprehensive resource for understanding transcriptional programs underlying the age-dependent utilization of diverse types of RNAs from different organellar genomes. . Antisense transcriptomes during the leaf lifespan. A, Nuclear chromosomal maps of Arabidopsis antisense transcriptomes throughout the leaf lifespan. Genome-wide views of annotated reference genes (TAIR 10, gray) and genes expressed in the antisense direction (red) for all five chromosomes (Chr1-Chr5) are shown. 32,440 is the number of annotated genes, which exclude rRNAs, tRNAs, snoRNAs, and primary miRNAs from the total annotated genes (33,325) in the nuclear genome. In total, 2,226 nuclear genes were expressed in the antisense direction during the leaf lifespan. B, Chloroplast chromosomal maps of Arabidopsis leaf transcriptomes throughout the lifespan. The circular map shows the names of annotated genes (outer line), annotated mRNAs (orange boxes), and genes expressed in the antisense direction at 4, 12, 20, and 28 d from outer to inner tracks. 88 is the number of annotated mRNAs from the total annotated genes (133) in the chloroplast genome. 83 is the number of chloroplast genes expressed in the antisense direction. C, Proportions of genes whose transcripts were found in sense (blue), antisense (red), or both (green) directions to annotated genes in the nuclear, chloroplast, and mitochondrial genomes. D, GOCC enrichment analysis of nuclear genes with antisense transcripts. Scale bar, enrichment P-values. E, K-means clustering of antisense transcript profiles at G-to-M and M-to-S (top) and enriched GOBPs and GOCCs of the genes generating DE-ATs at M-to-S only (bottom). Scale bar, enrichment P-values of the genes with DE-ATs obtained from the hypergeometric tests.
Our multidimensional transcriptome data highlight unexpected and important findings. One of the most notable features that we identified was that senescing leaves showed more coordinated temporal changes in transcriptomes at multiple levels than did growing leaves. A higher degree of coexpression among enriched biological processes (Fig. 2D) as well as among regulatory smRNAs (Fig.  5C) was evident at the late developmental stage.
This finding is intriguing in considering the result that in animals, the integrity of coexpression networks decreases with age (Southworth et al., 2009;Soltow et al., 2010;Viñuela et al., 2010). The difference in the coordination properties of gene expression along aging between plants and animals may reflect that the aging process in plants has evolved to maximize nutrient recycling from old to young or newly developing tissues, whereas aging in animals . smRNA transcriptomes of Arabidopsis leaf during its lifespan. A, Proportion of smRNAs generated from each RNA-type. N, nucleus; C, chloroplasts; intergenic, reads mapped to intergenic regions. B, Proportions of smRNAs with age-dependent abundance changes among smRNA species. C, Changes in the distribution of RNA types of smRNAs during the leaf lifespan. Stacked bar plots illustrate normalized mapped read counts from smRNA-seq data, which were categorized by RNA types, at seven time points. Full scale (left); zooming into the 0%-10% range on the y axis to view less-prominent RNA-type representations (right).  Figure 2, A and B. E, smRNA-TF regulatory networks at G-to-M and M-to-S. miRNAs, tasiRNAs, 21-nt smRNAs, and TFs are shown as triangular, diamond, arrowhead, and circular nodes, respectively. The color of a node represents a cluster to which a smRNA or TF belongs. Black edges indicate the predicted smRNA-target interactions.
is a degenerative process leading to the end of the lifespan.
Our transcriptome analyses further enabled us to identify multilayered regulatory programs that are utilized to temporally coordinate gene expression during the leaf lifespan (Figs. 2 and 5), including transcriptional regulation by dynamic expression of TFs and posttranscriptional regulation by various types of small ncRNAs. In senescing leaves, TFmediated transcriptional regulatory programs were more actively involved in the regulation of age-associated biological processes (Fig. 2E). Various small ncRNAs, including miRNAs and 21-nt smRNAs, and their targets, such as TF (Fig. 5), also showed dynamic changes, suggesting that RNA-level regulation forms complex regulatory networks in coordinating various developmental processes, including senescence, during the leaf lifespan. These findings strongly suggest that plants actively utilize diverse molecular programs coordinating the disassembly and degeneration processes.
A highly resolved and multidimensional transcriptome map generated from our RNA-seq data reveals distinctive structural features of the Arabidopsis leaf transcriptome. Substantial proportions of genes in the chloroplast and mitochondrial genomes, but only 6.9% of nuclear genes, generated antisense transcripts (Fig. 3C). This finding indicates that the organelles, which are derived from prokaryotic ancestors, possess distinct mechanisms for the generation of antisense transcripts. Notably, nuclear genes with antisense transcripts preferentially encode chloroplast-targeted proteins (Fig. 3D), implying the involvement of antisense transcripts in the regulation of chloroplast function.
Changes in transcriptional activity in the chloroplast during the leaf lifespan are fundamental features in the Arabidopsis leaf (Zoschke et al., 2007;Breeze et al., 2011). However, no transcriptome resource currently allows the investigation of agedependent changes in the transcriptional programs of chloroplast genes. Our findings that one-third of the leaf transcriptome originated in the chloroplast genome and that the chloroplast transcriptome showed the most dramatic change during the leaf lifespan (Fig. 6, A-D) provide strong clues supporting that chloroplast transcriptome is a key constituent in the leaf lifespan programs. Interestingly, transcriptional programs for chloroplast function and their regulatory units were tightly coordinated during the leaf lifespan (Fig. 6E). This observation strongly supports the proposition that the cellular degeneration process during leaf aging begins with the decay of chloroplasts, but that the nucleus and mitochondria remain intact until the final stages of the leaf lifespan (Lim et al., 2007). This is an intriguing finding, because it is different from what was observed in animal aging processes that is tightly associated with loss of the mitochondrial function (Finley and Haigis, 2009;Gomes et al., 2013). We also provided evidence supporting the idea that chloroplast functionality is determined by tight coordination between the nuclear and chloroplast transcriptomes. Genes whose corresponding proteins function as protein complexes in chloroplasts, but that are encoded in both nuclear and chloroplast genomes, showed shared age-dependent expression changes during the leaf lifespan (Fig. 6,  F-H). This finding implies that plants have evolved a system for transcriptional coordination as a means of interorganellar communication between the nucleus and chloroplasts and utilize coordinated transcriptional programs to optimize photosynthetic activity and translation activity during the leaf lifespan.
The high-resolution expression dynamics of multiple types of RNAs in the Arabidopsis leaf transcriptome data set lead to hypotheses that can extend current models of age-dependent transitions of key processes (such as photosynthesis, carbon and nitrogen metabolism, and hormonal regulation) during the leaf lifespan. The functional roles of protein-coding RNAs (TFs and differential splicing variants) and regulatory smRNAs (miRNAs, tasiRNAs, and 21-nt smRNAs) with age-dependent expression changes in the initiation and progression of senescence can be investigated by detailed functional experiments, as previously demonstrated for miR164, that can extend the current models of senescence networks Kim 2014a). Likewise, multidimensional dynamic transcriptomes can be used to extend the current understanding of biological networks underlying transitions involved in other key developmental processes.
Various aspects of the comprehensive resource in the Arabidopsis leaf transcriptome data can be analyzed to address a broad spectrum of questions on leaf lifespan. For example, leaf aging is the metabolic turning point from nutrient assimilation to nutrient remobilization. Thus, a thorough understanding of biological networks that delineates leaf aging is essential to the development of strategies to maximize the production of various types of biomass. Other global datasets, such as protein and genetic interactions or protein localizations, can be generated and then integrated with our multidimensional dynamic transcriptomes. This effort would provide spatiotemporal networks that further extend the understanding of the integrated operation of ageassociated networks in space and time. Finally, not just plants, but virtually all organisms undergo a series of developmental changes throughout their lifespan, ending in death. Life history, senescence, and death are fundamental and even philosophical matters for human beings. This study, through the creation of new paradigms in understanding the molecular bases of life history and senescence, may provide a fundamental breakthrough in perceiving the age-dependent senescence and death processes occurring in nature.

CONCLUSION
The multidimensional transcriptome dataset of the Arabidopsis leaf allow us to discern dynamics of transcriptional programs with age-dependent utilization of diverse types of RNAs from different organellar genomes along the leaf lifespan. In-depth analysis of the transcriptome data provides a cornerstone to understand fundamental transcriptional programs underlying age-dependent developmental decisions and functional and regulatory shifts during the lifespan of Arabidopsis leaves.
We conclude that, in contrast to animal aging, leaf senescence proceeds with tight temporal and distinctive interorganellar coordination of various transcriptomes, which has been evolved to maximize nutrient recycling from old tissues to young or newly developing tissues.

Preparation of Arabidopsis Leaves
Fourteen developmental stages of the Arabidopsis (Arabidopsis thaliana) fourth rosette leaf were chosen to explore the leaf transcriptome throughout its lifespan. The fourth rosette leaves of wild-type Columbia plants grown under long-day conditions (16 h light: 8 h dark) were harvested at 4 h after light-on at 2-d intervals from 4 d after emergence to 30 d. Fourth rosette leaves reached a diameter of approximately 3 to 5 mm at 4 d and were fully grown at 16 d; they started yellowing from the tip of the leaf at 20 d and lost approximately 50% of their chlorophyll content at 28 d.

Library Construction and RNA-seq
For strand-specific total RNA-seq, total RNA was extracted from the Arabidopsis leaves using WelPrep (WelGENE, Korea). Total RNA integrity was confirmed using the Agilent Technologies 2100 BioAnalyzer (RNA integrity no. . 8.5). Strand-specific cDNA libraries were generated from total RNA (150 ng) without mRNA selection using the Illumina TruSeq Stranded Total RNA Library Prep Kit with TruSeq Universal Adapter and indexed TruSeq adapters, and Ribo-Minus (Life Technologies) treatment was performed for rRNA depletion. Two independent sample sets were used as biological replicates, and the generated library was subjected to paired-end 101-nt sequencing using an Illumina Hiseq 2000 system. In the case of strand-specific smRNA-seq, total RNA was extracted from seven developmental stages of Arabidopsis leaf using WelPrep. The small RNA library was prepared from a total RNA (1 mg) sample ligated with Illumina RNA 59 adapters and Illumina RNA 39 adapters using the Illumina small RNA preparation kit protocol, according to the manufacturer's instructions. Two independent sample sets were used as biological replicates, and the generated library was then subjected to single-end 51-nt sequencing using an Illumina Hiseq 2000 system.

Alignment of Sequenced Reads to the Arabidopsis Genome
Sequenced reads generated from the strand-specific total RNA-seq were initially submitted to FastQC (Andrews, 2010) to check the quality of raw sequences. The adapter sequences were then trimmed from raw sequences using cutadapt software (Martin, 2011). These trimmed reads were aligned to the Arabidopsis genome (TAIR10; Lamesch et al., 2012) using TopHat  with two mismatches and up to 10 multiple alignments allowed. For strand-specific smRNA-seq, sequenced reads were also submitted to FastQC to check the quality of raw sequences. We then trimmed the adapter sequence from the reads produced from the smRNA-seq. SmRNA-seq reads, with the adaptor sequence trimmed, were retained if they were longer than 17 nt in length. Finally, the remaining reads were aligned to the Arabidopsis genome (TAIR10) using Bowtie (Langmead et al., 2009) with no mismatches and all possible multiple alignments allowed. The alignment results of the strandspecific total RNA-seq and strand-specific smRNA-seq data are summarized in Supplemental Table S1.

Transcript Quantification from RNA-seq Data
We first assembled the aligned reads to transcripts and then used Cufflinks to calculate the fragments per kilobase per million mapped reads (FPKMs) values of the assembled transcripts from the strand-specific total RNA-seq data (Trapnell et al., 2010). Gene annotations were obtained from TAIR10, and we further added up-to-date annotation information on the regions of two 25S rRNAs. Two biological replicates were independently quantitated and then averaged at each time point. Note that FPKM values from only one of the two replicates at 4 d were used because of the poor data quality of the other replicate. From the smRNA-seq alignment results, for each smRNA species, we identified the reads that were aligned to the annotated regions of smRNAs, based on the annotation file obtained from TAIR10. The reads per million mapped reads (RPMs) value was then calculated using the count of the identified reads as the relative abundance of each smRNA. In the case of miRNAs, we used the annotation information for both premiRNAs and mature miRNAs in miRBase20 (Kozomara and Griffiths-Jones, 2011), because TAIR10 provides only premiRNA annotation. Using this annotation information, we calculated the read counts and then RPMs for all the annotated premiRNAs and mature miRNAs separately. Two biological replicates were independently quantitated and then averaged at each time point. As mentioned above, RPM values from only one of the two replicates at 4 d were used because of the poor data quality of the other replicate. We determined a cutoff FPKM value for strand-specific total RNA-seq reads by calculating the number of reads that could cover 50% of the total RNA transcript length. For smRNA transcriptome, we determined that a smRNA species was expressed when the smRNA sequence included more than four smRNA reads at two or more time points among the seven time points.

Identification and Clustering of Differentially Expressed Genes
FPKM values for all RNA transcripts were transformed to log 2 -values after adding 1 to avoid negative infinity, and the log 2 (FPKM) values were smoothed, as done in several previous studies (Del Sorbo et al., 2013;Vardi et al., 2014), before further analyses. The smoothing of time-course gene expression data enables us to clearly analyze temporal expression patterns in clustering and also provides stable statistical distributions based on which we can reliably select the transcripts showing significant correlations between two time-course gene expression profiles and DETs from time-course gene expression data. Among various smoothing methods, the Savitzky-Golay method (with polynomial order = 2 and frame size = 5) was employed to have the same number of time points after smoothing, and thus the smoothed data at a time point can be directly mapped the original data in the corresponding time point (Savitzky and Golay, 1964). To focus on the transcripts showing reproducible temporal gene expression profiles, the transcripts showing significant correlation between the two replicates of time-course gene expression profiles (Wieland et al., 2004;Tsai et al., 2007) was used in the following analyses. To determine a cutoff for the correlation coefficient used to select the transcripts showing reproducible temporal gene expression profiles, we performed the following statistical procedure: (1) for each set of time-course gene expression profiles, we randomly permuted the columns of the n 3 m time-course data matrix (X) for n transcripts and m time points; (2) we calculated the Pearson correlation coefficient for each transcript using the randomized time-course data matrix (X rand ); (3) steps 1 to 2 were repeated 1,000 times; (4) we estimated an empirical distribution for the null hypothesis of the correlation coefficient (a transcript showed no correlation between the two time-course gene expression profiles) using the correlation coefficients obtained from the 1,000 random permutations (Supplemental Fig. S7); and (5) the cutoff was determined as 0.70 to be sufficiently larger than the 95th percentile (0.63) of the null distribution of the correlation coefficient. Log 2 (fold-changes) were calculated by subtracting the log 2 (FPKM) value at 4 d for G-to-M stage and 16 d for M-to-S stage from those at the other time points. Of the transcripts selected using the correlation cutoff, the ones with maximum absolute log 2 (fold-changes) $ 1 were considered to be differentially expressed. The DETs were then clustered based on their age-dependent expression patterns using correlation-based k-means clustering with k = 25. The 25 clusters were further grouped into 'up', 'down', and 'others' groups by performing a hierarchical clustering to their means obtained from the k-means clustering (linkage method = complete, and dissimilarity measure = Euclidean). For clustering of differentially expressed smRNAs, RPM values for all miRNAs and 21-nt long siRNAs including tasiRNAs detected from the smRNA-seq were processed the same as above, but k-means clustering for smRNA reads was performed with k = 10. Finally, the 10 clusters were further grouped into five groups using the hierarchical clustering method as described above.

Enrichment Analysis of GOBPs
GOBPs enriched by the genes in the 12 clusters identified from the clustering of the DETs were identified using the hypergeometric test. For each cluster, of the resulting GOBPs obtained from the enrichment test, we selected the GOBPs that included more than five genes and had enrichment P-values # 0.05. The P-values for all the selected GOBPs in the 12 clusters were presented in an enrichment heat map. In addition, we identified GOBPs enriched by the predicted targets of the small ncRNAs that showed age-dependent expression changes by selecting terms with more than three genes and had the enrichment P-values # 0.05.

Calculation of Degree of Coexpression across GOBPs
To calculate the degree of coexpression between two GOBPs (GOBPs A and B) among 38 GOBPs shown in Figure 2C, a Pearson correlation coefficient for each pair of DETs (one from GOBP A and the other from GOBP B) at G-to-M or M-to-S was computed. The pairs of DETs with correlation coefficients $ 0.99 were then counted. The hypergeometric test was used to evaluate the significance of the number of pairs with correlation coefficients $0.99 as the probability of obtaining the observed number by chance. Significantly coassociated GOBPs were selected as pairs of GOBPs with hypergeometric P-values # 0.05.

Transcriptional Regulatory Circular Map
Protein-DNA interactions (PDIs) of DE-TFs belonging to six clusters (U1-3 and D1-3) at G-to-M or M-to-S stage with their target DETs belonging to 38 GOBPs shown at Figure 2C were visualized using Circos (Krzywinski et al., 2009). For the representation, PDIs confirmed by microarray or qPCR experiments were curated. Six DE-TF clusters and 38 GOBPs were represented by the separate bands (e.g. TF bands for six DE-TF clusters indicated by the dashed arrow in Fig. 2F) on the circumference of Circos. The length of each band is proportional to the count of genes belonging to six DE-TF clusters or 38 GOBPs. Each arrowed link between DE-TF and GOBP bands on the Circos represents that the GOBP included target genes of the TFs in the DE-TF clusters. To evaluate the statistical significance of the number of PDI links between DE-TFs and their target DETs on the circular map at G-to-M or M-to-S stage, permutation tests were performed based on the transcriptional regulatory circular map shown in Figure 2F. For the evaluation, we randomly sampled the same numbers of DE-TFs and DETs at each stage. Then, we counted the number of the curated PDIs between the randomly sampled DE-TFs and DETs. Using the counts from 10 7 random sampling experiments, we estimated the null distribution of the number of PDIs. We then calculated a P-value for the observed number of PDIs on the circular map by one-sided test using the null distribution. For the evaluation of the statistical significance of the number of links of DE-TFs to target DETs included in GOBPs enriched by the DETs at the corresponding stage (G-to-M or M-to-S; Group 1; nongray (red/purple/blue) labeled GOBPs at each stage in Fig. 2F) and also to target DETs included in GOBPs not enriched by the DETs at the corresponding stage (Group 2; gray labeled GOBPs at each stage in Fig. 2F), the same procedure of permutation test was used.

Analysis of Antisense Transcripts
We reversed the strand information in the reference GFF obtained from TAIR10 and generated a new GFF file with the reversed strand information to assemble antisense transcripts. Using this new GFF file, antisense transcripts were assembled using the same alignment result files and quantified as described above, with the reversed strand information. Finally, we determined which antisense transcripts were expressed using the cutoff FPKM values that were used to determine expressed sense transcripts from the strand-specific total RNA-seq data. Determination of DE-ATs and clustering of DE-ATs were performed by the same way with the clustering of total RNA transcripts.

Genome-Wide Visualization of Expressed Antisense Transcripts
The expressed nuclear antisense transcripts, based on TAIR10, were visualized on the five nuclear chromosomes. For reference, we first visualized all the annotated transcripts in TAIR10 in black (intergenic regions in white) over the five individual chromosomes. We then visualized antisense transcripts detected from the strand-specific total RNA-seq in red, within the five chromosomes. The genome-wide visualization of reads aligned to annotated chloroplast transcripts, based on the TAIR10 database, was generated using Circos (Krzywinski et al., 2009). At each base position, we computed the read count by piling up all aligned reads. The read counts greater than the 95th percentile for the entire read counts were set to those of the 95th percentile.

Pile-Up Visualization of Aligned Reads for Selected Chloroplast Genomic Loci
We selected two chloroplast genomic loci as examples that showed agedependent changes in antisense transcripts. For each of the selected loci, piled reads at individual base positions in the sense and antisense directions were separately counted using 'mpileup' in SAMtools (Li et al., 2009). The sense and antisense read counts at individual base positions were plotted at each of the 14 time points along the lifespan, but only the selected 7 time points were shown for conciseness (Supplemental Fig. S2). The structures of the example transcripts (transcript models in gray) are also shown, based on TAIR10 annotation.

GOCC Enrichment Analysis for the Nuclear Genes with Antisense Transcripts
Using BiNGO (Maere et al., 2005), a Cytoscape plugin (Smoot et al., 2011), GOCC enrichment analysis was performed to determine GOCCs enriched by nuclear genes with antisense transcripts. We then selected the GOCCs with the enrichment P-values # 0.05 after Benjamini and Hochberg false discovery rate correction.

Construction of Small ncRNA Regulatory Networks
For the small ncRNAs showing age-dependent expression changes, prediction of their target genes was performed based on the following three target prediction methods using the indicated parameters: (1) psRNA target tool (Dai and Zhao, 2011) using a score #4; (2) TAPIR method (Bonnet et al., 2010) using a score #5 and ratio of the free energy of the duplex to the free energy of the same duplex having only perfect matches $0.7; and (3) Target Finder method (Allen et al., 2005) using a score #5. Among the predicted target candidates, the targets of miRNAs and tasiRNAs were selected as the genes that show agedependent expression changes and have anticorrelation with small ncRNAs (correlation coefficient # -0.7). The same algorithms were used to predict the targets of 21-nt smRNAs. We further filtered 21-nt smRNA target genes by checking the evidence of cleavage based on the previously reported PARE data (GSM1330569, GSM1330570, GSM1330571, GSM1330572, GSM1330573, and GSM1330574 in GEO database). Finally, the targets of 21-nt smRNAs were selected as the ones with matched (mismatch score #4.5) Arabidopsis degradome species from the published PARE data (read count $2) as previously described (Folkes et al., 2012;Thatcher et al., 2015). Pairs of small ncRNA-target genes identified were used to construct small ncRNA-target regulatory networks. Among the small ncRNA-target pairs in Supplemental Table S8, we then selected the pairs with the targets involved in the GOBPs enriched by the targets of miRNAs or 21-nt siRNAs at G-to-M or M-to-S.

Comparison of Degrees of Coexpression between G-to-M and M-to-S Stages
The degree of coexpression for small ncRNAs was evaluated using Pearson correlation coefficients between their expression profiles. To this end, DE-smRNAs were first selected at each stage (G-to-M or M-to-S) as mentioned above. For all possible pairs of DE-smRNAs at each stage, we then computed Pearson correlation coefficients between the expression profiles of individual pairs of DE-smRNAs at the stage. The distributions of the correlation coefficients were represented as boxplots for the two stages, and Wilcoxon rank sum test was performed to determine whether the distributions at the two stages were different from each other. The degree of coexpression for genes encoding photosynthetic complex proteins at G-to-M and M-to-S was also evaluated using the same method.
RP genes localized in the three organelles were first selected. The degrees of coexpression between the selected genes were compared as described above. Pearson correlation coefficients were first calculated for all the pairs of the genes in each group. The distribution of correlation coefficients in each group was then represented using a boxplot as shown in Figure 6F. Finally, one-way ANOVA was performed with Bonferroni correction as a post hoc test to determine whether the distributions between the following pairs are different: (1) mitochondria and cytoplasm, (2) cytoplasm and chloroplasts, and (3) mitochondria and chloroplasts.

Accession Numbers
All of these processed sequence data have been deposited in NCBI's Gene Expression Omnibus (Edgar et al., 2002) with the GEO series accession number of GSE43616. Raw data can be found in NCBI's Sequence Read Archive (Leinonen et al., 2011) with the SRA study accession SRP018034.

Supplemental Data
The following supplemental materials are available.
Supplemental Figure S1. Cellular compartment enrichment analysis for the nuclear genes with antisense transcripts.
Supplemental Figure S2. Visualization of chloroplast loci that express antisense transcripts.
Supplemental Figure S3. GOBPs enriched by identified targets of small ncRNAs.
Supplemental Figure S4. Comparison of 21-nt smRNA regulatory network between G-to-M and M-to-S.
Supplemental Figure S5. Expression of transcripts encoding metabolism related chloroplast-targeted proteins along the leaf lifespan.
Supplemental Figure S6. Coordinated expression changes of transcripts encoding photosynthetic complex proteins along leaf lifespan.
Supplemental Figure S7. The null distribution of correlation coefficients between two replicates of time-course gene expression profiles.
Supplemental Table S1. Mapping result of RNA-seq reads.
Supplemental Table S2. Number of detected transcripts and differentially expressed transcripts.
Supplemental Table S3. List of differentially expressed transcripts.
Supplemental Table S4. Degree of coexpression between all pairs of GOBPs.
Supplemental Table S5. TF-target interactions of transcriptional regulatory circular map.
Supplemental Table S6. List of the genes expressing antisense transcripts.
Supplemental Table S7. List of differentially expressed small RNAs.
Supplemental Table S8. List of predicted smRNA-target pairs.