Flavonoid Biosynthesis Is Likely More Susceptible to Elevation and Tree Age Than Other Branch Pathways Involved in Phenylpropanoid Biosynthesis in Ginkgo Leaves

Ginkgo leaves are always resources for flavonoids pharmaceutical industry. However, the effect of the elevation and tree age changes on flavonoid biosynthesis have not been detailly explored in Ginkgo leaves. In addition, whether these environmental pressures have similar effects on the biosynthesis of other non-flavonoids polyphenolics in phenylpropanoid biosynthesis is not known at present. In this research, de novo transcriptome sequencing of Ginkgo leaves was performed coupled with ultra-performance liquid chromatography/quadrupole time-of-flight mass spectrometry analyses to obtain a comprehensive understanding of the influence of elevation and tree age on phenylpropanoid biosynthesis. A total of 557,659,530 clean reads were assembled into 188,155 unigenes, of which 135,102 (71.80%) were successfully annotated in seven public databases. The putative DFRs, LARs, and ANRs were significantly up-regulated with the increase of elevation in young Ginkgo tree leaves. The relative concentration of flavonoid derivatives with high parent ion intensity was likely to imply that the elevation increase promoted the biosynthesis of flavonoids. Complex gene variations involved in flavonoid biosynthesis were observed with the tree age increase. However, flavonoid derivatives analysis predicted that the rise of tree age was more likely to be detrimental to the flavonoids manufacture. Otherwise, multiple genes implicated in the synthesis of hydroxycinnamates, lignin, and lignan exhibited fluctuations with the elevation increase. Significantly up-regulated CADs and down-regulated PRDs potentially led to the accumulation of p-Coumaryl alcohol, one of the lignin monomers, and might inhibit further lignification. Overall, the putative DFRs seemed to show more considerable variability toward these stress, and appeared to be the main regulatory point in the flavonoid biosynthesis. Light enhancement caused by elevation increase may be the main reason for flavonoids accumulation. Flavonoid biosynthesis exhibited a greater degree of perturbation than that of hydroxycinnamates, lignins and lignans, potentially suggesting that flavonoid biosynthesis might be more susceptible than other branch pathways involved in phenylpropanoid biosynthesis. This research effectively expanded the functional genomic library and provide new insights into phenylpropanoid biosynthesis in Ginkgo.

Ginkgo leaves are always resources for flavonoids pharmaceutical industry. However, the effect of the elevation and tree age changes on flavonoid biosynthesis have not been detailly explored in Ginkgo leaves. In addition, whether these environmental pressures have similar effects on the biosynthesis of other non-flavonoids polyphenolics in phenylpropanoid biosynthesis is not known at present. In this research, de novo transcriptome sequencing of Ginkgo leaves was performed coupled with ultraperformance liquid chromatography/quadrupole time-of-flight mass spectrometry analyses to obtain a comprehensive understanding of the influence of elevation and tree age on phenylpropanoid biosynthesis. A total of 557,659,530 clean reads were assembled into 188,155 unigenes, of which 135,102 (71.80%) were successfully annotated in seven public databases. The putative DFRs, LARs, and ANRs were significantly up-regulated with the increase of elevation in young Ginkgo tree leaves. The relative concentration of flavonoid derivatives with high parent ion intensity was likely to imply that the elevation increase promoted the biosynthesis of flavonoids. Complex gene variations involved in flavonoid biosynthesis were observed with the tree age increase. However, flavonoid derivatives analysis predicted that the rise of tree age was more likely to be detrimental to the flavonoids manufacture. Otherwise, multiple genes implicated in the synthesis of hydroxycinnamates, lignin, and lignan exhibited fluctuations with the elevation increase. Significantly up-regulated CADs and down-regulated PRDs potentially led to the accumulation of p-Coumaryl alcohol, one of the lignin monomers, and might inhibit further lignification. Overall, the putative DFRs seemed to show more considerable variability toward these stress, and appeared to be the main regulatory point in the flavonoid biosynthesis. Light enhancement caused

INTRODUCTION
Phenylpropanoids biosynthesis pathway functions as a natural factory to produce a variety of secondary metabolites in response to biological and abiotic stimuli (Deng and Lu, 2017). Derived from phenylalanine in most plants or tyrosine in partial monocots, the central phenylpropanoids mainly include flavonoids, monolignols, hydroxycinnamates (HCAs), lignins, and lignans, which acting as components of cell walls, protectants against UV radiation, signaling molecules phytoalexins against herbivores and pathogens (Vogt, 2010;Deng and Lu, 2017).
Generally, the first enzymatic step of phenylpropanoid biosynthesis begins with the deamination of phenylalanine by phenylalanine ammonia lyase (PAL) to yield cinnamic acid (Barros et al., 2016). In the second step, cinnamate 4-hydroxylase (C4H) hydroxylates and transforms the cinnamic acid into p-coumaric acid. Then, 4-coumaroyl CoA ligase (4CL) catalyzes p-coumaric acid into p-coumaroyl-CoA, which is a crucial branch point leading to the generation of flavonoids, lignins and lignans (Vogt, 2010).
Plants have evolved a variety of antioxidant systems to cope with reactive oxygen species (ROS) generated under ambient pressures (Brunetti et al., 2015). As an important antioxidant system, phenylpropanoid biosynthesis is regulated by many factors, including light, irrigation, temperature, and fertilization (Kim et al., 2015). Of which, light appears to be the most important factor that affects the dynamic synthesis of multiple polyphenolic substances. For instance, the light of different wavelengths has been widely used to accumulate phenolic compounds in lettuce (Johkan et al., 2010), chrysanthemum (Jeong et al., 2012), and buckwheat (Thwe et al., 2014). Besides, it has been widely reported that the enhancement of light can significantly increase the flavonoids content in fruits (Awad et al., 2001;Azuma et al., 2012;Zoratti et al., 2014). Studies of circadian rhythms show that light signals promote the synthesis of lignin in the night in Arabidopsis (Rogers et al., 2005). Lignans were observed a notable increase under continuous white fluorescent and monochromatic red or blue light emitting diode (LED) light (Hata et al., 2013). However, the light-induced fluctuation of HCA biosynthesis shows uneven regularity (Kolb et al., 2001;Hemm et al., 2004).
Ginkgo biloba (Ginkgo) is widely used as traditional Chinese medicine for centuries, especially on asthma and cardiovascular disease (Kleijnen and Knipschild, 1992;Birks and Grimley Evans, 2007). The flavonoids and terpene trilactones are currently considered the most prominent pharmacological components of Ginkgo (van Beek, 2002). The studies on the phenylpropanoid biosynthesis of Ginkgo mainly focused on flavonoids, especially the influences of populations (Rimkiene et al., 2017), seasonality (Sati et al., 2013), and growth stages (Mikašauskaitė et al., 2013).
In general, the elevation increase always contributes to light enhancement and a decrease in the average temperature. Little evidence indicated the full impact of variable elevation and tree age on the flavonoid biosynthesis of Ginkgo and its possible reasons. The non-flavonoids polyphenolics need further discussion as well. Moreover, few discussions have focused on the effect of tree age on the phenylpropanoids biosynthesis in woody plant leaves. And whether different branches of phenylpropanoid synthesis will show consistent responses to the same environmental stress need more discussion. Therefore, in this study, wild Ginkgo trees were sampled to investigate the specific profile changes of phenylpropanoids biosynthesis induced by the variation of elevation and tree age via non-referenced transcriptome analysis coupled with ultraperformance liquid chromatography/quadrupole time-of-flight mass spectrometry.

Ginkgo Leaves
The Ginkgo trees grow wildly in Pot bottles of mountain Nature Reserve in Hunan Province, China. A total of 12 trees (26 • 55 13 N to 30 • 6 44 N, 110 • 36 16 E to 110 • 49 2 E) were chosen to collect mature leaf samples (Supplementary Table S1). According to the formula (age of tree = diameter * growth factor) proposed by the International Society of Arboriculture in the United States for estimating the age of trees, it can be inferred that the diameter at breast height (DBH) positively correlates with the age of the trees in the same (or similar) environment. So, the DBH (converted from trunk circumference) was used as a measure of the relative extent of tree ages. The elevation and the DBH were determined by Global Positioning System (GPS) instrument and tape measure respectively. All samples were divided into four groups -low elevation and young age (LY), high elevation and young age (HY), low elevation and older age (LO), high elevation and older age (HO). Each group involved three tree individuals. Each sample was mixed by thirty leaves of similarly size and health, which were cut off from the sunny side 5 m above the ground of the same tree (Figure 1). After harvesting and shorttime surface cleaning by 75% ethanol and sterile water, these leaves were frozen in liquid nitrogen immediately and stored at FIGURE 1 | A schematic diagram of sample collection and data measurement. Small black circles indicated the sampling sites.
−80 • C until used. The same copies of each leaf sample were separately used for RNA extracting and UPLC-MS.

RNA Extraction, Library Preparation, and Sequencing
Total RNA extraction was carried out by E.Z.N.A. R Plant RNA Kit (Omega Bio-Tek, Inc.) according to the manufacturer's instruction coupled with genomic DNA digesting with DNase I. Then 1% agarose gel electrophoresis was used to monitor RNA degradation and contamination. RNA purity, concentration and integrity were measured by NanoPhotometer R spectrophotometer (Implen, Westlake Village, CA, United States), Qubit 2.0 Fluorometer (Life Technologies, Foster City, CA, United States) and Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, United States). Dynabeads Oligo (dT) (Invitrogen, United States) was used to isolate mRNA from total RNA. Then the NEB Next R Ultra RNA Library Prep Kit for Illumina R (NEB, United States) was used to construct 12 cDNA libraries following manufacturer's recommendations. After some necessary processes (index-codes adding, mRNA purification, fragmentation, cDNA synthesis, exonuclease/polymerase cleavage, adenylation of 3 ends of cDNA fragments, adaptor ligation, size selection, PCR and purification) and the terminal quality assessment on the Agilent Bioanalyzer 2100 system, the library preparations were sequenced on an Illumina Hiseq PE150 platform and paired-end reads were generated.

Quality Control, de novo Assembly, and Annotation
Raw data in fastq format (Cock et al., 2010) were processed through in-house Perl scripts. Clean data were obtained by removing reads containing adapter, reads containing ploy-N and low-quality reads from raw data. At the same time, Q20, Q30, GC-content and sequence duplication level of the clean data were calculated. Transcriptome assembly was accomplished by Trinity (v2.4.0) (Grabherr et al., 2011;Haas et al., 2013) with min_kmer_cov set to 1 and all other parameters set default. All transcripts were clustered into unigenes in Corset (v1.05) with all parameters default as recommended (Davidson and Oshlack, 2014). Then, these unigenes were searched against seven databases for function annotation. NCBI blast 2.2.28+ (Altschul et al., 1997) was used for the alignments of unigenes to Nt database with an E-value threshold of 1E-5. The program diamond (v0.8.22) (Buchfink et al., 2015) was selected to perform the comparison against Nr (E-value 1E-5), KOG/COG (E-value 1E-3), and Swiss-Prot (E-value 1E-5) databases. The hmmscan in HMMER 3.0 was operated to search Pfam (Finn et al., 2008) with an E-value threshold of 0.01. The GO (Gene Ontology) annotations were carried out in Blast2GO (v2.5) (Gotz et al., 2008) with an E-value threshold of 1E-6 based on the Nr and Pfam annotations. KAAS 1 (Moriya et al., 2007) was used for the KEGG annotations (E-value 1E-6). Transcription factors (TFs) was predicted in iTAK 1.2 with all parameters set default (Perez-Rodriguez et al., 2010;Jin et al., 2014;Zheng et al., 2016).

Gene Expression Levels and Differentially Expressed Genes Identification
The assembled transcriptome spliced by Trinity was set as the reference sequence (ref) (Haas et al., 2013). The clean reads of each sample were mapped back to this ref in RSEM (v1.2.15) (Li and Dewey, 2011) with the Bowtie2 mismatch set 0 as default. The mapping results including read counts of each sample were normalized by calculating the FPKM (expected number of Fragments Per Kilobase of transcript sequence per Millions of base pairs sequenced) to obtain relative expression levels of unigenes (Trapnell et al., 2010) (Supplementary Datasheet S1). An FPKM > 0.3 was defined as the threshold of significant gene expression (Ramskold et al., 2009). Differentially expressed genes (DEGs) identification between different groups was performed by DESeq2 R package (v1.6.3) (Love et al., 2014) after FPKM tables were normalized as recommended (Dillies et al., 2013). The unigenes with an adjusted p-value < 0.05 were assigned as DEGs. To figure out the major functional groups of DEGs, GOseq R package was used to implement Gene Ontology (GO) enrichment analysis of the DEGs based on the Wallenius noncentral hyper-geometric distribution (Young et al., 2010). KOBAS (Mao et al., 2005) was used to test the statistical enrichment of DEGs in KEGG pathways.

Quantitative Real-Time PCR Validation
Total RNA was reverse transcribed into cDNA after residual genomic DNA degradation by ReverTra Ace R qPCR RTMaster Mix with gDNA remover Kit (TOYOBO, Japan). Then quantitative PCR experiments were performed on Bio-Rad iQ5 Optical System (BIO-RAD Laboratories, Inc., United States) with KOD SYBR R qPCR Mix Kit (TOYOBO, Japan). Primers used in this process were designed in Primer Premier 6 (United States), and the primer sequences are given in Supplementary Table S2. The relative expression levels were calculated according to the absolute quantification of mRNA based on standard curve line as commonly recommended (Bustin, 2000). Glyceraldehyde 3phosphate dehydrogenase (GAPDH) gene was used as an internal standard (Meng et al., 2016;Meng and Fricke, 2017). As designed, each group included three biological replicates and each cDNA sample was carried out with three independent technical repetitions.

LC-MS Samples Preparation
To minimize the error, andrographolide was dissolved in acetonitrile as an internal standard solution with the terminal concentration of 133 µM. Of each leaf sample, 200 mg liquid nitrogen-grinded powder was homogenized in 30 ml 70% ethanol solution (v: v = 70: 30) followed by 1 min vortex and 1 h ultrasound extraction as previously described with minor modifications (Yu et al., 2003;Tohge et al., 2005;Zhou et al., 2014). After a 13000 rpm centrifugation for 10 min at 4 • C, 1 ml supernatant of the solution was transferred and evaporated to dryness under nitrogen gas at 37 • C . The residue was re-dissolved in 1 ml acetonitrile and centrifuged at 13000 rpm for 10 min at 4 • C. 950 µL supernatant was transferred to mix with 50 µL internal standard solution to be the final sample for further UPLC-QTOF/MS analysis.
AB SCIEX TripleTOF TM 5600+ system (AB SCIEX Technologies, United States), equipped with an electrospray ionization (ESI) source, was coupled to the UPLC system and used to scan parent ion molecular weight from 100 to 1500. Other MS parameters were set as below: electrospray ionization temperature (

Quality Control
To monitor the stability and quality of the non-targeted bioanalysis, the quality control (QC) samples were prepared by pooling 0.2 g powder of each sample and extracted in the same method . Linearity was assessed by calculating coefficients (R 2 ) of the diluting of QC samples in series (1,2,5,8,10,20,50,100,200, and 500 times). Repeatability was evaluated by six replicates of QC samples prepared and extracted parallelly as above. Since all tests are completed in 1 day, precision was investigated by intra-day variability, which was inspected by injecting the QC sample six times along with experimental samples detection in 1 day.

UPLC-MS Data Preprocessing and Biomarkers Identification
The raw data outputted from LC-MS was pretreated by MarkerView (version 1.2.1.1, AB SCIEX Technologies, United States), including peak recognition (retention time 2-28 min, noise threshold 100), alignment, calibration of the internal standard, filtering and normalization to total area. A three-dimensional data set contained sample information, peak retention time (RT), peak relative intensities and mass-tocharge ratio (m/z) was obtained to perform a series of statistical analysis. PeakView (version 1.2.0.3, AB SCIEX Technologies, United States) was recommended to visualize raw data of target components in two-stage mass-to-charge ratio map. Then, based on fragment ion information, such components were identified by comparing to HMDB 2 , PubChem 3 , NIST 4 , MassBank 5 , and METLIN 6 databases (Linstrom and Mallard, 2001;Horai et al., 2010;Domingo-Almenara et al., 2018;Kim et al., 2018;Wishart et al., 2018).

RESULTS
Sequencing and de novo Assembly 83.64 Gb high-quality sequences were obtained from sequencing and pretreatment (Supplementary Table S3) with 0.01-0.02% error rates. The average Q20 and Q30 were 96.94 and 92.37%. The average GC content was 46.51%. All the clean reads were pooled together for de novo assembly in Trinity (Grabherr et al., 2011). A total of 188,155 unigenes with a mean length of 1,278 bp and an N50 of 1,855 bp were generated, among which 81,241 unigenes (43.18%) were longer than 1,000 bp. The detailed length distribution of the transcripts and unigenes is shown in Figure 2.

Functional Annotation and Classification
There were 135,102 unigenes, accounting for 71.8%, annotated in at least one database, and 17,500 (9.3%) unigenes were observed matches in all seven databases (Supplementary Table S4 and Post-translational modification, protein turnover, chaperones, followed by (J) Translation, ribosomal structure and biogenesis, and (R) General function prediction only. The KEGG pathway database was always used to predict the interactions of the genes and the products. Overall, 51,137 unigenes were functionally classified into 132 pathways corresponding to multiple secondary hierarchies (Supplementary Figure S1c).

Gene Expression Statistics and DEGs Analysis
All clean reads of each sample were analyzed back to the transcriptome assembly. As listed in Supplementary  Table S5, over 70% of the sequences per sample were matched, which indicated reliable data processing and quality. Density distributions of gene expressions in each group were shown in Supplementary Figure S2, where the difference between LO and HO was the most significant. Significant transcriptional differences between all groups were reflected in the PCA plot (Supplementary Figure S3).
HY vs. LY and HO vs. LO were used to predict the effect of elevation on flavonoid biosynthesis in Ginkgo leaves at two plant age levels. Parallelly, LO vs. LY and HO vs. HY were selected to evaluate the effect of plant age on flavonoid biosynthesis in Ginkgo leaves at two elevation levels. Figure S4 revealed that 1,264 DEGs were observed when the elevation raised in HY vs. LY, including 547 up-regulated and 717 down-regulated unigenes (Supplementary Table S6). And in HO vs. LO, 2,046 DEGs were observed, 186 up-regulated and 1,860 down-regulated. Venn diagrams (Supplementary Figure S5) showed that 46 DEGs were identified in both HY vs. LY and HO vs. LO, of which 10 DEGs were up-regulated and 34 DEGs were down-regulated in both comparisons. In terms of quantity, it seemed that the high age of Ginkgo trees significantly affected the response of genes expression to elevation change, mainly reflected in the significant increase of down-regulated unigenes.

Volcano plots in Supplementary
As to the influence of plant age, 6,438 DEGs were filtered out from LO vs. LY, 4,915 of which were up-regulated and 1,523 of which were down-regulated. However, only a total of 139 DEGs were found in HO vs. HY, including 79 DEGs up-regulated and 60 DEGs down-regulated. Thereby, there were only nine DEGs observed both in LO vs. LY and HO vs. HY with two up-regulated and two down-regulated synchronously. Obviously, the higher elevation appeared to make gene expression insensitive to the change of Ginkgo trees age.

GO Enrichment of DEGs
Thirty most enriched GO categories of each comparison were plotted in Supplementary Figure S6. The results showed that the DEGs in HY vs. LY were significantly enriched in "obsolete peroxidase reaction" (GO:0006804) of biological process, "peroxidase activity" (GO:0004601) and "oxidoreductase activity, acting on peroxide as acceptor" (GO:0016684) of molecular function. The independent enrichment of upregulated unigenes indicated a significant enhancement of hormone regulation and metabolism due to the rise of elevation in HY vs. LY, especially in "regulation of hormone levels" (GO:0010817), "cellular hormone metabolic process" (GO:0034754) and "hormone metabolic process" (GO:0042445), which was hypothesized to play an essential role in the regulation of flavonoid biosynthesis (Brunetti et al., 2018). However, in HO vs. LO, there were 65 enriched categories, of which "catalytic activity" (GO:0003824) was enriched the most in molecular function, and "carboxylic acid metabolic process" (GO:0019752) was enriched the most in biological process, while "tricarboxylic acid cycle enzyme complex" (GO:0045239) was the only enriched category in cellular component. The down-regulated unigenes were significantly enriched in decentralized GO terms, including multiple terms relate to the biological process of nucleic acid. It seemed that the older Ginkgo trees need more nucleic acids process to cope with the abiotic stress caused by elevation change.
In LO vs. LY, 45 categories in biological process, 10 categories in cellular component and 49 categories in molecular function were observed significant enrichment. "Catalytic activity" (GO:0003824) was the most enriched category in molecular function, followed by "pyrophosphatase activity" (GO:0016462) and "hydrolase activity, acting on acid anhydrides, in phosphorus-containing anhydrides" (GO:0016818). "Localization" (GO:0051179) was the most enriched category in biological process, followed by "transport" (GO:0006810) and "establishment of localization" (GO:0051234). "Holliday junction helicase complex" (GO:0009379) was the most enriched category in cellular component, followed by "Holliday junction resolvase complex" (GO:0048476) and "DNA helicase complex" (GO:0033202). Further, the up-regulated unigenes were enriched in 120 categories while the downregulated ones enriched in nine categories. Nevertheless, no categories were significant enriched in HO vs. HY, although "toxic substance binding" (GO:0015643) of molecular function and "immunoglobulin complex" (GO:0019814) of cellular component possessed the lowest corrected p-values. Moreover, "binding" (GO:0005488) was the most abundant category with 78 DEGs, followed by "cellular process" (GO:0009987) and "metabolic process" (GO:0008152).

KEGG Enrichment of DEGs
All DEGs, as well as the up-and down-regulated DEGs separately, were mapped to the KEGG database to identify genes involved in signal transduction or metabolic pathways. Twenty pathways with the most enrichment degree of each comparison were shown in bar plots.
In brief, Ginkgo trees of low age showed the disturbance of benzene propane, flavonoid biosynthesis, and circadian rhythm in the face of elevation rise. Flavonoid biosynthesis showed a significant up-regulation at elevated altitude. Ginkgo individuals of different ages at low altitudes show major differences in amino acid synthesis and metabolism. When they became older, significantly up-regulated unigenes are associated with chlorophyll and photosynthesis. Also, fatty acid elongation, phenylpropanoid synthesis, and plant hormone signal transduction were downgraded to varying degrees.
However, in HO vs. LO, there were only one putative DFR unigene up-regulated with a log 2 (FC) value of 6.72 and one putative FLS unigene down-regulated with a log 2 (FC) value of −2.08. For LO vs. LY, one putative F3M unigene, three putative DFR unigenes and two putative ANR unigenes were identified different-level up-regulation when the plant age became older. Meanwhile, one putative C3 H unigene, one putative F3H unigene, one putative F3M unigene, one putative DFR unigene and one putative FLS unigene were down-regulated despite several |log 2 (FC)| values were less than one. Unfortunately, no DEGs were filtered in flavonoid biosynthesis in HO vs. HY.
In summary, a higher elevation led to more unigenes upregulated than down-regulated in flavonoid biosynthesis in HY vs. LY at a young plant age level. However, when plant age became older in HO vs. LO, only one putative DFR unigene and one putative FLS unigene were both down-regulated with the increase of elevation. The situation was much more complicated when the influence of plant age was considered in LO vs. LY at a low elevation level. Two putative unigenes of ANR was observed up-regulated, and the filtered putative unigenes of C3 H, F3H, and FLS were all down-regulated with the decrease of plant age. Additionally, F3M possessed one up-regulated and one down-regulated putative unigenes. For DFR, three putative unigenes were up-regulated and one was down-regulated as plant age increased. Unexpectedly, no unigenes related to flavonoid biosynthesis were found significant variation when elevation increased in HO vs. HY. Seemingly, the flavonoid biosynthesis pathway was more sensitive and easier to be influenced when wild Ginkgo individuals were at low elevation or young age level. DFRs may be the key regulatory point for flavonoid biosynthesis disturbance in Ginkgo leaves.

RT-qPCR Validation
RT-qPCR was used to verify the expression data of RNA-Seq and the subsequent splicing, assembly, and expression calculation. To exclude the influence of alternative splicing of gene unigenes, primers were designed in the homology regions of the unigenes which were transcribed from the same gene and retained highsimilarity annotation. The log 2 (FC) results ( Figure 4C) indicated the similarity between the RT-qPCR and the transcriptome situation of the filtered DEGs which held the highest highsimilarity annotation of each gene (black boxed in Figure 4B).

Flavonoid Derivatives Seemed to Show an Upward Trend With an Elevation Increase
The overlapping typical total ion chromatograms (TICs) of six QC samples obtained from LC-MS in negative mode (Supplementary Figure S8) demonstrated the acceptable error during the whole data collection. Meanwhile, six extracted ions (Supplementary Figure S9 and Supplementary Tables S7-1, S7-2, S7-3) were chosen to assess the repeatability of sampling and processing, as well as the stability and linearity of detection system by relative standard deviations (RSDs). The RTs varied from 0.03-1.40% on repeatability while 0.02-1.69% on intraday variability and 0.04-1.89% on inter-day variability. For peak area, RSDs presented 6.13-15.68, 1.06-3.38, and 7.43-16.83% on repeatability, intra-day variability and inter-day variability, respectively. The mass accuracy possessed the RSDs from 8.77E-5 to 2.52E-4%. The linearity results showed the reliability of the experimental system with all R 2 s more than 0.997.
According to the non-targeted metabolites detection and identification by UPLC/qTOF MS method, A total of 13 ingredients involved in flavonoid biosynthesis pathway and 11 flavonoid downstream derivatives that possessed high parent ion signal intensity in MS were predicted and confirmed based on fragment information in multiple databases. The details of these ingredients, including RT, theoretical and experimental m/z, mass error and ions fragment information, were summarized in Table 1 and Supplementary Datasheet S3. The relative concentrations normalized from ion intensities were presented in the form of column charts with means and standard deviations in Figure 5. Unexpectedly, the differences in concentration of all ingredients between the four groups were not significant since Ginkgo trees grow in the wild.
Trans-cinnamic acid 2-hydroxylase (C2H), C4H, C3H, COMT, and F5H were the key enzyme in the process of converting cinnamic acid into HCAs. One putative C4H unigene and five putative COMT unigenes were observed different levels of fluctuation. However, the concentration statistics of p-Coumaric acid, Caffeic acid, Ferulic acid, and Sinapic acid indicated no significant variation between four comparisons. Seven putative bGLD unigenes were screened with notable differential gene expressions, which indicated that the rise of elevation might promote the conversion from cis-beta-D-Glucosyl-2-hydroxycinnamate to cis-2-Hydroxycinnamate while an older plant age suggesting a negative regulation on this process. Nevertheless, related metabolite detection is not significant enough to support this conclusion. Noteworthily, five up-regulated CAD unigenes and five downregulated PRD unigenes may be the main cause of a marked accumulation of p-Coumaryl alcohol when the elevation rose in HY vs. LY. They also have the potential possibility to promote the production of alcohols such as Caffeyl alcohol, Coniferyl alcohol, etc., although these metabolites not been identified. In this case, the up-regulation of CADs and the down-regulation of PRDs may lead to the inhibition of lignins (Guaiacyl lignin, Syringyl lignin, etc.) biosynthesis, which is more conducive to the synthesis of lignans. However, there was no significant evidence of PLRs and SIRDs, as well as the only concentration analysis of Secoisolariciresinol to support this speculation. When the age of Ginkgo trees increased in LO vs. LY, although some putative unigenes of CAD, PRD, and PLR showed an upward trend, the similarity of functional annotations was not high. In addition, there was no indication that the concentrations of p-Coumaryl and Secoisolariciresinol have increased significantly.
Hence, the elevation increase may lead to the accumulation of lignin monomers, like 4-Coumaryl alcohol, Coniferyl alcohol, 5-Hydroxyconiferyl alcohol and Sinapyl alcohol, but inhibit the lignification by up-regulating CADs and down-regulating PRDs. Thus, a hypothesis was presumed that the rise of elevation might potentially inhibit the synthesis of lignins, leading to an accumulation of Caffeyl alcohol, Coniferyl alcohol, etc. HCAs and lignans showed no obvious variational trend in this situation. Age variability caused related genetic perturbations, but there was insufficient evidence to prove the specific impact.

Light Induction May Be the Main Reason for Flavonoids Accumulation in the Case of Elevation Rise
Flavonoids have long been recognized as the primary antioxidants in plants which are useful in preventing UV-B radiation from penetrating into sensitive tissues and protecting against programmed cell death mediated by oxidative stress (Sadik et al., 2003;Michalak, 2006;Agati and Tattini, 2010;Agati et al., 2013). The rise of elevation directly leads to higher exposure to UV radiation and lower temperature, inducing the production of flavonoids (Jugran et al., 2016). However, the flavonols appear to accumulate with the light enhancement, while be little affected by temperature (Jaakola and Hohtola, 2010). In this research, an approximate 400 m increase of elevation in HY vs. LY could not cause a visible decrease in temperature, meanwhile implying that the excess light and the higher UV radiation may be the main reason for the accumulation of flavonoids. Mainly in the form of glycosides, the identified flavonoids derivatives (viz. Robinin, Quercitrin, Rhamnegin, etc.) presented a potentially positive correlation to the elevation increase, which seems to be inconsistent with the existing theories that high light irradiance accelerates the biosynthesis of dihydroxy B-ring-substituted flavonoids, but does not affect the biosynthesis of monohydroxy B-ring-substituted flavonoids (Ryan et al., 1998;Agati et al., 2011). Understandably, although B-ring dihydroxylated flavonoids display higher antioxidant activity toward ROS, the monohydroxylated ones exert a higher capacity of UV wavelengths absorption (Nabavi et al., 2018). In addition, the increase of elevation may cause an enhancement of other abiotic stresses, which may be the possible interpretation of an up-regulated tendency of some monohydroxy B-ring-substituted flavonoids, like Apigenin and Kaempferol. Notably, it has been confirmed that the synthesis of flavonoids is closely related to the regulation of auxin and abscisic acid (ABA) (Brunetti et al., 2018;Kumar et al., 2018). The ABA, reported to promote the accumulation of flavonols (Berli et al., 2010(Berli et al., , 2011, was observed to have an active integration with the light signaling (Bechtold et al., 2008;Wang et al., 2018), which may help to explain the enhancement of flavonol biosynthesis in leaves exposed to high light irradiance, even without ultraviolet radiation (Brunetti et al., 2018).

DFR May Be the Key Regulatory Point for Light-Induced Flavonoid Biosynthesis Disturbance in Ginkgo Leaves
Numerous evidence suggested that more than one homologous DFR unigenes existed in multiple plant species (Beld et al., 1989;Tanaka et al., 1995;Cassetta et al., 2013;He et al., 2015). Under light-induced conditions, the expression of DFRs was confirmed to change significantly, leading to the accumulation of flavonoid anthocyanins, condensed tannins and other substances (Paolocci et al., 2005;Shi et al., 2014). During flower color development, DFRs exhibited distinctive variability in circadian rhythms of Anthurium andraeanum, indicating its potentially feature as a critical regulatory point of anthocyanins biosynthesis (Collette et al., 2004). However, in the green leaves, especially of a woody plant, the potential regulatory status of DFR has not been fully discussed yet. Three DFR cDNA clones were isolated from Ginkgo biloba, of which GbDFR2 and GbDFR3 were related to anthocyanin accumulation in leaves while GbDFR1 seemed to be involved in environmental stress response (Cheng et al., 2013). In this research, the gene expression profile seemed to support a view that DFRs may be the key regulatory point for light-induced flavonoid biosynthesis disturbance in Ginkgo leaves. The expression profiles of DEGs ( Figure 4B) involved in flavonoid biosynthesis and the amino acid-based phylogenetic trees (Supplementary Figure S10) of related unigenes suggested the view that DFR disturbances were significantly stronger and more sensitive than other flavonoid synthesis-related genes when facing to the change of the elevation and plant age. The DFRs of high annotation similarity, especially Cluster-24089.39479 and Cluster-24089.54759 may potentially play a major role in the accumulation of antioxidants under light induction. Moreover, in the comparison of HY vs. LY, the significantly up-regulated LAR and ANR with the rise of various flavonoid glycosides may also be a good illustration of the effects of illumination enhancement. Although many cases indicated that CHS might produce significant responses to the UV stimulation (Chappell and Hahlbrock, 1984;Batschauer et al., 1996), no CHS unigenes were found to be notably altered in this study. Most of the researches in this filed were based on herbaceous plants, while woody plants need more comprehensive and in-depth insights to unveil the potential metabolic properties. In addition, DFRs may be more sensitive to the weaker changes in light intensity caused by small elevation variability in Ginkgo.

The Age Effects of Ginkgo Trees on the Biosynthesis of Flavonoids
Nowadays, the temporal effects on the flavonoid biosynthesis were mainly focused on seasonality, circadian rhythm and the growth stages of leaves, flowers or fruits in herbs (Zhang et al., 2003;Collette et al., 2004;Nakatsuka et al., 2005;Crecelius et al., 2017), whilst only a few discussions were carried out on the age effects of the perennial woody plant, especially some longevity species like Ginkgo biloba and Picea abies.
An accumulation of flavan-3-ols (especially proanthocyanidin and (−)-epicatechin) has been reported to be simultaneous with the age increase of Cistus clusii, a perennial evergreen bush (Hernández et al., 2011). However, in this study, no significant variability of (−)-epicatechin was observed with the age increase of the Ginkgo trees (LO vs. LY in Figure 5). Moreover, although the age increase caused disturbances in several genes' expression in flavonoid biosynthesis, such as the up-regulated DFRs and ANRs (Figure 4B), the relative concentration of most identified flavonoids or flavonols indicated a declining tendency of no statistically significance with the increase of Ginkgo age (LO vs. LY). On the other hand, although the senescence process has been reported to be correlated with the accumulation of phenolics in Pinus needles (Szymura, 2009;Kumar et al., 2018), no similar indication was observed existence in this study. Until now, no detailed discussions were taken out on the signs of Ginkgo senescence. Therefore, whether flavonoids play a vital role in the senescence of Gingko individuals requires further research.

The Transcriptomic and Metabolic Influence on Non-flavonoids Polyphenolics in the Phenylpropanoid Biosynthesis
There was no significant correlation between HCAs and elevation increase in this study. Increased altitude means increased illumination. The HCAs changes under light induction have long been questioned. Excessive sunlight exposure led to a decrease of HCAs:flavonoid ratio in Ligustrum vulgare leaves (Tattini et al., 2004), which held a similar ratio trend to this study. Inconsistently, a decrease of HCAs content was observed with the light enhancement in the roots of Arabidopsis (Hemm et al., 2004). However, high visible radiation could effectively stimulate the synthesis of hydroxycinnamic acid in grape leaves (Kolb et al., 2001). Different results may result from different species and tissues for sampling. The limited increase in average light intensity caused by elevation increase may be the potential reason why no significant changes have been observed in HCAs content.
Lignin biosynthesis, as a highly energy-consuming and irreversible process, responds to many environmental factors like light, sugar content, biological clocks (Zhao and Dixon, 2011). It has been reported that high-light stress could enhance lignin synthesis in Arabidopsis seedlings (Kimura et al., 2003) and hexose carbohydrates, the photosynthesis products, could promote lignin synthesis at night, relying on circadian rhythms (Rogers et al., 2005). However, in Ginkgo leaves, no significant signs were observed to reflect the light influence on lignin synthesis. The light-induced enhancement of lignin biosynthesis may be not fully applicable to the leaves of mature woody plants. Lignin biosynthesis genes are preferentially expressed in actively lignifying tissues, while many lignin repressors are preferentially expressed in non-lignifying tissues (Zhao and Dixon, 2011). For mature Ginkgo individuals, the lignification degree of the trunk and branches always are much higher than that of the leaves, so the expression of the lignification gene should be mainly concentrated in the branches rather than the leaves. This may also be the potential reason why the age variability of Ginkgo trees had not caused significant changes in the lignin synthesis pathway in leaves in this study.
The illumination change caused by elevation rise showed insignificant disturbances for lignans biosynthesis in this study. CPI-FK cells were reported to show a positive response to white fluorescent, red LED and blue LED light in the synthesis of lignans, indicating that light enhancement has the potential to improve lignans production (Morimoto et al., 2011;Satake et al., 2015). However, this phenomenon has only been verified in Forsythia koreana and Sesamum indicum, and verifications of more species need to proceed further.
In summary, the weak illumination enhancement caused by the elevation change may be the main reason for the significant accumulation of lignin monomers. However, there were no significant changes in HCAs, lignins and lignans. The potential reason was the mild variability of environmental factors, unlike the artificial designs of dramatic environmental difference. In addition, different physiological status caused by age changes did not cause significant fluctuations in the biosynthesis of nonflavonoids polyphenolics. As few insights were thrown to study the influence of individual age on its tissue metabolism, more focus and discussions should be paid to this field.

Flavonoid Biosynthesis May Be More Sensitive to Elevation and Tree Age Than Other Branch Pathways Involved in Phenylpropanoid Biosynthesis
An elevation increase of 300-400 m may result in a moderate enhancement in illumination intensity. In this situation, DEGs analysis seemed to show a stronger disturbance in flavonoid biosynthesis than in other non-flavonoids polyphenolics pathways, which indicated that flavonoid biosynthesis might be more sensitive to elevation and tree age than other branch pathways involved in phenylpropanoid biosynthesis. Metabolites profile seemed to support this conjecture. Otherwise, the synthesis of these polyphenolics is based on the same substrate, which will inevitably lead to the existence of competition. A lightinduced metabolism analysis in Arabidopsis roots appeared to support this inference (Hemm et al., 2004).

CONCLUSION
This research is the first insight into the effects of elevation and plant age on phenylpropanoid biosynthesis in wild Ginkgo trees, based on transcriptomics and metabolomics analysis.
The key genes involved in flavonoid biosynthesis, including DFR, LAR, and ANR were significantly up-regulated with the elevation increase. Mass spectrometry analysis indicated a similar potential trend that the content of most flavonoid derivatives increased with the elevation rise. The light fluctuation caused by elevation increase may be the main impetus in this process. DFRs appeared to exhibit the most sensitive responses in the face of environmental changes; therefore, might be the main regulatory point in the flavonoid synthesis. Although the age variability only affected the expression of related genes to a lesser extent, metabolites analysis speculated that the growth of Ginkgo age might be detrimental to flavonoid synthesis in leaves. In addition, multiple genes implicated in the synthesis of HCAs, lignin, and lignan exhibited fluctuations with the elevation changes. Notably, significantly up-regulated CADs and down-regulated PRDs potentially led to the accumulation of p-Coumaryl alcohol, one of the lignin monomers, and might inhibit further lignification, while the age increase did not observe a significant impact. Therefore, flavonoid biosynthesis exhibited a higher degree of perturbation than that of HCAs, lignins and lignans, potentially suggesting that flavonoid biosynthesis might be more susceptible than other branch pathways involved in phenylpropanoid biosynthesis. However, moderate environmental changes may mask fluctuations in insensitive components in this study. And the effect of Ginkgo age on the synthesis of phenylpropanoids in leaves needs further excavation. In conclusion, these findings can facilitate the functional genomic and metabolic research in Ginkgo, and provide a reference for artificial planting and industrial harvesting.