Transcriptome and Metabolome Combined to Analyze Quinoa Grain Quality Differences of Different Colors Cultivars

Quinoa (Chenopodium quinoa Wild.) has attracted considerable attention owing to its unique nutritional, economic, and medicinal values. Meanwhile, quinoa germplasm resources and grain colors are rich and diverse. In this study, we analyzed the composition of primary and secondary metabolites and the content of the grains of four different high-yield quinoa cultivars (black, red, white, and yellow) harvested 42 days after flowering. The grains were subjected to ultra-performance liquid chromatography-tandem mass spectrometry (UPLC-MS/MS) and transcriptome sequencing to identify the differentially expressed genes and metabolites. Analysis of candidate genes regulating the metabolic differences among cultivars found that the metabolite profiles differed between white and black quinoa, and that there were also clear differences between red and yellow quinoa. It also revealed significantly altered amino acid, alkaloid, tannin, phenolic acid, and lipid profiles among the four quinoa cultivars. Six common enrichment pathways, including phenylpropane biosynthesis, amino acid biosynthesis, and ABC transporter, were common to metabolites and genes. Moreover, we identified key genes highly correlated with specific metabolites and clarified the relationship between them. Our results provide theoretical and practical references for breeding novel quinoa cultivars with superior quality, yield, and stress tolerance. Furthermore, these findings introduce an original approach of integrating genomics and transcriptomics for screening target genes that regulate the desirable traits of quinoa grain.


Introduction
Quinoa (Chenopodium quinoa Wild.) is a dicotyledonous plant of the Amaranthaceae family, and is native to the Andes Mountains of South America. Global production and consumption of quinoa has exponentially increased over the past few decades [1,2] and the International Year of Quinoa was celebrated in 2013. Quinoa has attracted consumer attention as a "superfood" and is known for its cold, salt, and drought resistance; it is a C3 crop and is reputed as a "golden grain" [3][4][5]. Quinoa seeds are abundant in amino acids, proteins, vitamins, dietary fiber, carbohydrates, flavonoids, phenols, and saponins [6]. Quinoa has a generally higher protein content than most other cereals and has a complete and uniform essential amino acid profile [7]. The pharmacologically active constituents of quinoa might help lower the risks of cardiovascular and neurological disease, as well as diabetes [8]. Quinoa grain is small, which is traditionally used as a grain or mixed with flour, and is a low starch food. It is used in the preparation of various breads and baked goods. However, no existing quinoa variety is suitable for noodle or pasta preparation [2,9,10]. As quinoa is highly adaptable to widely divergent environments, it is now cultivated in France, Great Britain, the Netherlands, Italy, other European countries, and the United States. Quinoa was successfully introduced to China in 1987 and large-scale quinoa cultivation has been established there. However, the main production area for quinoa remains the Andes, in South America, where a wide range of varieties, genotypes, and wild ancestors have been found during some 5000 years of cultivation [11,12]. Quinoa cultivars differ in terms of their nutritional and functional characteristics. However, continuous quinoa cultivation on the same soil eventually degrades its nutritional characteristics and reduces grain yield [13]. Quinoa is a viable and favorable cereal alternative for countries suffering from food insecurity. It may help meet the objectives of increasing high-quality food production and accommodating the nutritional needs of an expanding global population despite severe climate change. Quinoa is highly nutritional and a rich source of bioactive ingredients with strong nutraceutical and market value. Ongoing research continues to clarify the functions and mechanisms of the abundant bioactive ingredients in quinoa. Nevertheless, quinoa yield and quality are influenced by heredity, environmental conditions, and biological factors. In particular, four quinoa varieties with different colors have differences. Therefore, it is preferable to cultivate high-quality, high-yield quinoa cultivars with stable biological characteristics. In this study, we analyzed the metabolomes and transcriptomes of red, white, yellow, and black quinoa. The results of our investigations could provide a reliable basis for breeding novel quinoa cultivars-including those conducive to noodle and pasta fabrication, and facilitate a better understanding of the factors that regulate quinoa grain quality, yield, and stress tolerance.

Grain Metabolism in Four Quinoa Cultivars
A total of 513 metabolites (Table S1) were detected in four sample groups. These included 74 amino acids and their derivatives, 82 phenolic acids, 42 nucleotides and their derivatives, 14 vitamins, 12 terpenoids, 13 tannins, 38 sugars and sugar alcohols, 38 alkaloids, 62 organic acids, and 105 lipids. There were 16 lignins and coumarins, 14 other metabolites, and three steroids. The principal component analysis (PCA) score and the heat map showed significant differences in grain metabolism among the four quinoa cultivars. The differences in grain metabolism were extensive between the white and other three quinoa cultivars, and narrow between the red and black quinoa cultivars ( Figure 1).  Table 1 shows that the differences in amino acids were smallest between B vs. R groups. They only differed in terms of 3,4-dihydroxy-L-phenylalanine and L-tyramine with log2FC for R vs. B being −7.05 and −1.69, respectively. There were obvious differences in 19 amino acids between B vs. W groups. Of these, L-threonine, L-asparagine, 3,4-dihydroxy-L-phenylalanine, and L-tyramine were downregulated and the maximum log2FC was −9.92, whereas the remaining 15 amino acids were upregulated. There were obvious differences in 19 amino acids between B vs. Y, of which L-pheny lalanyl-L-phenylalanine, homoarginine, and N-acetyl-L-phenylalanine were upregulated, whereas the other 16 amino acids were downregulated. Notably, 3,4-dihydroxy-L-phenylalanine showed the maximum difference among the groups (log2FC = 10.07). There were differences in 11 amino acids between R vs. Y. Nine of the amino acids were downregulated except homoarginine and N-acetyl-L-phenylalanine. There were differences in 20 amino acids between W vs. R. N-acetyl-L-tryptophan and L-threonine were upregulated, whereas all others were downregulated. A total of 28 amino acids differed between W vs. Y. Homoarginine, N-acetyl-L-phenylalanine, and N-acetyl-L-tryptophan were upregulated, whereas the other 25 amino acids were downregulated. Table S2 shows that 16 phenolic acids obviously differed between B vs. R groups; Syringin and feruloylmalic acid were upregulated, whereas all others were downregulated. The log2FC corresponding to p-coumaraldehyde was −16.19. A total of 35 phenolic acids obviously differed between B vs. W. A total of 23 phenolic acids, including 3,4-digalloylshikimic acid, tyrosol, 4-aminosalicylic acid, and especially p-coumaraldehyde (log2FC = −16.19), were obviously downregulated, while the other 12 phenolic acids were upregulated. A total of 38 phenolic acids obviously differed (eight upregulated and thirty downregulated) between B and Y groups. The log2FC of 2-(7-dihydroxyl)benzofuranylferulic acid and feruloyltartaric acid were 18.23 and 18.16, respectively, and those of 4-O-(6 -O-glucosylcafeylglucosyl)-4-hydroxybenzyl alcohol and p-coumarin were −15.33 and −16.19, respectively. There were obvious differences in 39 phenolic acids between R vs. Y, of which 11 phenolic acids, including 2-(7-dihydroxyl)-benzofuranylferulic acid, salvianolic acid A, and 6 -O-ferroloyl-D-sucrose, were upregulated, while 28 phenolic acids, such as p-coumaroylcafeoyltartaric acid, tyrosol, and 4-O-(6 -O-glucosylcaffeoylglucosyl)-4hydroxybenzoyl alcohol, were obviously downregulated. There were obvious differences in 29 phenolic acids between W vs. R groups. A total of 16 phenolic acids, such as 3,4-digalloylshikimic acid, tyrosol, and 4-aminosalicylic acid were upregulated, while 13 phenolic acids, including 4 -hydroxy-3 -methoxyacetophenone, sibiricose A5, and syringic acid-4-O-(6 -feruloyl)glucoside, were obviously downregulated. Moreover, there were obvious differences in 36 phenolic acids between W and Y groups. Of these, 12 were obviously upregulated, including 2-(7-dihydroxyl)-benzofuranylferritic acid, feruloyltartaric acid, and 3,4-diglloylshikimic acid, whereas 24 phenolic acids, such as 4-O-(6 -Oglucosylcaffeoylglucosyl)-4-hydroxybenzoyl alcohol, p-coumaroylcafeoyltartaric acid, and cimicifugic acid K, were obviously downregulated.

Analyses of Nucleotides and Their Derivatives
A total of 42 nucleotides and their derivatives were detected in 12 samples. Table S3 shows the differences in nucleotides and their derivatives among the four quinoa cultivars. There were no differences in the nucleotides and their derivatives between B vs. R groups. However, obvious differences were observed in 11 nucleotides and their derivatives between B vs. Y groups, where N7-methylguanosine and Isoxanthopterin were obviously upregulated, whereas the other three were obviously downregulated. Seven nucleotides and their derivatives obviously differed between R vs. Y groups, and five between W and R groups, all of which were downregulated. However, the difference in xanthine was the largest (log2FC = −6.39). A total of 12 nucleotides and their derivatives differed between W vs. Y groups, and the differences were larger than those between any other group pair. Isoxanthopterin was upregulated, whereas the other 11 were downregulated.
48 Note: Log2FC is the logarithm base 2 of fold change (FC) of the differential metabolite. If log2FC is positive, it means up regulation; if log2FC is negative, it means down regulation. "-" means no difference, the same below.

Analyses of Organic Acid-Related Metabolites
Only L-tartaric acid obviously differed between B vs. R groups, and was down regulated. There were obvious differences in 10 organic acids between B vs. W groups. Of these, seven organic acids, including 4,8-dihydroxyquinoline-2-carboxylic acid, 2-isopropylmalic acid, and 3-isopropylmalic acid, were obviously upregulated, while L-tartaric acid, 2hydroxy-4-methylpentanoic acid, and L-homoserine were obviously downregulated. There were obvious differences in 15 organic acids between B vs. Y groups. Among these, 4,8dihydroxyquinoline-2-carboxylic acid and 2-furanoic acid were obviously upregulated, whereas 13-hydroxybutyric acid, D-xylonic acid, and phosphoenolpyruvate were obviously downregulated. There were obvious differences in 12 organic acids between R vs. Y groups. Among these, 4,8-dihydroxyquinoline-2-carboxylic acid, 2-hydroxyethylphosphonic acid, 2-furanoic acid, and 2-hydroxyhexanoic acid were obviously upregulated, and the rest were obviously downregulated. There were obvious differences in 17 organic acids between W vs. R groups, of which L-tartaric acid, 2-furanoic acid, 2-hydroxy-4-methylpentanoic acid, and L-homoserine were obviously upregulated, while the rest were obviously downregulated. There were obvious differences in 24 organic acids between W vs. Y groups. Among these, L-tartaric acid and 2-furanoic acid were obviously upregulated, and the rest were downregulated (Table S5).

Analyses of Tannin-Related Metabolites
Among these, 1-O-galloyl-D-glucose was obviously upregulated, whereas the other 11 were obviously downregulated. Log2FC was −20.86 for both procyanidin B3 and procyanidin B2. There were obvious differences in 12 tannins between W vs. R groups. Of these, 3-O-methylgallic acid was obviously downregulated, whereas the others were obviously upregulated. Log2FC was 20.86 for both procyanidin B3 and procyanidin B2.

Analyses of Terpenoid-Related Metabolites
A total of 12 terpenoids were detected in the 12 samples. Table S9 shows the differences in terpenoid profiles among the four quinoa cultivars. Only camaldulenic acid obviously differed between B vs. R (log2FC = −8.47) groups. There were obvious differences in six terpenoids between B vs. W groups. Of these, 24,30-dihydroxy-12 (13)

Analyses of Vitamin-Related Metabolites
Table S10 shows the differences in vitamin profiles among the four quinoa cultivars. There were differences in two vitamins between B vs. R groups. L-ascorbic acid was upregulated (log2FC = 3.53), while biotin was obviously downregulated (log2FC = −1.02). There were five obvious differences in the vitamin content between B vs. W groups. Biotin, riboflavin, and L-ascorbic acid were upregulated, whereas nicotinamide and 4-pyridoxic acid were obviously downregulated. There were 11 obvious differences in vitamin content between B vs. Y groups. L-ascorbic acid, pyridoxal, biotin, and nicotinate-D-ribonucleoside were obviously upregulated, while pyridoxine, 4-pyridoxic acid-O-glucoside, and pyridoxine-5 -O-glucoside were obviously downregulated. There were eight obvious differences in vitamin content between R vs. Y groups. Pyridoxal, biotin, and nicotinate-D-ribonucleoside were upregulated, while pyridoxine, 4-pyridoxic acid-O-glucoside (log2FC = −10.95), and dehydroascorbic acid were obviously downregulated. There were five differences in vitamin content between W vs. R groups. Of these, 4pyridoxic acid, pyridoxine, and 4-pyridoxic acid-O-glucoside were obviously upregulated, whereas biotin and riboflavin were downregulated. There were six obvious differences in vitamin content between W vs. Y. Nicotinate-D-ribonucleoside was obviously upregulated, and the rest were obviously downregulated. For 4-pyridoxic acid-O-glucoside, log2FC was −9.52.

Analyses of Other Metabolites
Fourteen other metabolite classes were detected in the 12 samples. Table S11 shows the differences in other metabolite classes among the four quinoa cultivars. There were no obvious differences between B vs. R groups in terms of these metabolites. There were five obvious differences in other metabolites between B vs. W groups. Of these, 5,7-dihydroxychromone and D-(+)-melezitose-O-rhamnoside were obviously upregulated, whereas the rest were downregulated. The types and trends of six other metabolites did not obviously differ between B vs. W, or between R vs. Y groups. Among these, dihydrocharcone-4 -O-glucoside was upregulated, while 3-phospho-D-glycoric acid, N,Ndimethylformamide, palmitaldehyde, 4-methyl-5-thiazolethanol, and phenethylamine were downregulated. There were four obvious differences in other metabolites between W vs. R groups. Of these, dihydrocharcone-4 -O-glucoside and N,N-dimethylformamide were upregulated, whereas 5,7-dihydroxychromone and androsin were downregulated. There were five obvious differences in other metabolites between W vs. Y groups. Of these, dihydrocharcone-4 -O-glucoside was obviously upregulated, while 3-phospho-Dglyceric acid, 5,7-dihydroxychromone, palmitaldehyde, and phenethylamine were obviously downregulated.

Analyses of Correlations between Metabolome and Transcriptome among Different Quinoa Cultivars
In this study, we conducted canonical-correlation analyses (CCA) (Tables S12-S23) and network correlation analyses on differentially expressed genes and metabolites in the enrichment pathways for each comparison group. The enrichment pathways common to all six comparison groups were phenylalanine metabolism, phenylpropanoid biosynthesis, and amino acid, and ABC transporter biosynthesis. According to Figure 2 and Tables 3 and 4, thirty-eight and eight differentially expressed genes and metabolites, respectively, participated in the phenylalanine metabolic pathways of all six comparison groups. The correlation was highest between salicylic acid-2-O-glucoside and aspartate aminotransferase, cytoplasmic [EC:2.6.  Table 3 shows that gene-LOC110705616 expression was higher in black and red quinoa than it was in white and yellow quinoa. Table 4 shows that the relative salicylic acid content was higher in black and red quinoa than it was in white and yellow quinoa. Hence, gene-LOC110705616 induces salicylic acid-2-O-glucoside biosynthesis. Gene-LOC110718568 expression and relative cinnamic acid content followed the order black quinoa > red quinoa > white quinoa > yellow quinoa. Therefore, gene-LOC110718568 promoted cinnamic acid biosynthesis, gene-LOC110738220 inhibited N-acetyl-L-phenylalanine biosynthesis, gene-LOC110732717 and gene-LOC110728031 induced phenethylamine and 2-phenylethanol biosynthesis, gene-LOC110685697 promoted fumaric acid biosynthesis, and gene-LOC110735801 promoted p-hydroxyphenylacetic acid biosynthesis.   In the amino acid biosynthetic pathway, 45 and 18 differentially expressed genes and metabolites, respectively, had PCC > 0.9. Table 5 shows the differentially expressed genes that were the most highly correlated with the differentially expressed metabolites in the amino acid biosynthetic pathway. For phosphoenolpyruvate and gene-LOC110712605, PCC was 0.99, while for N-acetyl-L-glutamic acid, gene-LOC110725557, and gene-loc110683155, PCC was −0.97. Table 6 shows the differentially expressed genes and metabolites with PCC > 0.9 in the ABC transporter pathway. For D-(+)-trehalose and gene-LOC110689873, PCC was 0.96, whereas for lactobiose and gene-LOC110729558, PCC was −0.94. Figure 3 shows differentially expressed genes and metabolites with PCC > 0.8 in the phenylpropanoid biosynthetic pathway. Table S24 shows that in the phenylpropanoid biosynthetic pathway, forty and six differentially expressed genes and metabolites, respectively, had PCC > 0.9. Among these, scopoletin (7-hydroxy-5-methoxycoumarin) had the highest positive correlation with gene-LOC110701713 (PCC = 0.93) and the highest negative correlation with gene-LOC110715013 (PCC = −0.96). PPC was 0.9 between p-coumaric acid and gene-LOC110724195. For trans-5-O-(p-coumaroyl) shikimate and gene-LOC110724195, PCC was 0.92, whereas for trans-5-O-(p-coumaroyl)shikimate and gene-LOC110685879, PCC was −0.91. For caffeic aldehyde and gene-LOC110716610, PCC was 0.99, while that for caffeic aldehyde and gene-LOC110709413 was −1.00. For cinnamic acid and gene-LOC110685747, PCC was 0.92, whereas for cinnamic acid and gene-LOC110732640, PCC was −0.92. For syringin and gene-LOC110721053, PCC was 0.94, while for syringin and gene-LOC110707634, PCC was −0.90. Table 5. Correlation between differential genes and differential metabolites in amino acid biosynthesis pathway.  Table 6. Correlation between differential genes and differential metabolites in the ABC transporter pathway.

The Name of Metabolites Gene Name PCC
Lactobiose

Discussion
Quinoa germplasm resources and grain colors are rich and diverse. The colors of certain quinoa cultivars change during grain maturation; therefore, agronomic and quality traits widely differ among quinoa germplasms [1,14]. Meanwhile, there is broad diversity among quinoa cultivar in terms of antioxidant levels [15], nutritional characteristics [16], and dietary fiber content [17], and tannins containing good antibacterial and antioxidant effects [18]. Studies have shown that quinoa has antioxidant and bacteriostatic properties, lowers blood lipid levels, improves blood glucose, and could serve as a sustainable source of dietary supplements and functional components [19,20]. In this study, we selected seeds from three independent plants at the same growth stage and 42 d after flowering among the four cultivars (red, white, black, and yellow quinoa cultivars). After analyzing the metabolites of four different quinoa cultivars by UPLC-MS/MS, we determined the four quinoa cultivars significantly differed in terms of amino acid, alkaloid, organic acid, and tannin composition and content. Notably, they most obviously differed in terms of tannin content; the difference in tannin content was least between red and black quinoa, and the maximum tannin concentrations in red and black quinoa were far higher than those in white and yellow quinoa. The levels of the ubiquitous plant pigments procyanidin B3 and procyanidin B2 did not differ between red and black quinoa, or between white and yellow quinoa. Nevertheless, they were much higher in red and black quinoa than they were in white and yellow quinoa. Certain amino acids cannot be synthesized by the human body and must be acquired from food; the essential amino acid content is higher in quinoa than it is in rice, wheat, corn, or other cereals. Quinoa is abundant in eight amino acids essential for human health, as well as histidine required by infants and young children. Quinoa also contains lysine, threonine, and tryptophan, which are typically lacking in most other plant protein sources, and lysine is particularly deficient in cereals [21,22]. At the same time, quinoa can be made into soups, stews, biscuits, and various drinks to supplement human needs [23][24][25][26][27]. Here, we observed significant differences in amino acid and derivative composition, and content among the various quinoa cultivars. However, these differences were smallest between red and black quinoa. Our findings suggested that different quinoa cultivars will supplement different amino acids, depending on individual nutritional requirements. Moreover, the present study carefully analyzed the metabolite profiles of four different quinoa cultivars, and also compared the differences in their metabolite composition and content via metabolome and transcriptome correlation analyses. We believe that the results of the present study provide theoretical and practical references for the development and application of quinoa tannin, and that these findings lay theoretical and practical foundations for quinoa product development and enhancement, as well as the cultivation of novel quinoa cultivars with high grain yield and quality, and strong abiotic and biotic stress resistance.

Materials
Red quinoa (R), white quinoa (W), yellow quinoa (Y), and black quinoa (B) cultivars were harvested from plantations in in Xundian County, Kunming, China (102 • 41 E, 25 • 20 N) ( Figure 4). Individual plants flowering on the same day were marked. Seeds from three independent plants at the same growth stage were selected among the four cultivars, immediately frozen in liquid nitrogen 42 d after flowering, and stored at -80 • C until later use. We allocated R, W, Y and B to red quinoa, white quinoa, yellow quinoa, and black quinoa, respectively.

Widely Targeted Metabolome Detection and Analysis
Biological samples were placed in a vacuum freeze-dryer (Scientz-100F, Ningbo Scientz Biotechnology Co. Ltd., Ningbo, China) and ground at 30 Hz for 1.5 min (MM 400, Retsch, Haan, Deutschland) until they were powdered. Subsequently, 100 mg powder was weighed and resuspended in 1.2 mL of 70% (v/v) methanol extract. Each sample was refrigerated at 4 • C and was rotated six times overnight to improve extraction. Each sample was centrifuged at 10,000× g, 4 • C, for 10 min (ANPEL, Shanghai, China, http://www.anpel.com.cn/, accessed on 23 October 2021). Each supernatant was then filtered through a 0.22-µm microporous membrane and stored in a bottle until Ultra Performance Liquid Chromatography, UPLC, (SHIMADZU Nexera X2, https://www.shimadzu.com.cn/, accessed on 23 October 2021) and Tandem mass spectrometry, MS/MS (Applied Biosystems 4500 QTRAP, http://www. appliedbiosystems.com.cn/, accessed on 23 October 2021) (UPLC-MS/MS), analysis. The chromatographic column was an Agilent SB-C18 (1.8 µm; 2.1 mm × 100 mm; Agilent Technologies, Santa Clara, CA, USA). Mobile phase A was ultrapure water with 0.1% (v/v) formic acid. Mobile phase B phase was acetonitrile with 0.1% (v/v) formic acid. The elution gradient was as follows: 0.00 min, 5% mobile phase B; linear increase in mobile phase B to 95% within 9.00 min; 95% mobile phase B for 1 min; linear decrease in mobile phase B from 10.00 min to 11.10 min; and equilibration of mobile phase B to 5% until 14 min. The flow rate was 0.35 mL/min, the column temperature was 40 • C, and the injection volume was 4 µL. The MS conditions were: electrospray ionization (ESI) temperature, 550 • C; MS voltage, 5500 V in positive mode and −4500 V in negative mode; curtain gas (CUR) pressure, 25 psi; and collision-activated dissociation (CAD) parameter set to high. For the triple quadrupole (QQQ), each ion pair was scanned and detected according to its declustering potential (DP) and collision energy (CE) [28]. Multivariate statistical analysis was used to establish a reliable mathematical model summarizing the metabolic spectrum of the research object [29]. The data were scaled by unit variance and unsupervised principal component analysis (PCA) was performed using the prcomp function in R (www.r-project.org/, accessed on 23 October 2021). An orthogonal partial least squares-discriminant analysis (OPLS-DA) model was used to analyze the metabolome data, plot score, and permutation charts for each group, and reveal the differences among groups [30]. Variable influence on projection (VIP) values were extracted from the OPLS-DA results, which included generated score and permutation graphs. Tests were run on 200 permutations to avoid overfitting. Significantly differentially expressed metabolites were identified in each group. VIP ≥ 1, Foldchange ≥ 2 and ≤ 0.5 were used to screen differentially expressed metabolites for further analysis. The metabolites identified were annotated through the Kyoto Encyclopedia of Genes and Genomes (KEGG) (http://www.kegg.jp/kegg/compound/, accessed on 4 July 2021) [31] compound database, and then mapped to the KEGG pathway database (http://www.kegg.jp/kegg/pathway.htm, accessed on 4 July 2021).

Transcriptome Sequencing and Data Analysis
The experimental process of transcriptome sequencing includes RNA extraction, RNA detection, library construction, and computer sequencing. Sequencing and analysis were completed by Wuhan Metware Biotechnology Co., Ltd. (Wuhan, China. https://www. metware.cn, accessed on 18 April 2021). Total RNA was extracted from seeds of three independent plants at the same growth stage, selected among the four cultivars (red, white, yellow, and black quinoa cultivars), and the RNA was analyzed by agarose gel electrophoresis in RNA detection. Following the RNA quality inspection, for completeness and to determine whether there was DNA contamination, the library construction kit (NEBNext mRNA Sample Prep Reagent Set for Illumina, San Diego, CA, USA) was used to construct a sequencing library. For the library that met the requirements, the RNA concentration was measured with high accuracy using a qubit 2.0 fluorometer (Thermo Fisher Scientific, Waltham, MA, USA), and Agilent 2100 (Agilent Technologies, Santa Clara, CA, USA) was employed to detect the insert size of the library. Finally, 12 libraries representing three repeats and four grain samples of different color lines were constructed, and the transcriptome of the library on the illuminahiseq platform sequence was tested. After the gene expression level of each sample was obtained, the differentially expressed genes between samples were analyzed. DESeq2 (https://bioconductor.org/packages/ release/bioc/html/DESeq2.html, accessed on 10 November 2021) [32,33] was used for differential analysis to screen false discovery rate (FDR) < 0.05, | log2fold change | > = 1, and FDR < 0.05 was set as the threshold. After screening the differential genes according to the analysis purpose, a cluster heat map of different samples for functional annotation and enrichment analysis of differentially expressed genes, new gene analysis, variable splicing analysis, SNP, and indel analysis was produced.

Combined Transcriptome and Metabolome Analysis
The results of the differentially expressed metabolite (metabolome) analysis were combined with those of the transcriptome analysis. The genes showing altered transcriptomic, as well as metabolomic profiles, were mapped to the KEGG pathway chart, and histograms were plotted for them, showing the enrichment of pathways with both differential metabolites and differential genes. Correlation analyses were conducted on the differentially expressed genes and metabolites in each group. Pearson's correlation coefficients (PCC) of the genes and metabolites were calculated using the Cor program in R (www.r-project.org/, accessed on 3 November 2021). Genes and metabolites with PCC > 0.8 were selected to plot a network diagram representing the correlations among them. The overall correlations between indicators were reflected by the output of a canonical-correlation analysis (CCA) [34].

Conclusions
In this study, we applied a widely targeted metabolomics approach to analyze the primary and secondary metabolism of various quinoa cultivars. Thus, we found the metabolite profiles differed between white and black quinoa, and there are clear differences between red and yellow quinoa. The quinoa cultivars significantly differed, mainly in terms of their amino acid, tannin, alkaloid, and phenolic acid composition and content. Six common enrichment pathways, including phenylpropane biosynthesis, amino acid biosynthesis, and ABC transporter, were common to metabolites and genes. We identified key genes highly correlated with specific metabolites, and clarified the relationship between them. The results of this study provide a reliable basis for the cultivation of novel quinoa lines with superior yield, quality, and stress resistance characteristics. Moreover, this research empirically demonstrates the successful integration of metabolomics and transcriptomics for the comprehensive analysis of the metabolite profiles of quinoa and other crops, and the identification of the genes regulating these traits. However, the transgenic road of quinoa is still a problem that requires further attention.