De novo Characterization of the Platycladus orientalis Transcriptome and Analysis of Photosynthesis-Related Genes during Aging

: In China, Platycladus orientalis has a lifespan of thousands of years. The long lifespan of these trees may be relevant for the characterization of plant aging at the molecular level. However, the molecular mechanism of the aging process of P. orientalis is still unknown. To explore the relationship between age and growth of P. orientalis , we analyzed physiological changes during P. orientalis senescence. The malondialdehyde content was greater in 200-, 700-, and 1100-year-old ancient trees than in 20-year-old trees, whereas the peroxidase and superoxide dismutase activities, as well as the soluble protein content, exhibited the opposite trend. Furthermore, we performed a de novo transcriptome assembly using RNA-Seq and obtained 48,044 unigenes with an average length of 896 bp. A total of 418 di ﬀ erentially expressed genes were identiﬁed in di ﬀ erent stages of aging of P. orientalis . Clustering analysis revealed distinct timepoints at which the oxidation–reduction and photosynthesis pathways changed. Eight clusters with distinct expression patterns were identiﬁed. The expression levels of photosynthesis-, oxidation–reduction-, and transporter-related genes were down-regulated, whereas those of transcription-, signaling-, and senescence-related genes were up-regulated during aging. In addition, consistent with the most obviously down-regulated genes of photosynthesis-related genes, the photosynthetic indexes including chlorophyll a and b levels decreased steadily during P. orientalis aging. This study combined transcriptome with physiological and biochemical data, revealing potential candidate genes inﬂuencing senescence during P. orientalis aging. ,


Introduction
Several woody perennials, including Sequoia sempervirens, Sequoiadendron giganteum, Platycladus orientalis, and Pinus longaeva, can live for several hundred years or even millennia [1,2]. Plant biologists are currently interested in the physiological and molecular mechanisms mediating or preventing senescence [3,4]. Some physiological and metabolic changes during aging have been reported in conifers [4][5][6]. However, no significant age-related differences were observed between 4713-year-old which were used as controls. Additionally, all selected trees were located close to each other and were exposed to similar environmental conditions. Samples were collected as described by Zhang et al. [4]. The leaves were collected at 8:00 am in July of 2016. Specifically, new leaves were selected from the south-facing middle canopy layer of pest-and disease-free trees exposed to similar daily illumination times. Samples were flash frozen in liquid nitrogen.

Physiological Index Determination
The malondialdehyde (MDA) and soluble protein contents were measured as described by Cui et al. and Heath and Packer [27,28], respectively. The SOD and POD activities were detected as described by Dhindsa et al. [29]. The chlorophyll a and b and caretonoid contents were estimated as described by Breeze et al. [30]. Four separate experiments were conducted, differences were considered significant at the p < 0.05 and p < 0.01 levels. Data were analyzed by ANOVA in the SPSS software (IBM, Armonk, NY, USA).

Total RNA Extraction
A Column Plant RNAout kit (TIANDZ, Beijing, China) was used to isolate total RNA from the leaves of 20-, 200-,700-, and 1100-year-old P. orientalis trees, with three biological replicates. The RNA concentrations were determined with the NanoPhotometer ® spectrophotometer (Implen, Westlake Village, CA, USA), whereas RNA quality and integrity were assessed by 1.5% agarose gel electrophoresis and with a 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA).

Construction and Sequencing of cDNA Libraries
Oligo-(dT) magnetic beads were used to purify mRNA from 20 µg total RNA. Fragmentation buffer was added to the purified mRNA to generate shorter segments, which were used along with six-base random primers for first-strand cDNA synthesis. The second cDNA strand was synthesized with DNA polymerase I, dNTPs, and RNase H. Sequencing adapters were ligated to the short segments purified with the QIAquick PCR purification kit to discriminate among samples (three biological replicates). Segments were amplified by PCR and then separated by agarose gel electrophoresis to produce sequencing templates. The resulting cDNA libraries were sequenced with the GAIIx system at Novogene (Beijing, China). Sequencing datasets were deposited in the NCBI database (accession: SRX1757589).

Transcriptome Assembly and Functional Annotation
Raw reads were filtered by removing adapter sequences and low-quality reads (<Q20) [Q = −10log 10 (e), which indicates the base quality; e refers to the sequencing error rate]. Additionally, the Q20, Q30, and GC content of the clean data were calculated, after which the Trinity program was used for the de novo assembly of clean sequences [31]. All of the parameters were set to default values, and the expression level of the assembled transcripts was based on the fragments per kilobase per million mapped read (FPKM) values [32]. The abundance of the assembled transcripts was quantified by mapping all of the fragments onto a non-redundant set of transcripts. Following optimal assembly, the longest sequences were selected as unigenes. The unigenes were functionally annotated according to homology searches at the NCBI non-redundant protein (Nr), non-redundant nucleotide (Nt), and Swiss-Prot databases, with an E-value of 10E-5. The results based on the Nr BLAST algorithm were imported into Blast2GO [33] to obtain the relevant gene ontology (GO) terms, with E-values < 10E-6. The unigenes were also used as queries to screen the Clusters of Orthologous Groups (COG) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases to predict and classify functions and pathway assignments [34].

Quantification of Gene Expression Levels and Analyses of Differentially Expressed Genes
Gene expression levels were estimated with RSEM software [35]. Clean sequences were mapped back to the assembled data. Read counts for each gene were obtained from the mapping results and normalized against the FPKM values. The read counts were adjusted with the scaling normalization method of the edge R package before examining the differentially expressed genes (DEGs) in each sequenced library. The DEGs were analyzed with the DESeq2 R package. P-values were adjusted based on the q-value. The significant DEGs were identified based on a q-value < 0.05 and |log 2 fold-change| > 1. The sequences from the analysis were deposited in the NCBI database (accession: SRX1755981). To investigate the patterns of DEGs at various tree ages, K-means clustering was used to cluster the identified DEGs.

Quantitative Real-Time PCR
A PrimeScript™ RT reagent Kit (TaKaRa, Dalian, China) was used for reverse transcriptase reactions. A quantitative real-time PCR (RT-qPCR) assay was completed with SYBR ® Premix Ex Taq™ II (Tli RNase H Plus; TaKaRa). The P. orientalis αTUB gene was used as an internal reference. Primers were designed with the Primer 3 program (http://frodo.wi.mit.edu/primer3/). Relative gene expression levels (for at least three biological replicates) were calculated according to the 2 −∆∆Ct method [2] and are expressed as the mean ± standard deviation. Data were analyzed by ANOVA in the SPSS software with p < 0.05 set as the threshold for identifying significant differences.

Biochemical Changes during Senescence
Protein and ROS levels are often used as markers for the progression of senescence [3]. The soluble protein level was high in 20-year-old samples, and then gradually decreased to a minimum in 700-year-old samples ( Figure 2). The MDA content in 1100-year-old samples was 1.66-fold (p < 0.05) greater than that in 20-year-old samples. In contrast, SOD and POD activities gradually decreased with age. Thus, the extremely ancient P. orientalis trees may not be able to efficiently eliminate ROS, suggesting that they may be undergoing senescence.

De novo Assembly of the P. orientalis Transcriptome
Illumina HiSeq™ 2000 sequencing generates comprehensive and reliable data for plant gene expression analyses [36,37]. To characterize the aging events through DEGs, samples with three biological replicates for each age were sequenced. The resulting sequencing data comprised more than 6.22 Gb with an average GC content of 45.32% (Table S1a). Additionally, 62,184,973 clean reads were obtained (97.2%; Figure S1a). After trimming and assembling, 79,413 transcripts and 48,044 unigenes were generated ( Figure S1b; Table S1b). Their length distributions are shown in Figure S1c,d. The N50 value of the unigenes was 1714 bp with an average length of 896 bp (ranging from 201 to 21,855 bp; Table S1b), and 33.96% of the unigenes were longer than 300 bp (Figure S1b). The relative expression level of each unigene, which was estimated based on the FPKM approach, ranged from 0 to 14,810.83 FPKM with an average of 8.24 FPKM. Of the 48,044 unigenes, 38,202 (79.51%) were expressed at low levels (i.e., less than 10 FPKM) (Table S2).

Functional Classification of P. orientalis Unigenes
To verify and annotate the unigenes, all of the assembled sequences were initially used as queries in a BLASTx search of the Nr and Swiss-Prot protein databases. Among the 48,044 unigenes, 24,080 (50.12%) had significant hits in the Nr database, and 18,483 (38.47%) had significant matches to proteins in the Swiss-Prot database (Table 1). The annotated unigenes were grouped into the three main GO categories. A total of 20,029 unigenes were annotated according to the GO database ( Table 1). The enriched GO terms in the biological process category were cellular processes (11,580 unigenes), metabolic processes (11,322 unigenes), and single-organism processes (5710 unigenes). Only a few unigenes were involved in growth (52 unigenes), cell killing (22 unigenes), and rhythmic processes (8 unigenes). Regarding the cellular component category, the enriched GO terms were cell (6734 unigenes), cell part (6717 unigenes), and organelle (4910 unigenes). Only 1, 1, and 6 unigenes were annotated with the synapse, synapse part, and nucleoid GO terms, respectively. In terms of the molecular function category, the enriched GO terms were binding, catalytic, and transporter activities, with 11,929, 10,181, and 1333 unigenes, respectively. Only 4 and 8 unigenes had metallochaperone and receptor regulatory activities, respectively ( Figure 3).
All unigenes were classified based on the COG database, ultimately resulting in the classification of the putative proteins into at least 25 functional groups ( Figure 4). Specifically, the unigenes were categorized into the general function prediction (R, 1902 unigenes), post-translational modification (O, 1245 unigenes), signal transduction mechanisms (T, 770 unigenes), and translation (J, 651 unigenes) categories.  To further explore the involvement of the identified unigenes in metabolic pathways, we screened the KEGG database with the unigenes as queries. Based on functional annotations, the KEGG pathways were categorized into five main classes, among which the metabolism pathway included the most unigenes. A total of 835 and 610 unigenes were associated with carbohydrate metabolism and amino acid metabolism, respectively. Another 498 and 482 unigenes were related to energy metabolism and lipid metabolism, respectively ( Figure 5).

Identification of DEGs during Senescence
To investigate the molecular differences among the four analyzed tree ages, the gene expression levels were compared between ages to identify significant DEGs, with p < 0.05 as the threshold for significance (Figure 6a Table S3). Moreover, the comparison between 1100-and 20-year-old trees had the most DEGs that were specific to a particular comparison (165 of 214). Thus, aging may affect gene expression at the transcriptome level. Furthermore, the KEGG pathway enrichment analysis indicated that the photosynthesis (ko00190) and oxidative phosphorylation (ko00195) pathways were significantly enriched in the three comparisons (200-/20-year-old, 700-/20-year-old, and 1100-/20-year-old), but they were most highly enriched in the 1100-/20-year-old comparison (0.286 for the photosynthesis and 0.060 for the oxidative phosphorylation pathways). Additionally, metabolic pathways were specifically enriched in the 1100-/20-year-old comparison (Figure 6d).

Functional Classification of DEGs
To further investigate the functional classifications of the DEGs, we performed a K-means cluster analysis, which revealed eight clusters with distinct expression patterns (Figure 7a, Table S4). Cluster 4 was the largest subset (136 genes), and the expression levels of the genes in this cluster were considerably down-regulated in the 1100-year-old samples. The genes in Cluster 3 (32 genes) had expression patterns consistent with those of the Cluster 4 genes. The functional classification analysis indicated that most of the genes in Clusters 3 and 4 were related to redox reactions, photosynthesis, and stress responses (Figure 7b). Additionally, the expression levels of the Cluster 1 genes, 26 genes mainly associated with photosynthesis, gradually decreased along with tree age. The expression levels of the 79 Cluster 5 genes were down-regulated in the 200-and 700-year-old samples. These genes were mainly transport-related genes. The expression levels of the genes in Cluster 2 (55 genes) gradually increased in 20-to 1100-year-old samples. These genes were primarily associated with transcription and senescence. Increases in the expression of genes contributing to signal transduction and protein processing occurred mainly in Clusters 2, 6, 7, and 8. Thus, the DEGs are involved in various functions.

Expression Profiling of Differentially Regulated Genes by RT-qPCR
To validate the sequencing data and investigate the dynamic gene expression profiles, we analyzed the expression patterns of the following eight genes in a RT-qPCR assay: ATPase (ATP synthase subunit alpha), SNC1 (protein suppressor of npr1-1), MT2 (metallothionein-like protein class II), FBX (probable F-box protein), COX2 (cytochrome c oxidase subunit 2), ND2 (NADH dehydrogenase subunit 2), NDUFB6 (NADH-ubiquinone oxidoreductase chain 6), and TIR-NBS-LRR (putative truncated TIR-NBS-LRR protein-encoding gene) ( Table S5). The expression levels of the selected genes were consistent between the RT-qPCR and RNA-seq data, confirming the reliability of our DEG results ( Figure 8).

Down-Regulated Gene Clusters Resulted in Changes to Photosynthesis and the Redox Reaction
The expression levels of many photosynthesis-related genes were significantly down-regulated in Cluster 1 (Figure 9), as were the expression levels of many other genes encoding subunits of photosystem (PS) II complexes. In particular, the expression levels of the genes encoding PSII reaction center protein Z (PsbZ) and the PSII CP43 chlorophyll apoprotein decreased in a linear fashion along with tree age. In contrast, the expression levels of PSI-related genes, including PsaA and psaK were no substantial differences between 200-to 700-year-old samples in Cluster 3 and 4, illustrating that PSII is affected more than PSI during P. orientalis senescence. Moreover, the expression levels of genes involved in the redox reaction, such as COX1 (cytochrome c oxidase subunit 1) and COX2, as well as the gene encoding ribosomal protein S12, were significantly down-regulated in Clusters 1 and 3 (i.e., log 2 fold-change > 1 in ancient P. orientalis samples relative to the expression level in the 20-year-old samples). Furthermore, Cluster 4 was composed of 22 redox reaction-related genes with expression levels that significantly decreased from 20-to 1100-year-old samples (i.e., log 2 fold-change > 1). Genes involved in photosynthesis and the redox reaction, such as ATP94A1 (cytochrome P450 94A1) and DFR (dihydroflavonol-4-reductase), were expressed at lower levels in the 700-year-old samples than in the 20-year-old samples.

Up-Regulated Gene Clusters Resulted in Changes to Transcription, Signal Transduction, and Senescence
The expression levels of transcription-related genes involved with ethylene-responsive transcription factor 1A (EREBF) and DNAJ increased with the age of the P. orientalis trees in Clusters 2, 7, and 8, implying that these genes may encode proteins that function prior to senescence (Figure 8). Signal transduction-related genes involved with TMV (tobacco mosaic virus) resistance protein, disease resistance protein-and calmodulin-like protein (CML)-encoding genes in Clusters 2 and 7 were expressed at higher levels in young P. orientalis trees than in ancient trees (i.e., log 2 fold-change > 1). The expression levels of protein processing-related genes (e.g., genes encoding the RING (really interesting new gene)-H2 finger protein and aspartic proteinase nepenthesin-2) and a SIRK gene increased with age in P. orientalis (Clusters 2 and 7). The identified genes that were more highly expressed in the ancient trees than in the 20-year-old trees are consistent with the likelihood that the ancient trees are undergoing senescence.

Validation of the Expression Kinetics of the Photosynthesis-Related Genes during Senescence
To obtain reference data for photosynthesis-related genes that showed age-dependent transcript level profiles in P. orientalis, we analyzed the RbcS and YCF2 expression profiles. The RbcS and YCF2 expression levels decreased in P. orientalis with age ( Figure 10a; Table S6). These results suggest that photosynthesis is regulated by nuclear-and chloroplast-related genes.

Age-Related Chlorophyll Level Changes in P. orientalis
To clarify the expression-level changes of photosynthesis-related genes associated with age in P. orientalis, we analyzed the fluctuations of the chlorophyll and carotenoid levels. Chlorophyll a and b levels were highest in the 20-year-old samples (p < 0.05; Figure 10b) and decreased steadily with age. The carotenoids were highly accumulated in 20-year-old samples and had the lowest concentration in the 200-year-old samples. The trends in the changes to the chlorophyll levels during senescence implies that our DEG results are reliable.

Discussion
Senescence of trees can be affected by many external cues and internal factors, including heat stress, cold stress, and age [9]. In this study, the transcriptome results provide comprehensive and reasonable data resources for gene expression analyses of different-aged P. orientalis. Noticeably, the expression of genes in photosynthesis and oxidative phosphorylation pathways was down-regulated substantially along with age of P. orientalis. Similarly, the homologous of these genes were also down-regulated during senescence in Populus tremula and Arabidopsis thaliana, which is consistent with their involvement in the senescence of ancient trees [10,26]. This shows that the physiology of ancient P. orientalis is greatly influenced by age.
Chlorophyll level is an effective marker for the progression of senescence because it decreases during senescence in leaves of ancient P. orientalis, as well as in P. tremula and A. thaliana [10,26]. In parallel with chloroplast degeneration, the transcript levels of photosynthesis-related genes, including PsaA, PsbZ, RBcS, and ATPF1A, declined along with tree age. The similar results were reported for the senescent leaves of P. tremula and P. orientalis [4,38]. The PsaA gene encodes a PSI reaction center protein, and PsaA/PsaB polymers are necessary for the assembly of the PSI complex. In apple leaves, PsbZ is the main site of the photosynthetic apparatus that is damaged because of senescence [39]. Additionally, the light system is an important reactive site for the conversion of light energy into chemical energy [40]. Energy is required for ATP degradation and transport out of the cell during senescence. The ATPF1A gene encodes a key enzyme for energy metabolism-its down-regulation suggested significant physiological differences between ancient and young trees, similar to that observed in Pinus ponderosa [41]. Moreover, ATP is a significant factor in the regulation of mitochondrial redox activity and ROS production.
Decreased expression of photosynthesis-related genes will improve the reducibility of the electron transfer chain and increase the chances of electron leakage and ROS generation [42]. In this study, MDA content increased in leaves of ancient P. orientalis, whereas chlorophyll levels and SOD and POD activity declined (Figure 11), which is consistent with an earlier study of ancient pine [4]. These findings indicate that the ROS scavenging capabilities of ancient trees are weakened, resulting in ROS accumulation in ancient trees during aging [3]. As shown in Figure 10b, the leaves of 20-year-old trees have the highest carotenoid content, followed by 700-, 1100-and 200-year-old trees. Carotenoids comprise a large group of compounds, which may differ in the capacity to interact with free medicals and singlet oxygen and therefore act as effective antioxidants [43]. However, further studies are needed to confirm its role in P. orientalis senescence. In addition, Cys and aspartic proteases may play an important role during chloroplast degradation [26]. Our data indicate that most of the cellular proteins, such as aspartic proteinases, are highly expressed in ancient P. orientalis. The SIRK genes participate in apoptosis and programmed cell death [44]. Senescence may be followed by an increase in protein hydrolysis. The study of senescence-associated gene functional networks homologous to that of P. orientalis will increase our understanding of anti-stress activities in P. orientalis.
Senescence decreases the photosynthetic ability of P. orientalis to change redox reactions for active growth and development [45]. During P. orientalis aging, the genes involved in redox reactions are down-regulated. Most of these genes are related to oxidative phosphorylation, which is consistent with their roles during senescence in soybean leaves [46]. The expression levels of the genes encoding three main components of the COX cascade, COX1, COX2, and COX3, are down-regulated in response to low-temperature stress [47,48]. COX family plays a major role in the process of respiration, and our results suggest a reduction in the capability to convert chemical energy in ancient P. orientalis. Additionally, cytochrome P450 genes, including CYP716B1, CYP86B1, and CYP94A1, are involved with various oxidation reactions in plants [49]. They are also important for plant resistance to diseases and other stresses. In Vicia sativa, CYP94A1 might play a role in cutin biosynthesis, with implications for stress resistance [50]. These genes were highly expressed in 20-year-old trees, suggesting the young P. orientalis trees have a stronger resistance to disease than ancient trees. These observations imply that the inhibited redox reaction and abnormal energy metabolism results in various physiological metabolic changes during aging of slowly growing ancient trees.

Conclusions
In this study, the results of the physiological tests showed that the chlorophyll levels as well as the SOD and POD activities decline in leaves of ancient P. orientalis. De novo transcriptome analysis revealed 418 senescence-induced genes in P. orientalis trees. Age-specific differences in the transcriptome were most pronounced between 20-and 1100-year-old trees of P. orientalis. The expression levels of photosynthesis-, redox-, and transporter-related genes were down-regulated in ancient P. orientalis trees. The up-regulated genes are related to protein hydrolysis, senescence, transcription, and signal transduction, which imply that the ancient P. orientalis trees were undergoing senescence. In addition, the physiological data show that the chlorophyll a and b levels decrease steadily with P. orientalis age. The data presented herein suggest that the young P. orientalis trees are more resistant to long-term environmental stresses than the ancient trees. Through integrating RNA-Seq data with physiological and biochemical data, we provide insight into the molecular mechanisms regulating senescence and stress resistance in ancient P. orientalis trees. Consequently, the processes associated with senescence will need to be more comprehensively characterized in future studies.