Metabolome and Transcriptome Analyses of Red-colored Leaves From Eucommia Ulmoides ‘Huazhong No. 12’


 Background: Eucommia ulmoides ‘Oliver’ is an economically important tree species with highly medicinal and ecological values that is naturally distributed in China. E. ulmoides ‘Huazhong No. 12’ (H12) is the only red-leaf genotype of this species. In this study, the pigment contents of H12 and E. ulmoides ‘Huazhong No. 11’ (H11, green leaves) were determined. The differential metabolites in H12 and H11 were detected by UPLC-MS/MS, and the differentially expressed genes were screened by transcriptome. Then the key metabolites and corresponding gene regulation in anthocyanin related metabolic pathway were analyzed.Results: The chlorophyll a, chlorophyll b and carotenoid contents in H12 leaves were lower than H11, while the total anthocyanin content in H12 was 4.06 times higher than in H11. There were 96 up-regulated metabolites in H12, including anthocyanin, proanthocyanidins, flavonoid and flavonol. Among them, four differentially expressed anthocyanins were identified. A total of 8,368 differentially expressed genes were selected from transcriptome between H12 and H11. The flavone and flavonol biosynthesis pathway, anthocyanin pathway and photosynthetic pathway were analyzed. Finally, EuCHI, EuF3'H, EuF3'5'H, EuDFR and Eu3MaT1 were recommended as the key genes. The cyanidin, cyanidin 3-malonyl-glucoside and cyanidin 3, 5-glucoside were responsible for the H12 red leaves.Conclusion: This study revealed the metabolites and gene regulation of anthocyanin synthesis, and the potential function between the anthocyanin and photosynthetic gene expression in E. ulmoides red leaves. Notably, the results also provided a reference for the study of other plant leaf-coloring mechanisms.


Background
Anthocyanins are natural pigments in plant tissues. The plant leaves are red-purple because they are rich in pelargonidin and cyanidin, the basic unit of anthocyanin [1]. And anthocyanins are thought to provide photoprotection under stressful conditions [2]. On the one hand, red leaves had higher anthocyanin content, ltering out a portion of visible light. In addition, anthocyanin absorbed blue-green waveband to prevent excessive light damage to leaves. On the other hand, anthocyanin was thought to be able to eliminate the accumulation of reactive oxygen species (ROS) as strong antioxidant [3]. Therefore, the study of anthocyanin synthesis in leaves is of great signi cance for the analysis of plant ecological adaptability.
The anthocyanin metabolism and molecular regulation of the anthocyanin biosynthetic pathway have been well described in the fruit and owers of several plants, which share the same upstream enzymatic steps with proanthocyanidin and avonoid biosynthetic pathways [4][5][6]. In the pink owers of anthocyanin-rich tea (Camellia sinensis), cyanidin 3-O-glucoside and petunidin 3-O-glucoside were found to have direct impacts on anthocyanin biosynthesis via metabolome [7]. Peonidin, cyanidin and their glycoside derivatives were con rmed to be the main anthocyanins in the peels of purple asparagus [8].
Anthocyanin synthesis is involved in avonoid-related pathways.
Non-photochemical quenching (NPQ) is an important photoprotective mechanism which is activated in few seconds upon exposure to excess light. The PsbS protein of photosystem could dissipate excessive energy and maintain the e cient operation of photosynthesis [14]. However, this process produced  [2]. So, what was the relationship between anthocyanin and photoprotective function? Xu proposed a model that anthocyanin could eliminate excessive ROS and maintain photosynthetic capacity [3]. Moreover, the leaves with high content of anthocyanin and avonoids showed strong antioxidant capacity in Arabidopsis thaliana [15]. Hughes believed that anthocyanin in winter leaves act as light lters and may act as antioxidant pools [16].
However, Zhang hold an opinion that anthocyanin was mainly used as light lter rather than antioxidant in the leaves of Begonia semper orens [2].
Eucommia ulmoides is native to China and has been extensively cultivated across 27 provinces (24°50′N-41°50′N, 76°00′E-126°00′E) because of its high economic value and adaptability [17]. E. ulmoides leaves is a potential bioactive material, which contains chlorogenic acid, phenols, terpenes, avonoids [18]. Notably, Eucommia ulmoides 'Huazhong No. 12' is a special variety with red leaves, which expressed red color in summer and autumn. In addition, only palisade tissue was red in leaves. Previously, the study on 6-year-old red-leaves of E. ulmoides showed that anthocyanin was the main reason for the formation of leaf color, and the main factors affecting the content of anthocyanins were PAL enzyme, activity, temperature, pH value and light, while soluble sugar content has no signi cant effect on anthocyanin [19]. And the total avonoids content of red leaves E. ulmoides was signi cantly higher than that of green leaves [20]. Here, we combined the metabolome with the transcriptome to analyze the signi cantly different gene expression levels and substances present in H12 and E. ulmoides 'Huazhong No. 11' (H11, green leaves). In addition to identifying speci c anthocyanins in H12, we also revealed the accumulation of a set of avonoids, anthocyanins, procyanidin A and procyanidin B, accompanied by the regulation of structural genes in the phenylpropanoid and avonoid biosynthetic pathways.

Photosynthetic Pigment and anthocyanin Content of H11 and H12 Leaves
The detection revealed that chlorophyll a, b, a + b and carotenoid contents of the H12 leaves were lower than those of H11 leaves. And the chlorophyll a, b, a + b and carotenoid contents of the H12 were 0.54-, 0.71-, 0.58-, 0.54-fold compared with H11. However, the total anthocyanin content in H12 leaves was 4.06 times higher than that of H11. The results suggested that anthocyanin was responsible for the redness of H12 leaves (Fig. 2  proanthocyanidins are dimeric avonoids, and the contents of proanthocyanidin A1 and A2 in H12 were 4.15 and 89.46 times higher than in H11. The 69 differentially down-regulated metabolites mainly included 6 amino acids, 9 amino acid derivatives, 15 organic acid derivatives lipid-fatty acids and so on.
18 selected differential metabolites are listed in Table 1. In addition, several important mass spectrograms are shown in Figure S1.

KEGG Enrichment Analyses of Differential Metabolites in H11 and H12
The KEGG pathway enrichment analysis revealed that the biosynthesis of avone and avonol, avonoid, anthocyanin, as well as that of amino acids, changed signi cantly in H12 compared with H11. First, in the avone and avonol biosynthesis pathway (ko00941), 17 avone substances were enriched in this pathway and 10 metabolites were highly expressed in H12. Notably, dihydroquercetin was up-regulated, with its content being7.99 times higher in H12 than in H11. Most importantly, in the anthocyanin biosynthesis pathways (ko00942), the contents of cyanidins and cyanidin 3,5-glucoside in H12 were 10.51 and 45.71 times higher than in H11, respectively. Thus, cyanidins and cyanidin 3,5-glucoside appear to be substances involved in the purple-red appearance of H12 leaves. Moreover, in biosynthesis of amino acids pathway (ko01230), the contents of 7 metabolites including tryptophan and glutamine decreased signi cantly. But in the glutathione metabolism pathway (ko00480), the contents of glutathione disul de and glutathione were up-regulated in H12. The pathway maps mentioned above were shown in Figure S2.

Transcriptome Analysis
For the transcriptome sequencing, the average H12 sample produced 594,236,10 clean reads, the base error distribution rate was 0.02%, and the GC content was 46.57%. The average H11 sample produced 60,047,905 clean reads, the base error distribution rate was 0.02%, and the GC content was 46.76%. And the clean data have been submitted to SRA database of NCBI (Accession ID: SRR12569328 and SRR12569327). The sequenced reads were compared with the reference genome of E. ulmoides. The average rates of total mapping for H12 and H11 were 95.53% and 95.36%, respectively. There were 17,503 genes expressed in H12, 17,765 genes in H11, and 16,687 genes were expressed in both leaf types. DEGs met the criteria of log2 (fold change) ≥ 1 and corrected P ≤ 0.005. In total, 8,368 DEGs in H12 vs H11 were identi ed.
The GO analysis annotated 12,893 genes that mapped to the biological process, cell component and molecular functional classes. The COG annotation placed 24,851 genes into 25 COG categories. The KEGG database annotated 11,104 genes, the Swiss-Prot database annotated 20,382 genes, the Pfam database annotated 20,000 genes, and the Nr database annotated 25,233 genes. The number of genes annotated in the Nr database was the highest, followed by the COG database. The DEGs were subjected to a GO enrichment analysis ( Figure S3). The KEGG enrichment of DEGs was shown in Figure S4, which re ected the differential metabolic processes of H11 and H12 during the same period in leaf tissues.

Photosynthesis Biosynthetic Pathways
To gure out the relationship between anthocyanin and photoprotection, the expression of photosynthetic genes in E. ulmoides red leaves were analyzed. Photosynthesis is related to photosynthesis pathways (ko00195) and photosynthesis antenna proteins pathway (ko00196), which had 5 parts, including photosystem II, photosystem I, cytochrome b6/f complex, photosynthethic electron transport and F-type ATPase. They were shown in Figure S5. Thus, the 13 DEGs related to photosynthesis in H11 was signi cantly up-regulated than that in H12. And the 5 genes in H11 is down-regulated than that in H12 in Table 2. This result showed that most of photosynthesis DEGs in red leaves of H12 were down-regulated at that time.
Combined analysis of key metabolites and corresponding genes related to Flavonoid and Anthocyanidin

Candidate Transcription Factors
Among all of the transcripts, a total of 1931 transcripts were annotated as transcription factors (TFs) which were divided into 44 TF families. Among these TF families, 123 MYB transcripts was the largest, followed by 118 MYB-related and 100 bHLHs. And differential expressed transcripts were screen from MYBs, MYB-related family, bHLHs and WD40 TFs. Then, the phylogenetic analysis and alignment were performed to nd the candidate MYBs and bHLHs involved in anthocyanin biosynthesis. Moreover, the matched sequences were picked from published paper and NCBI.
The phylogenetic tree and alignment analysis results in Fig. 4 indicated EUC13670 (up-regulated, 2.05fold) was probably involved in the synthesis of anthocyanin and proanthocyanin [21,22]. Like other anthocyanin MYBs, EUC13670 contained the highly-conserved R2R3 domain and bHLH binding motifs at the N-terminus. The motif [D/E]Lx 2 [R/K] x 3 Lx 6 Lx 3 R is important for interaction with bHLH TFs [23]. In addition, the conservative [K/R]Pxxx[K/T][F/Y] sequence in the C-terminal was part of the signature for MYBs that positively regulated anthocyanin biosynthesis [24]. Intriguingly, EUC13670 was also likely to regulate procyanidins synthesis in E. ulmoides.
The deduced amino acid sequence of EUC06643, EUC01933, EUC18920 and EUC04483 were align with other bHLHs regulating anthocyanin biosynthesis. However, EUC04483 (down-regulated, -1.47-fold) was proved to be bHLH-MYC TF with helix-loop-helix DNA-binding domain and R2R3-MYB transcription factors in N-terminal. As shown in Fig. 5, the box 11, 18 and 13 domain at the N-terminus were key to the interaction of bHLH and MYB [25].

RT-qPCR Veri cation
To verify the RNA-Seq results, we selected 10 genes (2 avonoids, 4 anthocyanin and 3 photosynthetic metabolism pathway genes). The RT-qPCR results showed the up-regulated and down-regulated gene expression levels in the leaves of H11 and H12 were consistent with RNA-Seq (Fig. 6).

Discussion
The metabolite and gene regulation in leaf anthocyanin synthesis Most of the previous studies analyzed the physiological, cytological and molecular causes of leaf color formation. Then metabolome and transcriptome metabolism analysis were popular in the era of multi omics, showing the large-scale metabolites and gene regulation of anthocyanin synthesis in leaves. In this study, 4 anthocyanins, 2 proanthocyadins, 20 avones, 19 avonols, iso avones of E. ulmoides leaves were qualitatively and quantitatively analyzed based on UPLC-MS/MS. And the EuCHI, EuF3'H, EuF3'5'H, EuDFR and Eu3MaT1 were the key genes responsible for the red E. ulmoides leaves via the combined analysis of transcriptome and metabolome. Finally, cyanidin, cyanidin 3-malonyl-glucoside and cyanidin 3, 5-glucoside were important to red leaves of H12. The pattern of avonoid-related metabolites and gene expression was rst constructed in H12. Another similar study showed that the content of 7 anthocyanins in red leaf were higher than that in green leaf, and 19 key MYB regulators were co-expressed with the avonoid-anthocyanin biosynthetic genes in leaf color formation of Lagerstroemia indica cv. Ebony Embers [27]. The changes of avonoid related metabolites, enzyme activity and gene expression in the leaves color development of Cymbidium sinense 'Red Sun' were detected via metabolome. The above results indicate that transcriptome and metabolome have potential in nding the relation of avonoid-related metabolites and gene regulation in leaf color formation.

Structural genes and transcription factors involved in anthocyanin synthesis
The types and content of anthocyanins that were commonly colored could directly affect the color of plant tissue [28]. The combined analysis showed that the F3'5'H was down-regulated while F3'H was upregulated, which resulted in the activation of cyanidin pathway and the blocking of delphinidin pathway. Therefore, the reason for the red leaves of H12 may be the substrate competition between F3'5'H and F3'H. F3'5'H encodes a avonoid 3',5'-hydroxylase belonging to the cyp75a subfamily, while F3'H encodes a avanone 3-hydroxylase belonging to the cyp75b subfamily [29]. F3'5'H catalyzed naringenin and dihydrokaempferol into dihydromyricetin while F3'H catalyzed them into dihydroquercetin. Besides, there were reports that FLS could affect anthocyanin synthesis through substrate competition. FLS changed the ratio of avonol and anthocyanin, and nally changed the ower color of Petunia hybrida [30]. FLS competed with DFR, causing the white ower mutant of Muscari botryoides Mill to loss cyanidin compound compared with the blue-owered [31].
The MYB TF (EUC13670) and bHLH (EUC04483) were supposed to regulated anthocyanin synthesis in H12 leaves. Wang reported R2R3-MYB TF MYB6 of Populus tomentosa was mainly expressed in young leaves and the overexpression of MYB6 in transgenic poplar resulted in the signi cant up-regulation of anthocyanin and procyanidin accumulation [32]. In addition, PpGST together with PpMYB10.1 might regulate the anthocyanin accumulation [33]. The key TFs could regulate the expression of anthocyanin synthesis genes. For example, R2R3-MYB proteins MYB11, MYB12, and MYB111 participate in transcriptional activation of the early biosynthetic genes CHS, CHI and F3H [34], whereas R2R3-MYB TFs PAP1, PAP2, MYB113, MYB114 and TT2 activate the late biosynthesis enzymes in Arabidopsis seedling [35]. Furthermore, AcMYB123/AcbHLH42 complex regulated promoters of AcANS and AcF3GT1 that encode the dedicated enzymes for anthocyanin biosynthesis in the inner pericarp of the red-centered kiwifruit Actinidia chinensis cv. 'Hongyang' [4]. Thus, further experiments are needed to prove the function of EuMYB, EubHLH and EuGST in the process of leaf color formation, and predict their downstream target genes.

Potential photoprotective function of anthocyanins in E. ulmoides
Transcriptome showed that most of photosynthesis DEGs in H12 were down-regulated. In addition, the previous research reported the photosynthetic capacity differences between H11 and H12 leaves. The light compensation point of H12 was 24.00 µmol·m − 2 ·s − 1 while that of H11 was 26.00 µmol·m − 2 ·s − 1 .
And the light saturation point of H12 was 496.00 µmol·m − 2 ·s − 1 while that of H11 was 521.00 µmol·m − 2 ·s − 1 . the maximum net photosynthetic rate of H12 was 16.91µmol·m − 2· s − 1 while that of H11 was 18.50µmol·m − 2 ·s − 1 . Consequently, the capacity of red leaves to utilize strong light was weaker than green leaves in E. ulmoides. However, the maximum quantum yield of PSII (Fv/Fm) of H12 was signi cantly higher than that of H11 [20]. Gould and Vogelmann found that the red leaves had less photosynthesis than green leaves in Quintinia serrata [36]. The study about red leaves of Begonia semper orens showed the Fv/Fm in red leaves was signi cantly higher than that in green leaves during high-light stress. And red leaves were less sensitive to green light than green leaves because anthocyanin acted as light lters [2].
Additionally, anthocyanins were putative to be antioxidant and mediators of ROS-induced signaling cascades [37]. Manetas hold an opinion that anthocyanins could protect tissues from ROS produced by chloroplasts [38]. However, the other study con rmed that under high light stress, the O 2 − production rate and H 2 O 2 content of green leaf were 71.34% and 6.85% higher than those of red leaf in Begonia semper orens. Simultaneously, the non-enzymatic antioxidant activities of the reaction system determined by DPPH were higher in red leaves while the enzyme antioxidant activities in red leaves were lower than those in green leaves [2]. The Arabidopsis leaves with anthocyanins experiment suggested that the main function of anthocyanins was light attenuation rather than antioxidant in photoprotection. Thus, anthocyanin in the palisade tissue of E. ulmoides might play a role of light curtain. And more experiments are needed to supplement, such as the determination of antioxidant capacity and ROS content in H12 and H11 leaves under high light stress, as well as SOD, peroxidase, CAT enzyme activity.
Association analysis indicated that EuCHI, EuF3'H, EuF3'5'H, EuDFR and Eu3MaT1 were key genes and cyanidin, cyanidin 3-malonyl-glucoside and cyanidin 3, 5-glucoside were responsible for the H12 red leaves. Furthermore, that most of photosynthesis DEGs in red leaves of H12 were down-regulated. The results provide a reference for other plants to study the mechanism of leaf coloration at transcriptional and metabolic levels, and analyze the potential relationship between photosynthetic physiological function and gene regulation in red leaf plants.

Plant Materials and Treatments
The grafted seedlings of H12 (red leaves) and H11 (green leaves) were selected, and typical leaves were collected as materials on July 2018 (Fig. 1B)

Leaf Pigment Content Determination
The chlorophyll and carotenoid contents of H12 and H11 leaves were determined using the direct extraction method [39]. The leaves were weighed to 0.1 g (accurate to 0.01 g) after removing the veins and 80% acetone was added to 10 ml. Then, the supernatant was extracted after 72 h in darkness. The absorbance values at 445, 644 and 662 nm were determined using an ultraviolet spectrophotometer (JINGHUA Instruments 752, Shanghai Jinghua Technology Instruments Co., Ltd.). The anthocyanin content was determined. Samples of 0.1 g (accurate to 0.01 g) were ground in 10 ml 1% HCL methanol solution for 4 h and centrifuged for 15 min at 3,500 × g. The supernatant contents were determined at 530 nm and 600 nm. The calculated formula was described in the previous article [40].

Metabolite
The extraction and separation of freeze-dried leave samples, metabolite identi cation and quanti cation were done at MetWare Biotechnology Co., Ltd (Wuhan, China) and described in detail by Fu Wang [41]. The Ultra Performance Liquid Chromatography Mass Spectrometer (UPLC-MS) was connected to an electrospray ionization (ESI)-triple quadrupole-linear ion trap-mass spectrometer /MS system (Q TRAP-MS/MS, Applied Biosystems 6500 Q TRAP) in this experiment. The characteristic ions of each metabolite were screened through the triple quadrupole (QQQ) mass spectrometer. Then the metabolites were qualitatively analyzed based on the Metware database and public database of metabolite information.
The relative contents of corresponding metabolite were represented by the integrals of chromatographic peak area, which were calculated and corrected by MultiQuant (version 3.0.2) [42]. In addition, partial least squares discriminant statistical analysis was used to maximize the metabolite differences between two groups of samples. Variable importance in projection (VIP) ≥ 1 and fold change ≥ 2 or ≤ 0.5 were used as screening criteria [43].

RNA-Seq
Total RNA of each sample was extracted using a SuperPure plant polyRNA rapid extraction kit (GeneAnswer). RNA concentration and purity were determined using a NanoDrop 2000 spectrophotometer (Thremo Fisher Scienti c, Wilmington, U.S.A.) and the Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA), respectively. RNA integrity was validated using 1% agarose gel electrophoresis. Then, RNA-Seq libraries were prepared for sequencing at Majorbio Technology Co., Ltd (Shanghai, China). The Illumina Novaseq 6000 platform was used for the transcriptome sequencing. Raw reads were ltered to obtain high-quality reads. The clean reads were aligned with the reference genome using the software Hisat2 (version 2.1.0) [44]. After each sample was spliced using StringTie software (version 1.3.3b), the expression levels of genes and transcripts were quanti ed using the software RSEM (version 1.3.1), and quantitative index was expressed in Fragments Per Kilobase of exon model per Million mapped fragments (FPKM) [45]. All of the genes were then searched against protein databases, Nr, COG, KEGG and Swiss-Prot.
The TFs were annotated into the PlantTFDB database. The amino acid sequences of the MYB and bHLH TFs used for construct phylogenetic tree and sequence alignment were acquired from the NCBI database (http://www.ncbi.nlm.nih.gov). The phylogenetic tree image was generated using MEGA 7.0 (the neighbor-joining method with 1000 bootstrap replicates was used). The sequences alignment was conducted by DNAMAN version 9.

Real-Time Quantitative PCR (RT-qPCR) Validation
In total, 12 related genes were selected for RT-qPCR, and the primer list is shown in Table S1. The SYBR Premix Ex Taq™ kit (TaKaRa, Japan) was used, and the RT-qPCR expression assays were performed on a CFX96 platform (Bio-Rad, USA). The GADPH gene (EUC05716) was used as the relative quantitative reference gene. Three technical replicates were performed for each sample. The expression levels of the genes in the list were calculated using the 2 −△△Ct method [47]. The mean values ± standard deviations were determined, and an analysis of variance followed by Duncan's new multiple range test were performed for the signi cance analysis. Declarations problems in Henan Province (Grant Number 192102110169). The funding bodies were not involved in the design of study, collection, analysis, or interpretation of data, or manuscript writing.

Availability of data and materials
The transcriptome clean data have been submitted to SRA database of NCBI (Accession ID: SRR12569328 and SRR12569327).
Ethics approval and consent to participate Not applicable.