Integrated transcriptomic and metabolomic analyses reveals anthocyanin biosynthesis in leaf coloration of quinoa (Chenopodium quinoa Willd.)

Background Quinoa leaves demonstrate a diverse array of colors, offering a potential enhancement to landscape aesthetics and the development of leisure-oriented sightseeing agriculture in semi-arid regions. This study utilized integrated transcriptomic and metabolomic analyses to investigate the mechanisms underlying anthocyanin synthesis in both emerald green and pink quinoa leaves. Results Integrated transcriptomic and metabolomic analyses indicated that both flavonoid biosynthesis pathway (ko00941) and anthocyanin biosynthesis pathway (ko00942) were significantly associated with anthocyanin biosynthesis. Differentially expressed genes (DEGs) and differentially accumulated metabolites (DAMs) were analyzed between the two germplasms during different developmental periods. Ten DEGs were verified using qRT-PCR, and the results were consistent with those of the transcriptomic sequencing. The elevated expression of phenylalanine ammonia-lyase (PAL), chalcone synthase (CHS), 4-coumarate CoA ligase (4CL) and Hydroxycinnamoyltransferase (HCT), as well as the reduced expression of flavanone 3-hydroxylase (F3H) and Flavonol synthase (FLS), likely cause pink leaf formation. In addition, bHLH14, WRKY46, and TGA indirectly affected the activities of CHS and 4CL, collectively regulating the levels of cyanidin 3-O-(3’’, 6’’-O-dimalonyl) glucoside and naringenin. The diminished expression of PAL, 4CL, and HCT decreased the formation of cyanidin-3-O-(6”-O-malonyl-2”-O-glucuronyl) glucoside, leading to the emergence of emerald green leaves. Moreover, the lowered expression of TGA and WRKY46 indirectly regulated 4CL activity, serving as another important factor in maintaining the emerald green hue in leaves N1, N2, and N3. Conclusion These findings establish a foundation for elucidating the molecular regulatory mechanisms governing anthocyanin biosynthesis in quinoa leaves, and also provide some theoretical basis for the development of leisure and sightseeing agriculture. Supplementary Information The online version contains supplementary material available at 10.1186/s12870-024-04821-2.

Quinoa, esteemed for its high ornamental value, serves a dual purpose as both a food crop and a contributor to landscape greening, enhancing the development of leisure and sightseeing agriculture [5].Quinoa seedlings initially display emerald green leaves.As quinoa ripens, its leaves gradually develop a range of vibrant colors, such as pink, green, yellow, purple, red, and orange, contributing to the creation of an aesthetically pleasing landscape [6].The diversity in leaf coloration is closely related to anthocyanin biosynthesis.Anthocyanins, a type of flavonoid, are widespread in various plant leaves, fruits, and other tissues and organs and are known for their potent antioxidant properties [7,8].They are soluble in water, methanol, and other solvents [9] and can inhibit cancer cell growth [10], inflammation, and obesity [11], while also promoting vascular health [12], preserving eyesight, and contributing to beauty care.
The synthesis and metabolism of anthocyanins in quinoa is a complex, yet well-defined biosynthetic pathway.Zheng et al. [13] observed significant up-regulation of flavanone 3-hydroxylase (F3H), flavonoid 3'-5' hydroxylase (F3'5'H), and dihydroflavonol reductase (DFR), along with notable down-regulation of flavonol synthase (FLS), in red and purple flowers compared with white flowers of Nicotiana alata.Luo et al. [14] reported lower expression levels of phenylalanine ammonia-lyase (PAL), chalcone synthase (CHS), and chalcone isomerase (CHI) in green shamrock flowers than in purple shamrock flowers.In the case of tea (Camellia sinensis L.), DFR was up-regulated in pink flowers relative to white flowers, whereas FLS displayed higher expression levels in white flowers [15].In addition to these structural genes, several regulatory genes can affect the anthocyanin synthesis pathway.These genes operate indirectly by acting as transcription factors (TFs), thereby modulating the expression of structural genes.At the transcriptional level, the TF family governing anthocyanin biosynthesis predominantly comprises the v-myb avian myeloblastosis viral oncogene homolog (MYB), basic helix-loop-helix (bHLH), and WD40, which form complexes that directly bind to the promoters of anthocyanin biosynthesis structural genes [16].Nesi [17] and Baudry et al. [18] demonstrated that TESTA 2 (TT2) (MYB) can regulate late biosynthesis genes by binding to TT8 (bHLH42) and TTG1 (WD40), resulting in seed coat pigmentation.Yang et al. [19] found the repression of R2R3 MYB TF in white petals compared to pink petals.Meng et al. [20] found that MrMYB44-like1, MrMYB44-like2, and MrMYB44-like3 of the MYB TF family were highly expressed in green leaves of Malus crabapple (Malus Mill.).Jin et al. [21] identified that MYB28.1 and RL1 are associated with purple pigmentation of purple/green leaf.Additionally, Song et al. [22] discovered that MYB2 plays a positive regulatory role in determining the purple trait in purple head Chinese cabbage (Brassica rapa L.).
The molecular regulatory mechanisms governing leaf coloration in quinoa remain unexplored.In this study, germplasm materials exhibiting emerald and pink leaf hues were selected from the germplasm resources of quinoa at the Science and Technology Innovation Service Center in Hebei Province, China.Employing an integrated approach of transcriptomics and metabolomics, DEGs and DAMs were investigated for their association with anthocyanin biosynthesis in different leaf colors.The results of this study are anticipated to provide pivotal insights for comprehensive exploration of the molecular mechanisms underlying anthocyanin biosynthesis in quinoa.

Plant materials
The examined quinoa varieties, M202 and M146, displayed emerald green and pink leaves, respectively.These varieties were cultivated at the Zhangbei Experimental Station of Hebei Agricultural University on May 25, 2021 and May 26, 2022, with row spacing of 40 cm and plant spacing of 18 cm.
Leaf characteristics of M202 and M146 were recorded from 20 to 110 d of growth, with observations at 10-d intervals.The leaves from the middle of each plant were selected for evaluation.Three distinct stages were identified: color change initiation (70 days after germination, 70 DAG), partial color change (90 DAG), and full color change (110 DAG).The stages were labeled as N1, F1, N2, F2, and N3, F3 for M202 and M146, respectively, representing the initial, elongation, and maturity stages of leaf development.Each sample, ranging from 1 to 3 g, was used for three biological replicates.Samples were sealed in plastic bags with 5-10 grains of silica gel and stored at -80℃.

Quantification of total anthocyanin in quinoa leaves
The anthocyanin content was assessed using the method described by Serrano [23].Each 0.5 g sample (N1, F1, N2, F2, N3, and F3) was ground into powder in 0.1 mol/L hydrochloric methanol solution.The mixture was then shielded from light in tin foil and maintained at 4℃ until the tissue had a white appearance.Subsequently, centrifugation at 10,000 g for 10 min was performed, and the absorbance of the supernatant was measured at 530 and 657 nm.The difference in anthocyanin absorbance at 530 nm was calculated using the equation A 530 -0.25×A 657 .This value was normalized to the fresh weight to represent the total anthocyanin content.The total anthocyanin content of each sample was statistically analyzed by hypothesis test.P > 0.05 indicated no significant difference.0.01 < P < 0.05 indicated significant difference and marked with *.P < 0.01 and marked with **.
A657: Absorption value of anthocyanin at 657 nm.m: weight of sample (unit: g).

Transcriptomic sequencing indicated extremely significant difference
RNA was extracted from the samples using the TRIzol precipitation method in Rio [24].Subsequently, a cDNA library was prepared and the library check was completed.The samples were sequenced on the Illumina (http://www.illumina.com)platform with three biological replicates.The basic group mass values for the raw data were determined using Sequencing by Synthesis (SBS) techniques.High-quality clean data were obtained by removing reads containing joints, reads with an N ratio exceeding 10%, and reads with both a mass value Q ≤ 10 and a base number surpassing 50% of the entire read.High-quality clean data were compared to the tetraploid quinoa reference genome (https://www.ncbi.nlm.nih.gov/datasets/taxonomy/63459/) using an efficient HISAT2 [25] comparison system.Subsequently, StringTie [26] was employed for the assembly and quantification of the reads, generating the initial transcript.This transcript was further compared with various databases, including Nr (https://ftp.ncbi.nih.gov/blast/db/),Swiss-Prot (http://www.uniprot.org/),GO (http://www.geneontology.org/), and KEGG (http://www.genome.jp/kegg/)using BLAST.Annotation information was obtained by employing the HMMER [27] software in conjunction with the Pfam database.

Identification of DEGs
The transcript expression levels were normalized for accurate representation.StringTie [28] employed Fragments Per Kilobase Million (FPKM) as a metric to quantify transcript or gene expression levels using a maximum traffic algorithm.DESeq2 [29] software facilitated the analysis of DEGs.DEGs were selected based on the following criteria: |log2(FC)| ≥ 1 and False Discovery Rate (FDR) < 0.05, with Fold Change (FC) indicating differences in gene expression levels.DEGs were classified as up-regulated (log2(FC) > 1) or down-regulated (log2(FC) < 1).R software was used to generate volcano plots, visually representing the DEGs statistics for each comparison group compared to the FDR values.Gene overlap analyses were conducted for the selected DEGs and Venn diagrams were constructed to illustrate the unique and shared genes among the comparison groups.

GO and KEGG enrichment analyses of DEGs
Functional annotation of DEGs was performed using GO and KEGG databases.GO enrichment analysis was conducted using the GO database (http://www.geneontology.org/), which includes three distinct ontologies: molecular function (MF), cellular component (CC), and biological process (BP) [30].KEGG enrichment analysis was conducted using the KEGG (http://www.genome.jp/kegg/).The DEGs were systematically associated with terms in the GO database to calculate the number of genes within each GO terms.Hypergeometric tests were employed to identify the significant enrichment of GO terms within the DEGs in contrast to the entire genome background.In KEGG enrichment analysis, pathways were treated as the units of analysis.The hypergeometric test was applied to identify pathways exhibiting significant enrichment in the DEGs relative to the background gene set.For both GO and KEGG analyses, pathways were considered significantly enriched if they yielded a q-value (Benjamini -Hochberg-corrected p-value) of less than 0.05.

Metabolomic analysis
The metabolites were extracted and isolated according to the method of Huan et al. [31].Metabolites were analyzed using liquid chromatography-tandem mass spectrometry (LC-MS/MS) for both qualitative and quantitative assessments.Original data from the mass spectrometry analysis underwent baseline filtering.Peak identification, filtering, and alignment of each metabolite were performed using the XCMS package (v3.3.2) in R.This process yielded a data matrix of the mass-to-charge ratio (m/z), retention time, and peak intensity.Finally, the filtered data from all platforms were normalized based on the weight of samples used for extraction and according to sample median intensity group-wise, and the resulting data matrices were utilized for further analysis.Metabolite annotations were determined by referencing public mass spectrometry databases such as HMDB (http:// www.hmdb.ca/),Massbank (http://www.massbank.jp/),LipidMaps (http://www.lipidmaps.org), and METLIN [32] (http://metlin.scripps.edu/index.php).
All detected metabolites were annotated using the MetWare database.Metabolomic data were subjected to multivariate statistical analysis using MetaboAnalyst software (version 5.0).Unsupervised dimension reduction and principal component analysis were performed using the R package (http://www.r-project.org/) to visualize distinctions among the sample groups.Additionally, the orthogonal projection to latent structures-discriminant analysis (OPLS-DA) model was rigorously validated using cross-validation and permutation tests.Score and permutation plots were generated to illustrate the differences between comparison groups.Variable importance in projection (VIP) scores were used to describe the overall contribution of each variable to the model.The criteria for identifying DAMs were set at VIP > 1, |log2(FC)| ≥ 1, and P < 0.05.

Integrated transcriptomic and metabolomic analyses
The combined analysis of transcriptomic and metabolomic is based on the standard analysis results obtained from both omics approaches.By analysing the transcriptomics and metabolomics results, the simultaneous mapping of DEGs and DAMs binding onto the KEGG pathway map allows for a better explanation of the transcriptional regulatory mechanisms in metabolic pathways.DEGs and DAMs were subjected to joint analysis using the Pearson correlation coefficient (PCC) method.The threshold of PCC ≥ 0.6 and P < 0.05 was applied for this analysis.PCC ≥ 0.6 indicates a significant correlation between DEGs and DAMs.

Analysis of differentially expressed TFs
The selected DEGs were subjected to BLAST comparisons with the quinoa genome database and were subsequently combined with the PlantTFDB database (http://planttfdb.cbi.pku.edu.cn/) to identify differentially expressed TF genes.These selected TFs, along with their genes associated with anthocyanin biosynthesis, were further analyzed.The expression patterns of these differentially expressed TFs were analyzed according to the FPKM values from the transcriptomic data.

qRT-PCR verification
RNA was extracted and subjected to reverse transcription to synthesize cDNA for qRT-PCR verification.The reaction conditions comprised denaturation at 95℃ for 10 s, annealing at 60℃ for 20 s, and extension at 72℃ for 15 s for a total of 40 cycles.Using ACT7 [33] as an internal reference, primers were designed using Premier5.0software based on the CDS sequence corresponding to the target gene in the sequence library (Supplementary Table 1).Relative expression levels were determined using the 2 −ΔΔCt method [34].Three biological replicates were performed for each qRT-PCR reaction.

Phenotypic observation of leaf color in quinoa
The growth changes in quinoa leaves, spanning from 20 to 110 d, and the phenotypes at maturity were recorded (Fig. 1A).Leaves at 70, 90, and 110 d represented the stages of color transition, half coloration, and full coloration, respectively.The colorimeter readings for M202 revealed L values between 38 and 44, a values between − 12 and − 5, and b values between 15 and 22 (Fig. 1B).Moreover, M146 exhibited L values between 42 and 49, a values between 13 and 19, and b values between 4 and 8.Both M202 and M146 displayed darkening of leaf color as the developmental period was extended.The ΔE values, indicating color differences, increased with leaf growth.When N1 was the reference, ΔE values for N2 and N3 were 3.13 and 6.87, respectively.Similarly, using F1 as the standard, ΔE values for F2 and F3 were 2.87 and 5.95, respectively.

Transcriptomic analysis and screening of DEGs
Eighteen cDNA libraries were generated from the leaf samples at stages N1, F1, N2, F2, N3, and F3.Transcriptomic sequences were obtained after the removal of adapters and low-quality reads.A total of 140.52 Gb of high-quality clean data were obtained, with a Q20 base composition percentage exceeding 97.52% and a Q30 base composition percentage exceeding 92.96%.The GC content measured greater than 45.15% (Supplementary Table 2).

GO and KEGG enrichment analyses
GO enrichment analysis revealed the enrichment of GO terms in various comparisons.Specifically, in N1 vs. N2,   2).
Annotation analysis revealed that seven GO terms were closely associated with anthocyanin biosynthesis, primarily within BP and MF categories.For BP, the relevant GO terms included GO:0009812 (flavonoid metabolic process), GO:0010468 (regulation of gene expression), GO:0051553 (flavone biosynthetic process), GO:0009813 (flavonoid biosynthetic process), and GO:0009411 (response to UV).In the MF category, GO:0043169 (cation binding) and GO:0016703 (oxidoreductase activity) were identified.Furthermore, the differential expression of TFs WRKY22, HY5, and TGA was enriched in the GO:0010468 term (Fig. 4).

Analysis of transcription factors
In the TF analysis, 2645 DEGs spanning 55 TF families were identified.Notably, the bHLH, bZIP, and WRKY TF families, which are involved in anthocyanin regulation, were highlighted.Specifically, the bHLH and bZIP TF families exhibited significant enrichment of 201 and 91 DEGs, respectively.The WRKY TF family was enriched in 100 DEGs (Supplementary Fig. 4).Key TFs included MYC2 and bHLH14 of the bHLH family, HY5 and TGA of the bZIP family, and WRKY24 and WRKY46 of the WRKY family.These TFs were found to be expressed to varying degrees during quinoa leaf growth, and they were associated with four metabolic pathways: circadian rhythm-plant (ko04712), plant hormone signal transduction (ko04075), MAPK signaling pathway-plant (ko04016), and plant-pathogen interaction (ko04626) (Supplementary Table 4).

qRT-PCR verification of DEGs
To assess the authenticity and reliability of the transcriptome data and the extent of DEGs.Ten DEGs across the six enriched pathways were selected for validation via qRT-PCR, including 4CL (LOC110730923), C3'H (LOC110717575), CHI (LOC110704458), CHS (110,724,467), CYP75B1 (LOC110700687), F3H (LOC110724781), FG3 (LOC110719441), HCT (LOC110713661), PAL (LOC110724960) and CYP73A (LOC110684642).The qRT-PCR analysis revealed a successive linear down-regulation trend in the expression of 4CL, C3'H, CHI, CHS, F3H, FG3, and PAL from N1 to N3 in green leaves.This pattern was consistent with the findings obtained from transcriptomic analysis.The qRT-PCR analysis revealed a successive up-down regulation of 4CL, CHS, CYP73A, FG3, and HCT genes from F1 to F3 in pink leaves, which corroborated with the transcriptomics results.The qRT-PCR trends of CYP73A in the N2, CYP75B1 and HCT in the N3, CHI and F3H in the F1 and C3'H and PAL in the F2 were inconsistent with transcriptomics results.From F1 to F3, C3'H exhibited a u-shaped trend, and CHI showed a linear downward trend in the transcriptomics data.However, in the qRT-PCR results, C3'H displayed a linear upward trend, and CHI exhibited an inverted "V" trend.From N1 to N3, CYP73A showed down-regulated trend and CYP75B1 showed a V-shaped trend in transcriptomics data.However, in the qRT-PCR results, CYP73A displayed a V-shaped trend, and CYP75B1 showed a down-regulated trend (Fig. 8).In summary, the qRT-PCR verification results achieved the expected purpose of the experiment, although the expression of certain candidate genes deviated from the transcriptomic results at specific developmental stages.In addition to the individual structural genes controlling single physiological responses, gene expression patterns also involve gene interactions, pleiotropism, and multigenetic effects.In addition, gene expression is influenced by various factors such as the content of metabolic initiators, substrate concentration and TFs.

Phenotypic analysis of quinoa leaves
Phenotypic observations and analysis of anthocyanin content revealed higher levels in pink quinoa leaves than in emerald green leaves across all growth stages.This trend was consistent with findings in ornamental kale [35], purple and green perilla leaves [36], and quinoa grains [37], indicating a close relationship between the phenotype and anthocyanin content.Color difference analysis showed that ΔE values between adjacent periods were all below 4, indicating imperceptible color changes to the naked eye.However, at 70 DAG (N1, F1) and 110 DAG (N3, F3), the ΔE values exceeded 6, indicating noticeable color changes.Previous studies have demonstrated the accumulation of anthocyanins during senescence [38].According to the anthocyanin content measurements in this study, as quinoa leaves matured after discoloration, there was a progressive increase in anthocyanin content, resulting in a slight deepening of color at 110 DAG.This trend in anthocyanin content variation was consistent with the observed color difference.

Analysis of key DEGs
4CL is a pivotal regulator upstream of the flavonoid metabolic pathway, which regulates anthocyanin and flavonoid synthesis.4CL exhibits notably varied substrate preferences and specificities [39].CHI is the second key enzyme in flavonoid synthesis, and its up-regulation can increase the synthesis of flavonoid secondary metabolites in plants and catalyze intramolecular cyclization [40].Zhang et al. [41] demonstrated a positive correlation between the expression of 4CL and the formation of yellow peel coloration of melon.Ho et al. [42] demonstrated that low expression of CHI and F3H can promote hyperpigmentation.In this study, 4CL was significantly downregulated in the N3 period and significantly up-regulated in the F3 period, which indicating that the high expression of 4CL might promote the formation of pink leaves of quinoa, while the low expression might maintain the green leaves of quinoa.This phenomenon indicates that 4CL has substrate preference, which was consistent with the result of Zhang et al. 's study [41].Conversely, CHI expression displayed a consistent reduction across all periods, which can promote hyperpigmentation.CHI positively correlating with bHLH14.We speculated that bHLH14 regulates CHI during the discoloration phase and affect the color change of quinoa leaves.The precise regulatory mechanisms require further investigation.
Fukusaki et al. [43] employed RNAi to suppress CHS in Torenia lybrida, resulting in transgenic plants with white and grey flowers.Similar successes have been documented in the leaves of Gerbera hybrida [44], Petunia [45], and Gentian corolla [46].These results consistently indicated a positive correlation between CHS expression and vividness of plant color.In this study, CHS expression of pink quinoa leaves was significantly higher than emerald green leaves.This consistent with the results of previous studies.Hinderer et al. [47] demonstrated that CHS activity in carrots was inhibited by naringin and chalcone, influenced by downstream flavonoids in the metabolic pathway.In this study, the N3 vs. F3 comparison group displayed notable differential naringin expression for the first time.Moreover, CHS expression in the N3 and F3 periods was significantly lower than that in the N1, N2, F1 and F2 periods, implying that the increased naringin expression in quinoa leaves also exerted inhibitory effects on CHS activity.
PAL predominantly facilitates the catalysis of trans-cinnamic acid production from L-phenylalanine ammonia, serving as the pivotal entry-and rate-limiting enzyme within the phenylpropane pathway in plants [48].Xue et al. [49] reported the up-regulated expression of PAL can promote the pink testa in peanut.Mattus-Araya et al.  [50] reported high expression of PAL can speeds up the progress of color in developing F. chiloensis Fruit.In this study, PAL expression was significantly down-regulated during the N2 and N3 stages, with expression levels notably higher in the F3 stage than in N3 and s in pink quinoa leaves.This indicates that PAL may play a role in maintaining the persistent pink hue observed in quinoa leaves, whereas its down-regulated may contribute to the maintenance of green coloration.
F3H plays a pivotal role in catalyzing the hydroxylation of flavonoids, which is a crucial step in anthocyanin biosynthesis [51].Studies on alfalfa and carnations [52,53] have indicated that reduced F3H expression can lead to deeper coloration.In this study, both transcriptomic analysis and qRT-PCR validation consistently demonstrated a gradual decrease in F3H expression across distinct leaf growth periods within the same quinoa variety.This trend was consistent with the results of chromatic value in quinoa leaves.The results of this study are consistent with those of previous studies in this field and indicated that the low expression of F3H could promote the deepening of leaves color.However, the expression level of F3H in pink quinoa leaves was significantly higher than that in green quinoa leaves during the same period, suggesting that F3H and other genes jointly regulate the change of leaf color.The specific regulatory mechanism needs to be further study.
HCT has been identified as a key regulator of anthocyanin biosynthesis in quinoa leaves.HCT facilitates the conversion of 4-coumaroyl CoA into caffeoyl-CoA.In addition, Liu [35] elucidated the branching pathway of anthocyanin biosynthesis by inducing concurrent alterations in anthocyanin, chlorophyll, and carotenoid contents, resulting in the emergence of green spots within ornamental kale exhibiting pink leaves.In this way, HCT transforms 4-coumaroyl CoA into caffeoyl-CoA, which subsequently undergoes synthesis through CHS, F3H, DFR, and anthocyanidin synthase (ANS) to yield anthocyanins.This established anthocyanin synthesis pathway is consistent with our current findings.Studies have indicated that silencing HCT not only suppresses lignin biosynthesis but also leads to activation of CHS, resulting in the accumulation of various flavonol glycosides and acylated anthocyanins [54].Li et al. [55] showed that HCT exhibited substrate preference.Transcriptomic analysis in this study revealed continuous expression of HCT in emerald green leaves, with another up-regulation observed during the F3 period.This suggests that elevated HCT expression may facilitate the formation of pink leaves and promote substrate preference.Furthermore, in cultivar M146, there appears to be a propensity for caffeoyl-CoA formation, subsequently leading to the synthesis of various anthocyanins.However, in this study, no discernible relationship between HCT and CHS was observed, warranting further investigation to elucidate the specific regulatory mechanisms governing color change.
FLS catalyzes the conversion of dihydroflavonols to flavonols, a process that competes with DFR at a crucial branch point in the anthocyanin pathway.Ueyama et al. [56] demonstrated that up-regulation of FLS resulted in the production of white flowers in transgenic Nierembergia sp.In this study, significant down-regulation of FLS expression was observed in the F2 and F3 periods, with no notable difference detected in green quinoa leaves.This suggests that the down-regulated expression of FLS might be related to the formation of pink coloration in quinoa leaves.

Analysis of key DAMs
Metabolites within ko00941 and ko00942 pathways affect quinoa leaf pigmentation.Naringenin, which is derived from naringenin chalcone through the action of CHI, is a crucial intermediary in anthocyanin synthesis.Motallebi et al. [57] demonstrated significance of naringin as a polyphenolic flavonoid compound with notable pharmacological effects in plants.This study revealed a significant differential expression of naringin during the N3 and F3 periods, suggesting its pivotal role in the transition from green to pink pigmentation.Furthermore, although no substantial alterations were observed in cyanidin 3-O-(3'' ,6''-O-dimalonyl) glucoside across different periods within the same variety, distinct variations were observed between different varieties during the same period.This indicated that its high expression may play a key role in the formation of pink pigmentation.Cyanidin-3-O-(6"-O-malonyl-2"-O-glucuronyl) glucoside was significantly down-regulated in all phases of green quinoa leaves, whereas no significant difference was observed in pink quinoa leaves and across different varieties during the same period.We propose that cyanidin-3-O-(6"-Omalonyl-2"-O-glucuronyl) glucoside serves as the key metabolite for maintaining green quinoa leaves.Yifan et al. [58] demonstrated a direct correlation between changes in blueberry surface color and anthocyanin concentration.In this study, we observed significant differences in the content of pelargonidin, delphinidin, and cyanidin derivatives within quinoa leaves.This suggests a potential association between the color transformation of quinoa leaves and the levels of these derivatives.

Analysis of TFs
In this study, six pivotal TFs related to anthocyanin biosynthesis were identified: MYC2 and bHLH14 from the bHLH family; HY5 and TGA from the bZIP family; and WRKY24 and WRKY46 from the WRKY family.Based on transcriptomic-metabolomic data, an anthocyanin biosynthesis pathway specific to quinoa leaf coloration was elucidated (Fig. 9A).Pearson's coefficient analysis was used to determine correlations between 13 differentially expressed TFs and 18 DEGs (Fig. 9B).The results showed a positive correlation between bHLH14, WRKY46 and CHS.The high expression of CHS can enhance the expression of pink coloration in quinoa leaves.Therefore, bHLH14 and WRKY46 may work with CHS to positively regulate the formation of pink quinoa leaves and plays a crucial role in anthocyanin biosynthesis regulation.TGA and WRKY46 exhibited a positive correlation with 4CL, which showed TGA and WRKY46 may regulate the formation of pink quinoa leaves or maintain the green color of quinoa leaves.The TF regulatory mechanisms are intricate, necessitating further in-depth investigation in subsequent studies.Additionally, a comprehensive exploration of the functional roles of the bHLH and WRKY families is warranted.

Conclusions
In summary, fundamental insights into the anthocyanin synthesis pathways, related genes, and metabolites in the variously colored of quinoa were acquired through integrated transcriptomics-metabolomics analyses.The primary findings of this study are outlined below: (1) Integrated transcriptomic and metabolomic analyses indicated that both flavonoid biosynthesis pathway (ko00941) and anthocyanin biosynthesis pathway (ko00942) were significantly associated with anthocyanin biosynthesis.( 2

Fig. 3 Fig. 2
Fig. 3 Venn diagram of DEGs in the nine comparison groups.The Venn diagram represents the overlapping DEGs among the nine comparison groups, and the number represents the unique number of DEGs in each comparison group

Fig. 4 Fig. 6
Fig.4 Sankey and dot plots of DEGs enriched in GO analysis of quinoa leaves.The Sankey map on the left delineates the genes associated with anthocyanin biosynthesis enriched in each GO team.The horizontal coordinate of the dot plot on the right is the ratio of each GO team, the vertical coordinate represents the -log10 (P_value), and the size of each circle represents the number of genes enriched in each GO team

Fig. 7
Fig. 7 Sankey diagram.The left side represents the DEGs associated with the DAMs.The right side represents the DAMs associated with the DEGs

Fig. 8
Fig. 8 qRT-PCR validation of 10 key genes.The left Y-axis represents the relative expression of samples in qRT-PCR, and the right Y-axis represents the FPKM value in transcriptomics.The column represents the relative qRT-PCR expression of each component in each comparison group, corresponding to the left Y-axis.The black circle represents the FPKM value of each component in each comparison group, corresponding to the right Y-axis ) The diminished expression of PAL, 4CL, and HCT decreased the formation of Cyanidin-3-O-(6"-O-malonyl-2"-O-glucuronyl) glucoside, leading to the emergence of emerald green leaves.Moreover, the lowered expression of TGA and WRKY46 indirectly regulated 4CL activity, serving as another important factor in maintaining the emerald green hue in leaves N1, N2, and N3.(3) Cyanidin 3-O-(3'' ,6''-O-dimalonyl) glucoside and naringenin are identified as the key metabolites responsible for the distinct coloration observed in the leaves of