Physiological Characterization and Transcriptome Analysis of Camellia oleifera Abel. during Leaf Senescence

: Camellia (C.) oleifera Abel. is an evergreen small arbor with high economic value for producing edible oil that is well known for its high level of unsaturated fatty acids. The yield formation of tea oil extracted from fruit originates from the leaves, so leaf senescence, the ﬁnal stage of leaf development, is an important agronomic trait a ﬀ ecting the production and quality of tea oil. However, the physiological characteristics and molecular mechanism underlying leaf senescence of C . oleifera are poorly understood. In this study, we performed physiological observation and de novo transcriptome assembly for annual leaves and biennial leaves of C . oleifera . The physiological assays showed that the content of chlorophyll (Chl), soluble protein, and antioxidant enzymes including superoxide dismutase, peroxide dismutase, and catalase in senescing leaves decreased signiﬁcantly, while the proline and malondialdehyde concentration increased. By analyzing RNA-Seq data, we identiﬁed 4645 signiﬁcantly di ﬀ erentially expressed unigenes (DEGs) in biennial leaves with most associated with ﬂavonoid and phenylpropanoid biosynthesis and phenylalanine metabolism pathways. Among these DEGs, 77 senescence-associated genes (SAGs) including NOL , ATAF1 , MDAR , and SAG12 were classiﬁed to be related to Chl degradation, plant hormone, and oxidation pathways. The further analysis of the 77 SAGs based on the Spearman correlation algorithm showed that there was a signiﬁcant expression correlation between these SAGs, suggesting the potential connections between SAGs in jointly regulating leaf senescence. A total of 162 di ﬀ erentially expressed transcription factors (TFs) identiﬁed during leaf senescence were mostly distributed in MYB (myeloblastosis), ERF (Ethylene-responsive factor), WRKY, and NAC (NAM, ATAF1 / 2 and CUCU2) families. In addition, qRT-PCR analysis of 19 putative SAGs were in accordance with the RNA-Seq data, further conﬁrming the reliability and accuracy of the RNA-Seq. Collectively, we provide the ﬁrst report of the transcriptome analysis of C . oleifera leaves of two kinds of age and a basis for understanding the molecular mechanism of leaf senescence.


Introduction
Camellia (C.) oleifera Abel. is an evergreen small arbor within the genus Camellia of the family Theaceae. It has been cultivated and utilized in China for more than 2000 years and is one of the four major woody oil plants in the world [1]. The tea oil extracted from C. oleifera is mainly composed and analyzed by means of multiple bioinformatics tools. Our study provides an initial analysis for understanding the underlying molecular mechanisms of leaf senescence and is a valuable resource for further identification of genes related to leaf senescence in C. oleifera.

Plant Materials and Growth Conditions
The plant material used in this study was 6-year-old C. oleifera 'Cenxiruanzhi No. 3 . Trees were originally grown in the Fengsheng seedling field in Cenxi City, Guangxi Zhuang Autonomous Region, and transplanted into the greenhouse at Sanqingyuan, Beijing Forestry University, Beijing, with the maximum air temperature of 28 • C in the daytime and a minimum of 21 • C at night, and the humidity of 72%-85%. All plants were watered and fertilized under the same conditions. Seedlings displaying uniform growth and with no signs of disease or insects were selected for further research. The intact annual leaves (AL) and biennial leaves (BL) in the same position on the sunny side were collected in August 2019. All samples were immediately frozen in liquid nitrogen and stored at −80 • C for further analysis. More than three seedlings were pooled as one biological replicate, and three biological replicates were included for our collection. The samples were marked with serial numbers: AL_1, AL_2, AL_3, BL_1, BL_2, and BL_3.

Physiological Analysis of C. oleifera Leaves
The physiological changes differing between the two leaf types C. oleifera were determined before RNA-Seq. The total Chl content was determined by the ethanol-acetone method, with minor modifications [29]. The soluble protein content was measured by using the Coomassie brilliant blue G250 dye method [30]. Malondialdehyde (MDA) concentration was measured by the thiobarbituric acid method, as proposed by a previous study [31]. The proline content was determined according to the sulfosalicylic acid-acid ninhydrin method, with slight modifications [32]. The activities of superoxide dismutase (SOD) was measured by the Total Superoxide Dismutase Assay Kit (hydroxylamine method) and the activity of peroxidase (POD) was determined by the Peroxidase Assay Kit. The activity of catalase (CAT) was measured by the CATalase Assay Kit (Visible light). These assay kits were provided by the Jiancheng Bioengineering Institute (Nanjing, China).

Total RNA Extraction
C. oleifera leaves were ground into powder in liquid nitrogen with a sterilized mortar and pestle for RNA extraction. The total RNA was extracted from the leaves using FastPure Plant Total RNA Isolation Kits (polysaccharides and polyphenolics-rich) (Vazyme, Nanjing, China) according to the manufacturer's instructions. Agilent 2100 and NanoDrop were used to quantify and evaluate the purified RNA. Each assay included three biological replicates.

Construction of the cDNA Library and RNA Sequencing
The extracted RNA was sent to the Beijing Genomics Institute (BGI) (Shenzhen, China), where the library was constructed and sequenced. The mRNA was purified using magnetic beads attached to Oligo (dT) and then purified mRNA was fragmented into small pieces. The first-strand cDNA was synthesized by random hexamer-primed reverse transcription, and subsequently, second strand cDNA was generated. Afterward, the A-tail adaptor was added through incubation, the cDNA fragments obtained after end repair were amplified by PCR, and then were purified by Ampure XP Beads. The purified product was tested on an Agilent Technologies 2100 bioanalyzer. Subsequently, the double-stranded PCR products were denatured and circularized to obtain the final library. The constructed library was finally sequenced on the BGIseq500 platform (BGI, Shenzhen, China).

Analysis and Enrichment of Differentially Expressed Genes (DEGs)
The data were analyzed on the free online platform of Majorbio Cloud Platform (www.majorbio.com). RSEM (RNA-Seq by Expectation Maximization) was utilized to calculate gene expression, and TPM (Transcripts Per Million reads) was used to quantify the expression level of genes [43]. DESwq2 (http://bioconductor.org/packages/stats/bioc/DESeq2/) [44] was used for statistical analysis of raw counts, and the DEGs were defined with p-adjust < 0.05 and |log2FC| ≥ 1, p-adjust is the p-value after multiple test correction using BH (fdr correction with Benjamini/Hochberg). GO enrichment analysis of DEGs were performed by using Goatools (https//github.com/tanghaibao/GOatools) [45]. When the corrected p-value < 0.05, the DEGs were defined to be significantly enriched. The KEGG pathway enrichment analysis of DEGs was performed by KOBAS (http://kobas.cbi.pku.edu.cn/home.do) [46]. The DEGs were considered to be significantly enriched when the corrected p-value < 0.05.

Screening of Key Senescence-Associated Genes (SAGs)
We initially screened 77 SAGs on the Majorbio Cloud Platform. By using the BLAST (Basic Local Alignment Search Tool) program of the NCBI (National Center for Biotechnology Information) and LSD 3.0 (Leaf Senescence Database, https://bigd.big.ac.cn/lsd/index.php) [47], we selected some SAGs with high homology to important SAGs in other species for further analysis.

Analysis of Expression Correlation
The Spearman algorithm was used to obtain the correlation coefficient between genes, and the expression correlation was defined with a correlation coefficient ≥ 0.5 and p-adjust < 0.05 [48]. According to the expression correlation between genes, the visual network diagram was constructed.

Quantitative Real-Time PCR (qRT-PCR) Analysis of Gene Expression
To evaluate the quality of RNA-Seq data, 19 unigenes were selected for qRT-PCR analysis. The extracted RNA was used to synthesize first-strand cDNAs using a FastQuant cDNA First-Strand Synthesis Kit (Tiangen Biotechnology, Beijing, China). We used the StepOne Plus Real-Time PCR System (ABI, Vernon, CA, USA) for qRT-PCR, according to the manufacturer's instructions of 2 × Fast SYBR Green qPCR Master Mix (High ROX) (Servicebio, Wuhan, China). GAPDH was used as a reference gene. The primer sequences of selected genes are shown in Table S4. All qRT-PCR assays were performed for three biological replicates and three technical repeats.

Morphological and Physiological Observation of Leaves at Different Leaf Age
We observed the phenotype and measured physiological changes of two different kinds of C. oleifera leaves in the rapid fruit growth period. Our results showed that biennial leaves were significantly yellower than annual leaves ( Figure 1A). The visible morphological changes we observed are consistent with previous reports [49] where yellowing due to Chl degradation is one of the most obvious characteristics of leaf senescence. Accordingly, in contrast to annual leaves, the Chl content in biennial leaves was significantly reduced (by 34.6%), suggesting that the degradation rate of Chl in biennial C. oleifera leaves was higher than the synthesis rate ( Figure 1B). Additionally, the soluble protein content of biennial leaves was only 48.5% of that in annual leaves ( Figure 1C), suggesting that the hydrolysis of soluble protein dominates during the senescence process of C. oleifera. In addition, we also examined the accumulation of MDA and proline and found that the contents of MDA and proline in biennial leaves were 1.48-fold higher and 2.57-fold higher than that in annual leaves, respectively ( Figure 1D,E). Moreover, compared with annual leaves, the activity of SOD, POD, and CAT in senescent leaves decreased by 9.8%, 44.5%, and 49.8% ( Figure 1F-H), respectively, indicating that the activities of antioxidant enzymes were reduced in C. oleifera biennial leaves, thus repressing the scavenging ability of ROS (reactive oxygen species). observed are consistent with previous reports [49] where yellowing due to Chl degradation is one of the most obvious characteristics of leaf senescence. Accordingly, in contrast to annual leaves, the Chl content in biennial leaves was significantly reduced (by 34.6%), suggesting that the degradation rate of Chl in biennial C. oleifera leaves was higher than the synthesis rate ( Figure 1B). Additionally, the soluble protein content of biennial leaves was only 48.5% of that in annual leaves ( Figure 1C), suggesting that the hydrolysis of soluble protein dominates during the senescence process of C. oleifera. In addition, we also examined the accumulation of MDA and proline and found that the contents of MDA and proline in biennial leaves were 1.48-fold higher and 2.57-fold higher than that in annual leaves, respectively ( Figure 1D, E). Moreover, compared with annual leaves, the activity of SOD, POD, and CAT in senescent leaves decreased by 9.8%, 44.5%, and 49.8% ( Figure 1F-H), respectively, indicating that the activities of antioxidant enzymes were reduced in C. oleifera biennial leaves, thus repressing the scavenging ability of ROS (reactive oxygen species).

RNA Sequencing, De Novo Assembly, and Annotation
In order to explore the molecular mechanisms underlying the leaf senescence of C. oleifera at the transcriptomic level, six RNA libraries from C. oleifera leaves (AL1, AL2, AL3; BL1, BL2, and BL3) were sequenced and 37.16 Gb high-quality clean data were obtained. The clean data of each sample reached more than 6.08 Gb, and the clean reads of each sample ranged from 40,602,122 to 41,859,058 with the Q30 percentage over 90.32% (Table S1). After de novo assembly of the C. oleifera transcriptome, a total of 78,860 unigenes were generated with an N50 length of 1374 bp and an average length of 865.96 bp. The GC content was approximately 39.19% (Table 1).
The assembled unigenes were searched against the NR (NCBI non-redundant protein sequences), Swiss-Prot (Annotated protein sequence

RNA Sequencing, De Novo Assembly, and Annotation
In order to explore the molecular mechanisms underlying the leaf senescence of C. oleifera at the transcriptomic level, six RNA libraries from C. oleifera leaves (AL1, AL2, AL3; BL1, BL2, and BL3) were sequenced and 37.16 Gb high-quality clean data were obtained. The clean data of each sample reached more than 6.08 Gb, and the clean reads of each sample ranged from 40,602,122 to 41,859,058 with the Q30 percentage over 90.32% (Table S1). After de novo assembly of the C. oleifera transcriptome, a total of 78,860 unigenes were generated with an N50 length of 1374 bp and an average length of 865.96 bp. The GC content was approximately 39.19% (Table 1). The assembled unigenes were searched against the NR (NCBI non-redundant protein sequences), Swiss-Prot (Annotated protein sequence database), Pfam (Protein family), COG (Clusters of Genes), GO (Gene Ontology), and KEGG (Kyoto Encyclopedia of Genes and Genomes) databases, with 51.3% (40,425 of 78,860) of unigenes annotated. Among them, 26,486, 13,831, 29,317, 39,848, 24,214, and 23,566 annotated unigenes were obtained from the GO, KEGG, COG, NR, Swiss-Prot, and Pfam databases, respectively ( Figure 2). The results indicated that the NR database had the highest proportion of annotation and 39,848 unigenes were annotated. Due to the lack of a reference genome for C. oleifera, the NR database can provide us with information about the most similar matches to gene sequences from other species, and according to the NR database, the species with the greatest number of C. oleifera unigenes are Camellia sinensis (87.70%), followed by Actinidia chinensis (2.02%) and Vitis vinifera (1.22%) ( Figure S1). These results provide a basis for the future research of the potential function of these annotated unigenes. and Vitis vinifera (1.22%) ( Figure S1). These results provide a basis for the future research of the potential function of these annotated unigenes.

Analyses of Gene Expression and Identification of Differentially Expressed Genes (DEGs)
The TPM values were used to calculate the expression level of unigenes in each tissue, and unigenes with TPM ≥1 were defined as expressed [50]. The number of unigense expressed in annual leaves was 59,212, and that in biennial leaves was 47,700. In addition, 40,421 unigenes were co-

Analyses of Gene Expression and Identification of Differentially Expressed Genes (DEGs)
The TPM values were used to calculate the expression level of unigenes in each tissue, and unigenes with TPM ≥1 were defined as expressed [50]. The number of unigense expressed in annual leaves was 59,212, and that in biennial leaves was 47,700. In addition, 40,421 unigenes were co-expressed in both leaves ( Figure S2). Moreover, correlation analyses were performed to verify the consistency among the samples, and our results demonstrated that the correlation among three biological replicates was very high ( Figure S3). The DEGs were defined with the adjusted p-value < 0.05 and fold-change ≥ 2, and our results showed that there were 4645 DEGs between C. oleifera leaves of two types of ages. Taking the annual leaves as a control group, the number of upregulated genes and downregulated genes in biennial leaves were 2729 and 1916, respectively ( Figure 3A). To facilitate the visualization, a heatmap was produced on the basis of the clustering analysis of expression patterns for all DEGs ( Figure 3B). These findings demonstrate that a large number of genes are responsive to the senescence process. expression patterns for all DEGs ( Figure 3B). These findings demonstrate that a large number of genes are responsive to the senescence process.

Gene Ontology (GO) Enrichment Analysis of DEGs
To further elucidate the function of DEGs, a GO enrichment analysis was performed. In terms of all DEGs, there were 298 enriched GO terms (p-value ≤ 0.05). The most significantly enriched GO category was the extracellular region (GO ID: 0005576), followed by the intrinsic component of the membrane (GO ID: 0031224), nucleic acid binding transcription (GO ID: 0001071), transcription factor activity, sequence-specific DNA binding (GO ID: 0003700), and oxidoreductase activity (GO ID: 0016491) ( Figure 4A). In terms of upregulated DEGs, there were 292 significantly enriched GO terms including the extracellular region (GO ID: 0005576), dioxygenase activity (GO ID: 0051213), membrane part (GO ID: 0044425), etc ( Figure 4B). As for downregulated DEGs, there were 219 terms such as beta-amyrin synthase activity (GO ID: 0042300), oxidosqualene cyclase activity (GO ID: 0031559), and lanosterol synthase activity (GO ID: 0000250) ( Figure 4C). Each point represents an expressed unigene, the abscissa represents the fold change value of gene expression difference in the two samples, the vertical axis represents the statistical test value of the difference in gene expression, p-value. Compared with annual leaves, red shows that the gene expression was upregulated in biennial leaves, green represents the gene expression that was downregulated in biennial leaves, and grey represents non-differentially expressed genes. (B) Heat map of differentially expressed genes. AL_1, AL_2 and AL_3 represent three biological replicates of annual leaves. BL_1, BL_2, and BL_3 represent three biological replicates of biennial leaves.

Gene Ontology (GO) Enrichment Analysis of DEGs
To further elucidate the function of DEGs, a GO enrichment analysis was performed. In terms of all DEGs, there were 298 enriched GO terms (p-value ≤ 0.05). The most significantly enriched GO category was the extracellular region (GO ID: 0005576), followed by the intrinsic component of the membrane (GO ID: 0031224), nucleic acid binding transcription (GO ID: 0001071), transcription factor activity, sequence-specific DNA binding (GO ID: 0003700), and oxidoreductase activity (GO ID: 0016491) ( Figure 4A). In terms of upregulated DEGs, there were 292 significantly enriched GO terms including the extracellular region (GO ID: 0005576), dioxygenase activity (GO ID: 0051213), membrane part (GO ID: 0044425), etc ( Figure 4B). As for downregulated DEGs, there were 219 terms such as beta-amyrin synthase activity (GO ID: 0042300), oxidosqualene cyclase activity (GO ID: 0031559), and lanosterol synthase activity (GO ID: 0000250) ( Figure 4C).

Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway Analysis of DEGs
To further discover the leaf senescence pathway of C. oleifera leaves, KEGG analysis of DEGs was performed. A total of 126 KEGG pathways were significantly enriched, among which the phenylpropanoid biosynthesis (Pathway id: map00940), plant hormone signal transduction (Pathway id: map04075), and phenylalanine metabolism (Pathway id: map00360) pathways were the most obviously enriched pathways ( Figure 5A). The upregulated DEGs were assigned to the KEGG pathways involved in 125 pathways, among them, phenylpropanoid biosynthesis (Pathway id: map00940), phenylalanine metabolism (Pathway id: map00360), and arginine and proline metabolism (Pathway id: map00330) were the most highly represented ( Figure 5B). The downregulated DEGs were assigned to 90 KEGG pathways, among which the sesquiterpenoid and triterpenoid biosynthesis (Pathway id: map00909), plant hormone signal transduction (Pathway id: map04075), and ascorbate and aldarate metabolism (Pathway id: map00053) pathways were the most significantly enriched pathways ( Figure 5C).

SAGs Differentially Expressed in Two Kinds of Leaves
To identify the important gene expressions of C. oleifera leaves during leaf senescence, we analyzed DEGs associated with Chl degradation, plant hormones, oxidation pathways, and found 77 unigenes homologous to the key SAGs that have been reported in other species (Table S2). Many

SAGs Differentially Expressed in Two Kinds of Leaves
To identify the important gene expressions of C. oleifera leaves during leaf senescence, we analyzed DEGs associated with Chl degradation, plant hormones, oxidation pathways, and found 77 unigenes homologous to the key SAGs that have been reported in other species (Table S2). Many genes involved in the Chl degradation pathway were upregulated in biennial leaves including PAO, NOL, SGR, etc. ( Figure 6A, Figure S4A). In addition, genes related to ABA, ethylene, JA, and SA were mostly upregulated in biennial leaves such as MYBL, ATAF1, and PYL9 ( Figure S4B). Genes associated with auxin and cytokinin were also differentially expressed during leaf senescence. For example, the gene homologous to cytokinin receptor gene AHK2 in arabidopsis was downregulated, and the gene homologous to auxin-related gene SAUR36 was upregulated. During foliar senescence, evident changes for oxidation-related genes have been observed, for instance, WRKY75, encoding a transcription factor inhibiting the elimination of ROS, was upregulated, whereas MDAR, encoding a monodehydroascorbate reductase, was downregulated ( Figure 6C, Figure S4C). Furthermore, we found some other SAGs were also responsive to the senescence of C. oleifera leaves such as SAG12, which was evidently upregulated in biennial leaves ( Figure 6D). These findings suggest that there are multiple signaling pathways involved in the regulation of C. oleifera leaf senescence. genes involved in the Chl degradation pathway were upregulated in biennial leaves including PAO, NOL, SGR, etc. ( Figure 6A, Figure S4A). In addition, genes related to ABA, ethylene, JA, and SA were mostly upregulated in biennial leaves such as MYBL, ATAF1, and PYL9 ( Figure S4B). Genes associated with auxin and cytokinin were also differentially expressed during leaf senescence. For example, the gene homologous to cytokinin receptor gene AHK2 in arabidopsis was downregulated, and the gene homologous to auxin-related gene SAUR36 was upregulated. During foliar senescence, evident changes for oxidation-related genes have been observed, for instance, WRKY75, encoding a transcription factor inhibiting the elimination of ROS, was upregulated, whereas MDAR, encoding a monodehydroascorbate reductase, was downregulated ( Figure 6C, Figure S4C). Furthermore, we found some other SAGs were also responsive to the senescence of C. oleifera leaves such as SAG12, which was evidently upregulated in biennial leaves ( Figure 6D). These findings suggest that there are multiple signaling pathways involved in the regulation of C. oleifera leaf senescence.

Correlation Analysis of SAGs
Previous studies have shown that leaf senescence is a complicated physiological process regulated by multiple signal pathways. To better reveal the potential co-expression relationship between genes, we constructed an expression correlation network of the 77 SAGs on the basis of the Spearman correlation algorithm. The results indicated that there was a possible correlation among genes involved in different pathways such as Chl degradation, phytohormones, and oxidation pathways. Among these SAGs, CoWRKY75-2, CoPAO-2, CoSAG12, CoWRKY65-2, CoMYBL-1,

Correlation Analysis of SAGs
Previous studies have shown that leaf senescence is a complicated physiological process regulated by multiple signal pathways. To better reveal the potential co-expression relationship between genes, we constructed an expression correlation network of the 77 SAGs on the basis of the Spearman correlation algorithm. The results indicated that there was a possible correlation among genes involved in different pathways such as Chl degradation, phytohormones, and oxidation pathways. Among these SAGs, CoWRKY75-2, CoPAO-2, CoSAG12, CoWRKY65-2, CoMYBL-1, CoSAUR36-1, CoNOL-1, CoNAC59, CoGBF1, and CoS40-2 were significantly correlated with the expression of other SAGs (Figure 7). For instance, there is a significant expression correlation between CoSAG12 and many other SAGs such as CoNOL-2, CoORE1, CoMYBL-1, CoS40-2, and CoWRKY75-2, and these genes probably play important roles in Chl degradation, hormones, and oxidation pathways. CoSAUR36-1, CoNOL-1, CoNAC59, CoGBF1, and CoS40-2 were significantly correlated with the expression of other SAGs (Figure 7). For instance, there is a significant expression correlation between CoSAG12 and many other SAGs such as CoNOL-2, CoORE1, CoMYBL-1, CoS40-2, and CoWRKY75-2, and these genes probably play important roles in Chl degradation, hormones, and oxidation pathways. Figure 7. Expression correlation analysis of putative SAGs. Each node represents a unigene, and the connection between nodes represents the expression correlation of genes. Large node suggests that this gene has expression correlation with a large number of other genes.

Transcription Factors (TFs) Responding to Leaf Senescence
Transcription factors (TFs) exert vital function in the process of leaf senescence. Therefore, we analyzed the expression levels of TFs in C. oleifera leaves of two types of ages. Among all the DEGs, we identified 162 TFs (123 upregulated TFs and 39 downregulated TFs) (Table S3). Among these different TF families, the genes encoding MYB proteins (16.67%) accounted for the largest proportion of these genes, followed by the genes encoding AP2/ERF proteins (16.05%), WRKY (12.35%), and NAC proteins (9.3%) ( Figure 8). Notably, many TFs we have identified are homologous to SAGs in other species genome such as MYBL, ERF1, WRKY75, ATAF1, etc. Among these potential SAGs in C. oleifera, the genes that exhibited the most significant changes were selected for further qRT-PCR analysis. Figure 7. Expression correlation analysis of putative SAGs. Each node represents a unigene, and the connection between nodes represents the expression correlation of genes. Large node suggests that this gene has expression correlation with a large number of other genes.

Transcription Factors (TFs) Responding to Leaf Senescence
Transcription factors (TFs) exert vital function in the process of leaf senescence. Therefore, we analyzed the expression levels of TFs in C. oleifera leaves of two types of ages. Among all the DEGs, we identified 162 TFs (123 upregulated TFs and 39 downregulated TFs) (Table S3). Among these different TF families, the genes encoding MYB proteins (16.67%) accounted for the largest proportion of these genes, followed by the genes encoding AP2/ERF proteins (16.05%), WRKY (12.35%), and NAC proteins (9.3%) ( Figure 8). Notably, many TFs we have identified are homologous to SAGs in other species genome such as MYBL, ERF1, WRKY75, ATAF1, etc. Among these potential SAGs in C. oleifera, the genes that exhibited the most significant changes were selected for further qRT-PCR analysis. CoSAUR36-1, CoNOL-1, CoNAC59, CoGBF1, and CoS40-2 were significantly correlated with the expression of other SAGs (Figure 7). For instance, there is a significant expression correlation between CoSAG12 and many other SAGs such as CoNOL-2, CoORE1, CoMYBL-1, CoS40-2, and CoWRKY75-2, and these genes probably play important roles in Chl degradation, hormones, and oxidation pathways. Figure 7. Expression correlation analysis of putative SAGs. Each node represents a unigene, and the connection between nodes represents the expression correlation of genes. Large node suggests that this gene has expression correlation with a large number of other genes.

Transcription Factors (TFs) Responding to Leaf Senescence
Transcription factors (TFs) exert vital function in the process of leaf senescence. Therefore, we analyzed the expression levels of TFs in C. oleifera leaves of two types of ages. Among all the DEGs, we identified 162 TFs (123 upregulated TFs and 39 downregulated TFs) (Table S3). Among these different TF families, the genes encoding MYB proteins (16.67%) accounted for the largest proportion of these genes, followed by the genes encoding AP2/ERF proteins (16.05%), WRKY (12.35%), and NAC proteins (9.3%) ( Figure 8). Notably, many TFs we have identified are homologous to SAGs in other species genome such as MYBL, ERF1, WRKY75, ATAF1, etc. Among these potential SAGs in C. oleifera, the genes that exhibited the most significant changes were selected for further qRT-PCR analysis.

Validation of RNA-Seq Data by qRT-PCR
To verify the accuracy and reliability of the transcriptome data, 19 putative SAGs, which were speculated to be highly correlated with other SAGs by expression correlation analysis, were selected for further analysis using qRT-PCR. Our findings indicated that the results of qRT-PCR analysis were in accordance with the RNA-Seq data. Compared with annual leaves, the expression level of all predicted upregulated SAGs was significantly increased in biennial leaves, and two downregulated SAGs were correspondingly inhibited in senescent leaves (Figure 9). Therefore, these results suggest that our RNA-Seq data were reliable and consistent with the qRT-PCR analysis. Figure 8. Distribution of differentially expressed transcription factors (TFs). Compared with annual leaves, red represents up-regulated TFs in biennial leaves, blue represents down-regulated TFs in biennial leaves.

Validation of RNA-Seq Data by qRT-PCR
To verify the accuracy and reliability of the transcriptome data, 19 putative SAGs, which were speculated to be highly correlated with other SAGs by expression correlation analysis, were selected for further analysis using qRT-PCR. Our findings indicated that the results of qRT-PCR analysis were in accordance with the RNA-Seq data. Compared with annual leaves, the expression level of all predicted upregulated SAGs was significantly increased in biennial leaves, and two downregulated SAGs were correspondingly inhibited in senescent leaves (Figure 9). Therefore, these results suggest that our RNA-Seq data were reliable and consistent with the qRT-PCR analysis.

Discussion
In this study, we demonstrated that some important physiological differences occurred between younger and older leaves of C. oleifera. For example, compared with the younger annual leaves, the older biennial leaves showed less Chl concentration and antioxidant enzyme activities, and more accumulation of MDA. In fact, similar phenomena were also observed in other species. For instance, Gossypium hirsutum displays Chl breakdown and increased production of MDA during leaf senescence [51]. In Triticum aestivum, the activities of antioxidant enzymes such as SOD and CAT were inhibited in the process of leaf senescence, suggesting that the leaf senescence of the plant is generally accompanied by physiological changes such as the degradation of Chl and some macromolecules, and the reduction of antioxidant enzyme activity [52]. Our findings indicated that compared with the annual leaves, the biennial leaves showed obvious senescence phenotypes and

Discussion
In this study, we demonstrated that some important physiological differences occurred between younger and older leaves of C. oleifera. For example, compared with the younger annual leaves, the older biennial leaves showed less Chl concentration and antioxidant enzyme activities, and more accumulation of MDA. In fact, similar phenomena were also observed in other species. For instance, Gossypium hirsutum displays Chl breakdown and increased production of MDA during leaf senescence [51]. In Triticum aestivum, the activities of antioxidant enzymes such as SOD and CAT were inhibited in the process of leaf senescence, suggesting that the leaf senescence of the plant is generally accompanied by physiological changes such as the degradation of Chl and some macromolecules, and the reduction of antioxidant enzyme activity [52]. Our findings indicated that compared with the annual leaves, the biennial leaves showed obvious senescence phenotypes and physiological change, suggesting that these two kinds of leaves are suitable for the study of C. oleifera leaf senescence.
To gain more insight into the mechanism of leaf senescence in C. oleifera, a transcriptomic analysis was performed in both annual leaves and biennial leaves. In our study, a total of 4645 DEGs were identified, and we further interpreted the potential biological functions of DEGs from the gene function and signaling pathway through GO enrichment analysis and KEGG enrichment analysis, respectively. The results of GO enrichment showed that with the duration of leaf senescence of C. oleifera, 'extracellular region', 'intrinsic component of membrane', and 'oxidoreductase activity' were significantly enriched, which were consistent with the physiological changes we measured. During leaf senescence of C. oleifera, the activity of antioxidant enzyme and MDA concentration changed, and MDA is one of the most important products associated with membrane lipid peroxidation, and can further cause damage to cell membranes [53]. Therefore, on one hand, our results suggest that ROS plays a vital role in the process of C. oleifera leaf senescence, and on the other hand, these findings demonstrate that the physiological changes we determined and the analysis of transcriptome are reliable. Furthermore, based on KEGG enrichment, we found that 'phenylpropanoid biosynthesis' and 'plant hormone signal transduction' were significantly enriched. In fact, previous study also showed that during the leaf senescence of Sorghum bicolor, 'phenylpropanoid biosynthesis' and 'plant hormone signal transduction' exhibited the most DEGs [54].
Leaf senescence is a highly complex process, which is controlled by multiple pathways [49]. Based on previous studies and our findings on the physiological changes and transcriptomic analysis of C. oleifera, we speculated that leaf senescence may cause changes in pathways such as Chl degradation, plant hormones, and oxidation, etc. In order to reveal these changes, we performed differential expression statistics on the paramount pathways and found 77 essential SAGs in the process of C. oleifera leaf senescence, many of which are associated with Chl degradation, phytohormone signaling, and oxidation. For instance, compared with annual leaves, CoNOL and other Chl degradation related genes, ABA receptor gene CoPYL9 and senescence-related SAG12 were upregulated in biennial leaves, while cytokinin receptor CoAHK2 and antioxidant enzyme gene CoMDAR-2 were downregulated. Similar results were also observed in other species. For example, in Brassica napus, the SAG12 transcript evidently accumulates in older leaves [17]. In addition, in arabidopsis, during dark-induced senescence, the expression of SGR, NOL, and PAO also significantly increased to degrade Chl [55]. These results indicate that a similar senescence mechanism exists in the process of leaf senescence induced by environmental factor or age.
Notably, in this study, based on the analysis of these SAGs, we found that CoORE1 is simultaneously involved in Chl degradation and plant hormone signaling pathways, suggesting that SAGs may participate in more than one pathway to regulate leaf senescence. In fact, previous study in Brassica rapa demonstrated that BrNAC055 accelerates leaf senescence by activating the transcription of RBOH (ROS-producing enzymes respiratory burse oxidase homologue) and two Chl catabolic genes (BrNYC1 and BrNYE1), indicating that BrNAC055 is involved in Chl degradation and oxidation pathways [56]. To identify the potential correlation between these important pathways, we further constructed an expression correlation network for the 77 SAGs identified and found that there is an expression correlation between genes in different pathways. For example, CoORE1, which is involved in both Chl degradation and plant hormone pathways, correlated with CoWRKY75-2, a WRKY transcription factor that participates in both oxidation and hormone pathways. Furthermore, CoORE1 also correlated with other key SAGs such as CoSAG21-2 and CoSAG12. Expression correlation implies, to a certain extent, potential connections between genes. For instance, there is an expression correlation between senescence-related gene ONAC054 and ABA signaling associated gene OsABI5, and further experimental results indicated that ONAC054 directly upregulates the expression of OsABI5, thus regulating leaf senescence of rice [57]. Therefore, by constructing the expression correlation network, we predicted the possible connection between SAGs from a bioinformatics perspective. Nevertheless, the accurate regulatory relationship needs further experimental verification in future research.
TFs have important functions in transmitting and amplifying leaf senescence signals, thus activating a cascade of downstream genes [9]. Overexpression of Cucumis melo TF CmNAC60 in arabidopsis can promote leaf senescence of transgenic plants [58]. In this study, we performed statistical analysis on TFs and found that the differentially expressed TFs were mostly distributed in the MYB, WRKY, NAC families, etc. Expression correlation analysis showed that the functions of these senescence-associated TFs are complex. For instance, CoWRKY75-2 is highly correlated with the expression of several significant senescence-associated genes such as CoSAG12 and CoSAG21-1.
CoNAC72 not only participates in Chl degradation, but also plays an important role in the hormone pathway. Notably, there is an expression correlation between CoNAC72 and CoATAF1-2. It is reported that homodimerization and heterodimerization among NAC TFs are crucial mechanisms in controlling NAC transcription factor mediated plant developmental processes [59,60], suggesting that these NAC TFs may form homodimers or heterodimers to function in regulating leaf senescence of C. oleifera. This phenomenon can also be found in rice, in which ONAC020 and ONAC026 share similar expression patterns during seed development, and ONAC020 can interact with ONAC026 [61]. In fact, in addition to NAC TFs, there have been similar reports on the study of WRKY TFs. In arabidopsis, WRKY54 and WRKY70 have been identified as important in leaf senescence, and they share a similar expression pattern during leaf senescence and co-operate as negative regulators [62]. Furthermore, in the present study, by using qRT-PCR, we found that the expression levels of key transcription factors such as CoMYBL-1 and CoNAC59 changed significantly, which was in good accordance with the transcriptome data, suggesting that TFs play important roles in C. oleifera leaf senescence based on the transcriptional level. In future research, CoMYBL-1, CoNAC72, CoWRKY75-2, and other TFs can be chosen as important candidates for the potential association between TFs and other SAGs to elucidate the leaf senescence mechanism. In addition to participating in senescence, TFs such as NAC72, ATAF1, and WRKY75 have also been widely reported to play important roles in the process of abiotic stress response in plants such as drought and salinity [63][64][65]. Our physiological results demonstrated that the accumulation of ROS in the annual leaves and biennial leaves was significantly different. Previous reports have shown that it is a common function of TFs to participate in the regulation of plants in response to abiotic stress through ROS pathways. Our qRT-PCR analysis indicated that the stress-responsive genes were differentially expressed between the two leaf types. Therefore, we speculated that beyond senescence, the two types of leaves also have tolerance differences in response to drought and other abiotic stresses, and TFs such as CoNAC72 and CoATAF1 may also be involved in the process of C. oleifera in response to abiotic stress. However, future research should acquire more data to verify this information.

Conclusions
To summarize, our results present a transcriptome analysis that combined physiological data and phenotypic observation, which contributes to the understanding of gene expression profiling underlying leaf senescence of C. oleifera. Compared with the annual leaves, a large number of DEGs were identified in the biennial leaves including 162 differentially expressed TFs. The present study explored 77 putative SAGs, which may be involved in the regulation of leaf senescence through Chl degradation, plant hormones, and oxidation pathways. There was a significant expression correlation between these SAGs, suggesting that these SAGs may jointly regulate leaf senescence. qRT-PCR analysis of 19 putative SAGs were in accordance with the RNA-Seq data, and in future research, these SAGs may serve as candidate genes for delaying leaf senescence in molecular breeding programs, thus increasing the yield of C. oleifera.
Supplementary Materials: The following are available online at http://www.mdpi.com/1999-4907/11/8/812/s1, Figure S1: NR annotated homologous species distribution map, Figure S2: Number of expressed unigenes in two kinds of C. oleifera leaves, Figure S3: Correlation analysis between different samples, Figure S4: Heat map of SAGs associated with important pathways, Table S1: Summary of sequencing data, Table S2: The details of putative SAGs, Table S3: The detail of differentially expressed TFs, Table S4: Primer sequences used in the qRT-PCR experiment.