Metabolome profiling in different color period samples. Throughout the leaf color change cycle of A. mandshuricum, it undergoes a transition from green leaves (GL) to half red leaves (HRL) and then to red leaves (RL) (Fig. 1). Therefore, to detect the change in metabolite concentration (mainly focused on anthocyanins) associated with the leaf color change process in A. mandshuricum, three types of leaves during different leaf color change periods were used for experiments. We profiled the metabolome of the three samples using UPLC (ultraperformance liquid chromatography) and a widely targeted metabolomics approach. Then, 26 compounds involving anthocyanins were detected and grouped into 8 classes (Table S1 and Fig. 2A), including 5 cyanidins, 4 pelargonidins, 4 peonidins, 4 malvidins, 4 proanthocya-nidins, 2 delphinidins, 2 flavonols and 1 petunidin. Subsequently, the metabolite concentration data were used to perform a hier-archical heatmap clustering analysis of those samples. Among them, most metabolites, such as cyanidins (cyanidin 3,5-O-diglucoside, cyanidin 3-O-(6-O-malonyl-beta-D-glucoside), cyanidin 3-O-arabinoside, cyanidin 3-O-glucoside, cyanidin 3-O-rutinoside), pelargonidins (cyanidin 3-O-(6-O-malonyl-beta-D-glucoside), cyanidin 3-O-arabinoside, cyanidin 3-O-glucoside, cyanidin 3-O-rutinoside), peonidins (peonidin 3,5-O-diglucoside and peonidin 3-O-glucoside), proanthocyanidins (procyanidin B1, procyanidin B2, procyanidin B3, procyanidin C1), flavonols (quercetin 3-O-glucoside, fla_dihydrokaempferol) and petunidin 3-O-glucoside were identified as significantly upregulated during the red leaf period (RL) of A. mandshuricum (Fig. 2B). From Fig. S1, all the biological duplication groups were successfully grouped together, which indicated that the metabolic group data were highly reliable. In addition, it was obvious that the 9 samples were separated into two classes at the metabolite level, green class (GL) and deep-colored (HRL and RL) samples, and there were great differences in the patterns of metabolite accumulation between the two classes (Fig. S1).
To identify differential metabolites differentially accumulated metabolites (DAMs) between the different samples (GL vs HRL, GL vs HRL_RL, GL vs RL, HRL vs RL) in A. mandshuricum, the metabolites with VIP ≥ 1 and fold change ≥ 2 or ≤ 0.5 were selected, and the heatmap and Venn diagram were used to analyze the differences and expression levels of different metabolites more intuitively and quickly. The results showed that there were significant differences in metabolites between the groups in the three leaf color changes, including 8 (8 up-regulated), 8 (8 up-regulated), 10 (10 up-regulated) and 7 (7 up-regulated) DAMs in GL vs HRL, GL vs HRL, GL vs HRL-RL, GL vs RL, and HRL vs RL, respectively (Table S2). It is worth mentioning that DAMs in the anthocyanin biosynthesis pathway of the four comparative groups were significantly increased (Fig. 2C-F), so we proposed that they were related to the leaf color changes of A. mandshuricum, especially cyanidin 3,5-O-diglucoside, cyanidin 3-O-(6-O-malonyl-beta-D-glucoside), cyanidin 3-O-arabinoside, cyanidin 3-O-glucoside, pelargonidin 3-O-glucoside, peonidin 3-O-glucoside, and petunidin 3-O-glucoside. In the Venn diagram (Fig. 2G), 5 anthocyanins were the common metabolites among green-colored and deep-colored samples, and all of them were upregulated and enriched significantly in the anthocyanin biosynthesis pathway (ko00942) (Table S3). Therefore, we speculated that those DAMs produced in the biosynthetic pathways of anthocyanins might be key metabolites during the leaf color change periods of A. mandshuricum.
Analysis of transcriptome data. Throughout the leaf color change cycle of A. mandshuricum, Illumina HiSeq 2000 was used to construct and sequence nine cDNA libraries based on the total RNAs of green (GL), half-red (HRL) and red (RL) leaf samples of A. mandshuricum. A total of 378,613,872 raw reads and 368,399,930 clean reads were obtained, and in the clean reads, the average Q20, Q30 and GC content was 97.63%, 93.42% and 43.89%, respectively (Table S4). Then, gene function annotation was carried out, and we used BLAST software to compare unigenes to the KEGG, NR, Swiss-Prot, GO, KOG, Trembl and Pfam databases. A total of 201,642 unigenes were successfully annotated, of which 150946 unigenes (74.86%) were annotated in at least one database above (Table S5 and Fig. 3B). Among the 20162 unigenes, 96,552 (47.88%) unigenes were grouped into 25 KOG classifications; the most highly represented was group R (general function prediction only), and group C (energy production and conversion), group J (translation, ribosomal structure and biogenesis), group O (posttranslational modifi-cation, protein turnover, chaperones) and T (signal transduction mechanisms) also shared a high percentage of genes (Fig. 3C). A total of 116,336 unigenes were categorized into three GO categories (biological process, cellular component and molecular function) and further divided into 60 major functional groups (Fig. 3D). To reflect the expression level of transcripts, FPKM was used as an indicator to measure the level of gene expression. Principal component analysis (PCA) of the samples showed that the biological replications of GL samples and RL samples were obviously clustered together and significantly separated from each other, but the three HRL samples were dispersed (Fig. 3A). Similar to the metabolomic analysis, there was a significant separation of samples with different colors, indicating that changes in the accumulation of metabolites from green to red leaves were closely controlled by differences in gene expression.
Differentially expressed genes between different colored leaves. DESeq226 was used to analyze the differential expression between different sample groups, and the results showed that a total of 6988 (3106 down-regulated and 3882 up-regulated), 18,597 (11,442 down-regulated and 7155 up-regulated), and 11,474 (9433 down-regulated and 2041 up-regulated) DEGs were obtained in the GL vs HRL, GL vs RL, and HRL vs RL groups, respectively (Table S6). The differentially expressed genes (DEGs) were identified, and the results are shown in Fig. 4A. From the Venn diagram, 217 genes were the common DEGs in all comparison groups, suggesting that those genes across the three compared groups might be associated with the leaf color change of A. mandshuricum (Fig. 4B). Next, we performed GO and KEGG functional annotation and pathway analysis in the three groups and found that the top enriched KEGG terms in the GL vs HRL group were ko01100 (metabolic pathways), ko01110 (biosynthesis of secondary metabolites), ko03010 (ribosome), ko04626 (plant-pathogen interaction) and ko04075 (plant hormone signal transduction) (Fig. 4D). Among them, the metabolic pathways and biosynthesis of secondary metabolites were also abundant in the GL vs RL and HRL vs RL groups, indicating that these DEGs might play a key role in leaf color change pathways (Table S7). GO functional enrichment analysis of DEGs between GL vs RL and GL vs HRL revealed that flavonoid metabolic process (GO: 0009812) and flavonoid biosynthetic process (GO:0009813) were significantly enriched (Fig. 4C).
To gain insight into the dynamism of gene expression pattern changes over the leaf color change periods, k-means cluster analysis was performed using the FPKM value of GL, HRL, and RL transcriptome data. Six clusters of co-expressed genes with different expression trends during leaves of different colors were separated by cluster analysis, implying that each cluster may have similar molecular functions (Fig. 5A). Class 1 consisted of 6863 genes, and the expression trend first increased slowly and then decreased rapidly from green to red leaves. By contrast, class 3 contained 3169 genes with an opposite trend compared with class 1. Class 2 encompassed 3479 genes displaying an increase in expression levels from GL to HRL but a decrease from HRL to RL; this trend was similar to class 4 (1752 genes). A total of 6541 genes were identified in class 5, displaying weak decreases from GL to RL. The most special gene was class 6, which had 2899 genes and exhibited a significant increasing trend for expression level, indicating that they were positively associated with the color change from green to red in A. mandshuricum. Several transcription factors (TFs) were encoded by these genes in class 6, predominantly NAC, MYB-related, RWP-RK, GARP-G2-like, bHLH and C2H2, which might play significant roles in the regulation of gene expression levels involved in leaf coloration of A. mandshuricum (Fig. 5B).
Modulation of phenylpropanoid, flavonoid, and anthocyanidin biosynthesis pathway genes during leaf color change. Anthocyanin components determine the color of A. mandshuricum leaves. However, the key genes involved in A. mandshuricum anthocyanin synthesis have not been reported. In this study, we constructed phenylpropanoid, flavonoid, and anthocyanidin biosynthesis pathways in A. mandshuricum. Ninety DEGs were searched based on enrichment analysis (Table S8) and predicted the molecular mechanisms leading to red leaf coloration (Fig. 6). From the expression levels, most key pathway genes exhibited significant changes. Indeed, encoding PAL (11 DEGs), C4H (8 DEGs), 4CL (10 DEGs), CHS (4 DEGs), CHI (9 DEGs), DFR (9 DEGs), F3H (6 DEGs), F3’5’H (1 DEG), ANS (2 DEGs) and BZ1 (2 DEGs) enzymes were more highly expressed in HRL and RL than in GL, denoting that those DEGs were the key genes involved in the red pigmentation of A. mandshuricum leaves. Conversely, the expression level of one PAL gene (Cluster-20561.130819) was significantly down-regulated from the phenylalanine to cinnamic acid route when the leaves changed from green to red. Similarly, two 4CL genes (Cluster-20561.65497 and Cluster-20561.89667), one CHS gene (Cluster-20561.139874), one CHI gene (Cluster-20561.110302), four F3H genes (Cluster-20561.105759, Clus-ter-20561.110898, Cluster-20564.154341 and Cluster-20561.94199), one ANS gene (Cluster-20561.91533), two BZ1 genes (Cluster-20561.143.657 and Cluster-20561.81061) and one UGT75C1 gene (Cluster-20561.66602) were also significantly down-regulated in their key pathways, implying that silencing these genes was essential for the leaf change to red. In the meantime, some notable genes were also discovered, two C4H genes (Cluster-20561.26608 and Cluster-20561.35599), one 4CL gene (Cluster-20561.36167) and one BZ1 gene (Cluster-20561.99239) were activated in HRL stages, but with low expression in GL and RL. The discovery of these key differentially-expressed genes suggested that they could play a pivotal role in the process of A. mandshuricum leaf discoloration.
Comparing metabolite accumulation in anthocyanin synthesis pathways (ko00941) (Fig. 7), five cyanidins (up-regulated), two delphinidins (one up-regulated and one down-regulated), four malvidins (two up-regulated and two down-regulated), one pelnidin (up-regulated), four pelargonidins (up-regulated), and four peonidins (up-regulated) were determined. Thus, in summary, our results confirmed that the key genes identified in the biosynthesis pathways of flavonoids and anthocyanins could play a key role in A. mandshuricum leaves changing to red. Existing research has confirmed that MYB-related MYB, bHLH, bZIP and other transcription factors (TFs) are widely involved in the regulation of phenyl propane metabolic pathways as regulatory proteins27,28. Analysis of the expression levels of TFs such as bHLH, bZIP and MYB showed that some of them were down-regulated over the red periods, whereas others were significantly up-regulated, implying that these TFs positively regulate flavonoid and anthocyanidin biosynthesis pathway genes (Table S9). Overall, these TFs bHLH, bZIP, MYB and more coordinately regulate those key genes, leading to the accumulation of anthocyanin metabolites (cyanidin 3-O-arabinoside, cyanidin 3,5-O-diglucoside, cyanidin 3-O-glucoside, cyanidin 3-O-rutinoside, cyanidin 3-O-(6-O-malonyl-beta-D-glucoside), delphinidin 3-O-glucoside, malvidin 3,5-diglucoside, malvidin 3-O-rutinoside, pelnidin 3,5-O-diglucoside, pelargonidin 3-O-(6-O-malonyl-beta-D-glucoside), pelargonidin 3-O-galactoside, pelargonidin 3-O-glucoside, pelargonidin 3-O-rutinoside, peonidin 3,5-O-diglucoside, peonidin 3-O-galactoside, peonidin 3-O-glucoside, petunidin 3-O-glucoside), giving rise to the beautiful red coloration of A. mandshuricum leaves.
Analysis of transcriptome factors. In the biosynthesis process of anthocyanins, it is coregulated by many transcription factors29,30. In this study, we listed anthocyanidin biosynthesis-related TFs in A. mandshuricum, and a total of 341 (106 up-regulated and 235 down-regulated), 454 (103 up-regulated and 351 down-regulated) and 981 (346 up-regulated and 635 down-regulated) differentially expressed TFs were identified in GL vs HRL, HRL vs RL and GL vs RL, respectively (Table S10). Among them, most represented TFs were annotated as those encoding C2H2, zn-clus, bZIP, MYB-related, NAC, bHLH, C3H, SNF2, TRAF, WRKY, C2C2-GATA and MYB in the GL vs RL group. The most abundant TFs in GL vs HRL (46 up-regulated and 78 down-regulated) were C2H2, NAC, bHLH, MYB-related, C3H, B3-ARF, WRKY and MYB TFs. In the HRL vs RL comparison (9 up-regulated and 151 down-regulated), C2H2, AP2/ERF-ERF, bHLH, bZIP, C2C2-GATA, C3H, HMG, MYB and MYB-related TFs were the most annotated. Previous studies have shown that transcription factors, such as MYB, bHLH and bZIP, have a deeper effect on anthocyanin biosynthesis31,32, and we also found them in many transcription factors in A. mandshuricum. For example, most of the MYB TFs (such as Cluster-20561.80755, Cluster-20561.81808, Cluster-20561.110291, Cluster-20561.85332 and Cluster-20561.98223) were co-upregulated with most of the key genes in flavonoid biosynthesis during the discoloration of the leaves in A. mandshuricum, and the similar co-upregulated phenomenon was also found in bHLH (Cluster-20561.117348, Cluster-20561.26518, Clus-ter-20561.46232, Cluster-20561.92650, Cluster-20561.107315, and more) and bZIP (Cluster-20561.89129, Clus-ter-20561.98504, Cluster-20561.102544, Cluster-20561.141195, Cluster-20561.92791 and so on) TFs, indicating that those MYBs, bHLHs and bZIPs were important positive regulators in leaf color change periods of A. mandshuricum. In addition, it is worth noting that in the transcription factors of the three comparison groups, most bHLHs and bZIPs were downregulated when the leaves changed from green to red, indicating that some of the down-regulated genes mediated negative regulation in A. mandshuricum anthocyanin biosynthesis (Table S9).
Conjoint analysis between transcriptome and metabolome. To better understand the relationship between genes and metabolites, we mapped the DEGs and DAMs to the KEGG path map at the same time. From the results, many DEGs and DAMs were significantly enriched in the anthocyanin biosynthesis pathway. At the same time, we selected those results in which the Pearson correlation coefficient was greater than 0.8. A total of 14310 DEGs associated with 26 metabolites were found, and most of them jointly regulated the synthesis of one or more metabolites. They were significantly associated with Cya-3,5-O-diglu (6696 genes), Cya-3-O-(6-O-malonyl)-glu (6807 genes), Cya-3-O-ara (6884 genes), Cya-3-O-glu (6568 genes), Cya-3-O-rut (6872 genes), Del-3-O-rut (3789 genes), fla_dihydrokaempferol (6622 genes), fla_Quercetin-glu (6010 genes), Mal-3-O-(6-O-malonyl)-glu (3199 genes), Pel-3,5-O-diglu (6899 genes), Pel-3-O-(6-O-malonyl)-glu (6769 genes), Pel-3-O-glu (6821 genes), Peo-3,5-O-diglu (6755 genes), Peo-3-O-glu (6827 genes), Pet-3-O-glu (6536 genes) and others (Table S11).
From the results of Conjoint analysis between these DEGs and their correlated metabolites, we found two key structural genes ANS (Cluster-20561.86285) and BZ1 (Cluster-20561.99238) in anthocyanidin biosynthesis pathway were significantly up-regulated when the leaves change to red, indicating that they might enhance accumulation of the predominant metabolite cyanidin 3-O-glucoside which is their downstream metabolite, and contributed the red formation of A. mandshuricum leaves (Fig. 6).
qRT-PCR Validation of DEGs in transcriptome data. To further verify the reliability and accuracy of the transcriptome data, we initially screened 12 DEGs (Cluster-20561.96219, Cluster-20561.85494, Cluster-20561.85765, Cluster-20561.85766, Cluster-20561.102785, Cluster-20561.65565, Clus-ter-20561.86285, Cluster-20561.99238 and Cluster-20561.66602) which associated with phenylpropanoid, flavonoid, and anthocyanidin biosynthesis to detect expression levels in GL, HRL and RL with qRT-PCR. The relative expression trends of these candidate genes were similar to those of RNA-seq, indicating the high reliability and accuracy of transcriptome data (Fig. 8). In line with our expectations, the most pivotal genes ANS (Cluster-20561.86285) and BZ1 (Cluster-20561.99238) in anthocyanin biosynthesis of A. mandshuricum, were highly expressed when leaves change from green to red, suggesting that they are the candidate key genes underlying the leaf color formation in A. mandshuricum.