The 179 lines from the RIL population developed from a cross between A2 (PTGMS A2) and Zhangzagu No.3 (Line 1484-5 × A2) was used in this study. The population was evaluated in natural field conditions in the experimental farm of Zhangjiakou Academy of Agricultural Science, Zhangjiakou, China (144°88’N, 40°77’E). A randomized block design was used during 2016-2018 cropping seasons. Each line was replicated three times (individual plants) and planted in a two-row plot of 1.5 m length with the spacing of 0.1 m between plants and 0.2 m between rows. Field management, including irrigation, fertilization, weeding and pest control, followed the standard agricultural practices in millet production. For each line, leaves from three plants were taken at the five-leaf stage and pooled for one biological replicate. Leaf samples for each line were selected for DNA or chemical extraction. All the samples were harvested at 10:00-12:00 on that day, placed in liquid N₂ immediately and stored at -70°C until vacuum freeze-drying.

The 179 accessions used in this study were characterized by whole genome re-sequencing. DNA was isolated from young leaves using the CTAB method and sequencing libraries with short inserts were constructed following manufacturer’s instructions (Illumina). The samples were sequenced on an Illumina HiSeq 4000 platform. To retain reads of high quality, reads with fewer than 5% N (missing) bases and with fewer than 50% of bases of base quality < 5 were deemed as cleaned reads. All other reads were discarded.

All the chemicals were of analytical reagent grade. Gradient-grade methanol, acetonitrile and acetic acid were purchased from Merck Company, Germany (http://www.merck-chemicals.com). The water was doubly deionised with a Milli-Q water purification system (Millipore, Bedford, MA). Authentic standards were purchased from ANPEL, Shanghai, China (www.anpel.com.cn/cnw), BioBioPha Co., Ltd. (http://www.biobiopha.com/), and Sigma-Aldrich, USA (http://www.sigmaaldrich.com). Standards stock solutions were prepared using water, methanol and/or dimethyl sulfoxide (DMSO) as the solvent and stored at -20°C. Combined standard solutions of chemicals were prepared just before use by mixing individual stock solutions and diluting these mixtures with 70% aqueous methanol.

The freeze-dried leaves were crushed using a mixer mill (MM 400, Retsch) with a zirconia beads for 1.5 min at 30 Hz. A 100 mg mass of powder was weighted and extracted overnight at 4°C with 1.0 ml of 70% aqueous methanol. Following centrifugation at 10, 000 g for 10 min, the extracts were filtered (SCAA-104, 0.22 μm pore size; ANPEL, Shanghai, China, http://www.anpel.com.cn/) before liquid chromatography-mass spectrometry (LC-MS) analysis.

The sample extracts were analyzed using an LC-ESI-MS/MS system (HPLC, Shim-pack UFLC SHIMADZU CBM30A system, www.shimadzu.com.cn/; MS, Applied Biosystems 6500 Q TRAP, www.appliedbiosystems.com.cn/). The analytical conditions were as follows, HPLC: column, Waters ACQUITY UPLC HSS T3 C18 (1.8 µm, 2.1 mm*100 mm); solvent system, water (0.04% acetic acid): acetonitrile (0.04% acetic acid); gradient program, 100:0 V/V at 0 min, 5:95V/V at 10.0 min, 5:95V/V at 11.0 min, 95:5 V/V at 11.1 min, 95:5 V/V at 15.0 min; flow rate, 0.35 ml/min; temperature, 40°C; and injection volume: 5 μl. The effluent was alternatively connected to an ESI-triple quadrupole-linear ion trap (Q TRAP)-MS.

LIT and triple quadrupole (QQQ) scans were acquired on a triple quadrupole-linear ion trap mass spectrometer (Q TRAP) using an API 6500 Q TRAP LC/MS/MS System, which was equipped with an ESI Turbo Ion-Spray interface operated in a positive ion mode and controlled by Analyst 1.6.3 software (AB Sciex). The ESI source operation parameters were as follows: ion source, turbo spray; source temperature 550°C; ion spray voltage (IS) 5,500 V; ion source gas I (GSI), gas II (GSII), curtain gas (CUR) were set at 55, 60, and 30.0 psi, respectively; and the collision gas (CAD) was high. Instrument tuning and mass calibration were performed with 10 and 100 μmol/L polypropylene glycol solutions in QQQ and LIT modes, respectively. The QQQ scans were acquired as MRM experiments with the collision gas (nitrogen) set to 5 psi. The DP and CE for individual MRM transitions were performed with further DP and CE optimization. A specific set of MRM transitions was monitored for each period according to the metabolites that were eluted within this period.

The metabolite data were log2-transformed for statistical analysis to improve normality. Broad-sense heritability (H²) was calculated using the following formula: H^{2 =} 1-10^{-2LOD/ n} , where n is the sample size (Arends et al., 2010). The values of the coefficient of variation (CV) were calculated for each metabolite and agronomic trait (three-year data are calculated separately) expressed as S/A, where S and A represent the standard deviation and the average of metabolite and agronomic trait in the population, respectively. Pearson’s correlation and the statistical significance between traits were estimated using programs housed in R (http://www.r-project.org/). Visualization correlation networks were constructed using Cytoscape 3.7.0 (Smoot et al., 2010).

A high-density genetic map was constructed for the RILs (Zhang et al., 2017). The QTL analysis of each trait was performed using the R package qtl version 1.46.2 (https://rqtl.org/), with a scanning step of 0.1 cM and PIN (probability in stepwise regression) of 0.01 (Li et al., 2007). The LOD threshold was set to 2.5 for both metabolites and agronomic traits. The confidence interval for each QTL was assigned as a 1.5-LOD drop of the peak. The additive effect and percentage of phenotypic variance associated with a QTL (contribution) were estimated using the same program. For metabolic QTL (mQTL), the QTL intervals of the same metabolite overlapped in two replicates, will be selected for research. If the phenotypic variance was greater than 15%, it was considered a major QTL (Salvi and Tuberosa, 2005).

The whole genome was divided into 3 Mb partitions, and the number of mQTL per partition was counted. Using 1,000 permutation tests, each mQTL was randomly assigned to a 3 Mb interval, and the number of mQTLs obtained in each interval was counted. The cut-off number of mQTLs per 3 Mb by chance alone was 14 in mature seeds with P < 0.05, respectively. A larger number was regarded as a mQTL hotspot (Gong et al., 2013).

The amino acid sequences of reported genes were obtained from NCBI according to their accession numbers (http://www.ncbi.nlm.nih.gov/). Candidate gene information in this study was obtained from the draft assembly of the millet genome (Bennetzen et al., 2012). The alignment of amino acid sequences was performed using ClustalW bundled in MEGA 5, and neighbor-joining trees were constructed using MEGA 5 software with all default parameters. The reliability of the reconstructed tree was evaluated using a bootstrap test with 1000 replicates.

A total of 63 agronomic traits were determined in 2016-2018 (Zhangjiakou) and the detailed information were shown in the Supplementary table 1 . Most phenotypes were measured using common methods such as plant height, thousand-grain weight and so on (Zhang et al., 2017; Fan et al., 2019). Other traits, such as anther color, which was observed at the time of floret dehiscence at anthesis, was scored on a scale of 1 point for white, 3 points for yellow and 5 points for brown; spike neck shape was to observe the bending degree and posture of stem node under spike, the morphology of these traits were evaluated and scored by experienced personnel. Briefly, genomic (2,202 bins of 33,579 SNPs integrated), metabolomic (3,452 metabolic signals) and multi-omics data (genomic and metabolomic data integration) were used to predict 63 agronomic traits using the BLUP method in R (rrBLUP) (Endelman, 2011) and LASSO methods in R (glmnet) (Hastie et al., 2017). The predictability was measured using a 10-fold cross-validation method. The 179 RILs were then randomly divided into 10 groups, 9 of which were used to construct the model. The remaining RILs were predicted. The predictive power (predictability) is defined as the Pearson’s correlation coefficient between the phenotypic observations and the predicted values (Shi et al., 2020).

Using widely targeted liquid chromatography-tandem mass spectrometry (LC-MS/MS)-based metabolic profiling method (Chen et al., 2013), we performed metabolic profiling with leaves at the five-leaf stage from 179 accessions derived from a cross between two elite foxtail millet varieties, A2 and Zhangzagu No.3. A total of 3,452 reproducible metabolite signals were detected, of which 381 metabolites were qualitatively analyzed through standard comparison and putatively annotated ( Figure 1A ). These annotated metabolites achieved a coverage of multiple important metabolic pathways for plants, including Flavonoids, Lipids, Phenolic acids, Amino acids and derivatives, Organic acids, Nucleotide and derivatives, Alkaloids, Anthocyanins, Lignans and Coumarins ( Figure 1B , Supplementary Table 2 ).

Pearson’s correlation coefficient of 381 annotated metabolites in 179 millets were calculated to explore the links between different metabolic pathways ( Figure 1C , Supplementary Table 3 ). We found that the number of metabolites with significant positive correlation (r > 0.5, p < 0.05, shown in red) were far more than those with significant negative correlation (r < -0.5, p < 0.05, shown in blue). We also found some closely related areas in correlation heatmap, especially for lipid metabolites, which indicated that these lipid metabolic pathways were relatively independent (Cluster 3, shown in blue box) of other types of metabolites. It is worth noting that flavonoids has formed two different clusters (Cluster 1, Cluster 2), and the two clusters showed significant negative correlation with each other. Cluster 1 (shown in red box) contained more flavonoid metabolites mainly derived from the Apigenin - Luteolin - Chrysoeriol - Selgin - Tricin- flavonoid lignin pathway ( Figure 1D ), while Cluster 2 (shown in green box) was mainly a group of flavonoids with the same modification that might have a relationship of metabolic substrate competition.

Coefficient of variation (CV) is often used to assess the extent of metabolic differences between populations. There were 76.64% annotated metabolites with a CV above 30%, ( Figure 1E , Supplementary Table 2 ), which suggested great variations of metabolites in different foxtail millet varieties. Moreover, the CV of primary metabolites were lower (t-test, p-value=1.63E-7) than that of secondary metabolites ( Figure 1F , Supplementary Table 4 ), indicating that secondary metabolites have greater variation than primary metabolites. Generally, secondary metabolites with greater variation in the population have important functions during plant growth and facilitate the mining of related genetic factors, such as anthocyanins and flavonoids, whose CV were 413.03% and 125.91%, respectively

The high-density genetic linkage map for “Zhangzagu No.03 × A2” RIL population used in this study was constructed as previously (Zhang et al., 2017). The linkage map consisted of 2,202 bin and 33,579 SNPs from all 9 chromosomes of foxtail millet. The genetic map spanned 1,934.6 cM of the foxtail millet genome, with average 0.96 cM per bin. Based on the high-density genetic linkage map, a total of 1,049 mQTLs (LOD > 2.5) from 992 metabolites were mapped, among which 114 annotated metabolites were mapped to 130 mQTLs ( Supplementary Table 5 ). According to the mapping results of mQTL, the proportion of mQTLs identified by secondary metabolites was higher than that in primary metabolites ( Supplementary Figures 1A, B ).

To have a more detailed understanding of mQTL results, we have also analyzed the distribution of mQTLs, a total of 11 mQTL hotspots on seven chromosomes were identified, indicating that some regulated genes of multiple metabolites may located in these regions ( Figure 2 ). It’s worth noting that there were 59 mQTLs from 55 annotated metabolites co-localized on the hotspot_5 ( Supplementary Table 5, 6 ). Most of these metabolites were involved in the phenylpropane and alkaloid-putrescine metabolic pathway, thus we believed that there may be regulatory genes of multiple metabolic pathways in hotspot_5. However, fewer mQTLs were detected on other chromosomes than expected, especially on chromosome 4 ( Supplementary Table 7 ).

Stable high-quality mQTLs can help us to discover potential candidate genes. Based on the structural characteristics of the metabolites and metabolic pathways have been reported, we screened 28 metabolite-related candidate genes from 14 mQTLs ( Table 1 ). The most significant mQTL for GZ2498 (Chrysoeriol-O-hexoside-malonyl-hexoside) was a segment on chromosome 3 in 49.35 Mb (241.3-251.4 cM, LOD = 3.18, PVE = 6.44%). As a candidate gene at this locus, Si021852m had high homology to the reported Arabidopsis gene AtDTX41 (Blastp E-value=0.0, Identity: 60.9%), which encodes a MATE efflux family protein involved in the pathway flavonoid biosynthesis (Zhao and Dixon, 2009). QTL of GZ2770 (Tricin O-sinapoylhexoside) was mapped to the 2.44-3.98 Mb on Chromosome 9 (28.3-45.8 cM, LOD=2.68 and PVE=6.44), Si006338m, Si008637m and Si006185m in the interval were homologous with a anthocyanidin 3-O-glucosyltransferase gene ZmBZ1 (Blastp E-value=0.0, Identity: 81.9%) in maize, a flavonoid 3-dioxygenase 3 gene OsF3’H-3 (Blastp E-value=4.4E-64, Identity: 38.0%) in rice and a quercetin 3-O-glucosltransferse gene AtUGT73B4 (Blastp E-value=0.0, Identity: 67.0%) in Arabidopsis, respectively (Roth et al., 1991; Lim et al., 2004; Kim et al., 2008; Byeon and Back, 2015). Besides, at the QTL of GZ2525 (apigenin-C-rutinoside) in 21.00-36.05 Mb on Chr 9 (112.6-132.7 cM, LOD=2.74 and PVE=4.02), we screened two genes Si035224m and Si035935m, they are highly homologous with the OsF3’H gene (Si035224m: Blastp E-value=0.0, Identity: 81.1%; Si035935m: Blastp E-value=0.0, Identity: 80.6%). In Rice, OsF3’H had been reported as flavonoid 3-hydroxylase which catalyzes the 3’-hydroxylation of the flavonoid B-ring to the 3’,4’-hydroxylated state (Lam et al., 2015). These all demonstrated the high quality and reliability of the mQTL of this study.

We noticed that multiple phenolic acids at the upstream of the phenylpropane metabolic pathway were mapped to hotspot_2 with high LOD and PVE on chromosome 1 ( Figure 3A ). Through the annotation of metabolite pathway and gene function, we had unearthed 6 candidate genes in this locus, including a Hydroxycinnamoyltransferase (HCT) Si016926m and 5 redundant Phenylalanine ammonia-lyase (PAL) ( Figure 3B ). Orthologs of these candidate genes had been reported to play an important role in the regulation of phenolic acid synthesis in the phenylpropane metabolic pathway ( Figure 3C ). Phylogenetic analysis showed that the PAL genes were clustering with the reported PAL from monocotyledonous cereal crops, like rice and triticum, rather than the PAL from dicotyledonous plants such as Arabidopsis (Zhu et al., 1995; Liao et al., 1996) ( Figure 3D ). Si016926m clustering with the HCT/HQT genes that had been reported to transfer hydroxycinnamate to shikimate or quinine, and is involved in the formation of phenolic acid (Chao et al., 2021) ( Figure 3E ).

In order to further dissect the function variation of candidate genes, we analyzed the SNPs within the candidate genes. A total of 30 SNPs in these 6 candidate genes were related to at least one metabolite in the pathway ( Supplementary Table 8 ). Among them, a non-synonymous mutation SNP (SNP1:3064324, C/T) in Si016926m caused the conversion between the basic amino acid lysine (Lys, K) and the neutral amino acid aspartyl (Asn, N). There were significant differences in metabolites content among different SNP haplotypes (p < 1E-03). As a upstream metabolite of phenylpropane metabolic pathway, GZ0329 (L-Phenylalanine) content in millet varieties with C is lower than those millet varieties with T or heterozygous in this site. On the contrary, the content of other downstream metabolites in millet with C had higher content than those millet varieties with T or heterozygous in the same SNP ( Figure 3F ). Moreover, a non-synonymous mutation SNP (SNP1:3176740, T/G) was found in the coding region of Si016467m, which can lead to the conversion between non-polar alanine (Ala, A) and polar threonine (Thr, T)The content of GZ0329 in the millet with T was significantly higher than that of the G, while other downstream substances such as GZ0548 (1-O-p-Coumaroyl quinic acid) and GZ2292 (3-O-Feruloyl quinic acid), showed the opposite trend ( Figure 3G ). Therefore, SNP1:3064324 in Si016926m and SNP1:3176740 in Si016467m did affect the synthesis of metabolites in phenolic acids pathway.

As the parents of the RIL population, Zhangzagu No.3 and A2 had significant differences not only in plant height (PH) and heading stage (HS), but also in number of grain (NG), length of ear neck (LEN/NL), gross weight (GW), ear length (EL), top second leaf length (TSL) and stem diameter (SD) through evaluation (t-test, p<0.05) ( Supplementary Figures 2A, B ). Further, we evaluated the important agronomic traits of millet throughout the growth period in the millet RIL population for three consecutive years (2016-2018) ( Supplementary Table 1 ). The descriptive statistics of each agronomic traits of the population were shown in Supplementary Table 9 . The agronomic traits of millet RIL populations were quite different, and the average phenotypic coefficient of variation reached 24.05%. By analyzing the correlation of agronomic traits, we tagged two cluster in the correlation heatmap ( Supplementary Figure 2C ). In the first cluster shown in blue box, there were significant positive correlations among 12 agronomic traits, such as ear shape (ES) and anther color (AC). Another cluster shown in green box containing five yield-related traits such as panicle weight (PW) and grain weight (GraW), also showed highly positive correlation.

We used the three-year phenotypic data to perform single-environment QTL analysis, and a total of 132 pQTLs (LOD > 2.5) were identified ( Figure 4A ; Supplementary Table 10 ). In the 2016 pQTL results, most of the traits were only mapped to one pQTL, of which the bristle color (BC) mapped a maximum of three pQTLs, while LEN/NL_2017 and PN_2018 can mapped five pQTLs, respectively. This confirmed that agronomic traits were closely linked with environmental factors and affected the pQTL results. Thus, stable pQTLs that had been repeatedly identified in more than one year for the same trait can more accurately reflect genetic characteristics of phenotypes. We identified 21 stable pQTLs across all chromosomes, and more than half of them explained a phenotypic variation of greater than 10% , especially qBc4-1 with the highest LOD and PVE ( Figure 4B ; Supplementary Table 11 ). Similar to the results of mQTL, the distribution of pQTL on the chromosome was not evenly ( Figure 4C ).

Based on the genome annotation and reported functional genes, we searched for potential candidate genes in the pQTL interval. A total of 39 candidate genes were mined from 25 pQTLs loci with 12 agronomic traits related to crop growth, development and yield formation ( Table 2 ). The stable pQTL qHs9-1 (LOD=4.99-11.45, PVE=10.33-20.85%) for HS located within 0.39-1.06 Mb on chromosome 9, where we found a transcription factor Si039184m which encodes a GHD protein ( Supplementary Figure 3A, B ). The homologous gene OsGHD7 of Si039184m (Blastp E-Value=8.19E-59, Identity: 49.8%) in rice had been reported to play a major role in regulating the flowering time under long day (LD) conditions by negatively regulating the expression of EHD1 and HD3A ( Supplementary Figure 3C ) (Hori et al., 2013). Another pQTL qHs9-2 (LOD=3.36, PVE=5.62%) located in the 57.82-58.44 Mb on chromosome 9 was also identified, in this loucs, we found a candidate gene Si034749m encoding EHD4 protein domain ( Supplementary Figures 3A, B ). The function of its homologous gene OsEHD4 (Blastp E-Value=5.86E-93, Identity: 40.6%) in rice had been proved to be related to flowering as a downstream gene of GDH7 (Gao et al., 2013). Phylogenetic analysis revealed that both Si039184m and Si034749m were grouped together with the reported GDH genes and EHD genes from monocotyledonous cereal crops, respectively, which indicated the consistent function of them ( Supplementary Figure 3D ). At the locus of qTne9-1 (57.07-57.82 Mb, LOD=2.74, PVE=4.41%) near qHs9-2, a yield-related gene Si035172m attracted attention. Its homologous gene OsTE in rice (Blastp E-Value=0.0, identity: 93.0%) encodes a Cdh1 protein belonging to a co-activator of APC/C. The APC/C-TE complex has special effect on regulating the lateral branch and tillering of the embryo, which is an important factor to determine the plant type and grain yield ( Supplementary Figures 3E, F, G ) (Netty and Yossi, 2017). The utilization of this segment will be beneficial to cultivate the planting resources of early-maturing and high-yielding millet.

The yield-related pQTLs were simultaneously mapped in the 33.43 Mb-37.15 Mb interval of chromosome 4, including qTgw4-1, qPn4-1, qPw4-1, qGraW4-1 ( Supplementary Figure 3H ). The co-location pQTLs of the similly traits was also beneficial for candidate gene mining, and the gene Si008391m with annotated function of cytokinin dehydrogenase was excavated in this locus ( Supplementary Figure 3I ). It was found that OsCKX2 gene (E-Value=6.09E-164, identity: 52.7%) is a homologous gene of Si008391m in rice, and OsCKX2 affected rice yield by regulating the content of cytokinin ( Supplementary Figure 3J ) (Li et al., 2013). This provided important evidence and resources for genetic breeding related to yield in millet.

In order to explore the relationship between phenotypic characters and metabolome, we focused on co-localized metabolites that were highly correlated with the phenotypic traits. We built a metabolite-agronomic trait association network based on the correlation between 381 annotated metabolites and 63 phenotypes ( Figure 5A ). A total of 498 significant correlations were screened (|r|>0.3, p<0.05), which suggested these metabolites may be involved in the formation of agronomic traits ( Supplementary Table 12 ). Among them, most of the metabolites were in the phenylpropane metabolic pathway. Besides, there were 104 mQTLs of annotated metabolites co-located with pQTLs ( Figure 5B ; Supplementary Table 13 ). A total of 48 metabolites were significantly correlated with their corresponding phenotypes, indicating that these genetic loci may affect the phenotype by regulating the content of metabolites. Notably, most of the annotated metabolites with co-localization were in the phenylpropane metabolic pathway, such as GZ0852 (apigenin 7-rutinoside), which was co-located with SD (2018), and the Pearson’s correlation coefficient was 0.98 ( Supplementary Table 12 ).

Since the detection of metabolites and the investigation of agronomic traits were carried out in different periods, the high correlation between metabolites and agronomic traits suggested the presence of a genetic relationship. We found that multiple metabolites and agronomic traits were mapped on the same locus 41.33-44.15 Mb of chromosome 5 with high correlation ( Figures 6A, B ). Among them, the metabolites of the flavonoid-lignin synthesis pathway were positively correlated with architecture-related traits, and negatively correlated with yield-related traits ( Figures 6C, D ). Through functional annotation analysis, we found that the Si001573m in this interval was the homologous gene of the Sd1 known as the green revolution gene in Rice (Blastp E-Value=0.0, identity: 82.8%), which can affect plant height by regulating GA synthesis (Sasaki et al., 2002; Spielmeyer et al., 2002; Peng et al., 2021). Phylogenetic analysis showed that Si001573m grouped with Gibberellin 20 oxidase of other crops, including Sd1 (OsGA20ox2) in rice ( Figure 6E ). Interstingly, we screened four candidate genes of flavonoid-lignin synthsis pathway nearly the Si001573m ( Figure 6B ). Among those, Si000845m was annotated as an anthocyanin regulatory Lc protein, its homologous gene in maize has been confirmed to regulate the glycosylation of flavonoids (Blastp E-Value=0.0, identity: 71.0%) (Ludwig et al., 1989; Nesi et al., 2000). Phylogenetic analysis revealed that Si000845m was also highly homologous to Lc transcription factors in other crops ( Figure 6F ). In addition, Si001301m and Si001332m encoded glycosyltransferase proteins. Their homologous genes have been reported to play a role in glycosylation modification of flavonoids in strawberry (Blastp E-Value=5.2E-129, identity: 43.8% and Blastp E-Value=5.8E-132, identity: 45.3%, respectively). The catalytic function of Si001301m and Si001332m for rhamnosyl modification of flavonoids was further determined by phylogenetic analysis ( Figure 6G ) (Lunkenbein et al., 2006). Si001026m has been annotated as a flavonoid hydroxyl ligase, and its homologous genes can be used for hydroxylation modification of flavonoid in licorice (Blastp E-Value=4.98E-137, identity: 57.6%). Phylogenetic analysis showed that Si001026m was clustered with F2’H reported in other species ( Figure 6E ) (Akashi et al., 1998).

There is still a lack of reliable reports or evidences to verify the relationship between metabolites and phenotypes, despite the presence of co-localizations and strong correlation in millet. There were five genes located close together within an ~1.81 Mb region (Chr5: 41.79-43.60 Mb). We speculated that there might be a close linkage between them, causing the non-separation of the same genotype in the hybrid offspring. Therefore, plant height and flavonoid-lignin were highly correlated. In order to confirm this hypothesis, we focus on SNP5: 43169135 near the gene Si001573m that formed two different genotypes of C and T in the population. Further analysis revealed that there were significant differences in PH value of millet under the two genotypes (p=3.33E-19), and the content of flavonoid- lignin pathway metabolites such as GZ2761 (3’,4’,5’-Tricetin O-rutinoside) and GZ3060 (Tricin 4’-O-syringic acid) also showed the same trend ( Figure 6H ), which provided important resources for the development of molecular markers related to functional metabolites and agronomic traits at the same time.

In this study, we used widely-targeted LC-MS/MS-based metabolic profiling method to detect 3,452 metabolite signals in millet seedling leaves, and elucidated 381 metabolites through standard comparison and chromatography-mass spectrometry data, including primary metabolites and secondary metabolites. Additionally, the coefficient of variation (CV) of primary metabolites was significantly lower than that of the secondary metabolites, which was in line with the results in other cereal crops such as rice, wheat and barley. Interestingly, flavonoids in millet and barley, especially anthocyanins, have the largest CV. However, in wheat, the CV of phenolamide alkaloids was the largest, and the CV of flavonoids and anthocyanins was relatively small (Chen et al., 2013; Chen et al., 2020a; Shi et al., 2020; Zeng et al., 2020).

The correlation between metabolites not only reflect the relationship among metabolites in synthesis pathways and genetic regulation, but also help to predict the structure of unknown metabolites. We observed a significant association between metabolites of the same class, such as flavonoids, phenolic acids, lipids, amino acids and nucleotide. Metabolites in adjacent pathways from different classes still have significant associations, such as flavonoids - anthocyanins, and flavonoids - phenolic acids, which were consistent with the results in wheat and rice (Gong et al., 2013; Shi et al., 2020). Apart from this, some metabolites were not closely connected in metabolic pathways were also highly correlated. For example, the Pearson’s correlation coefficient of GZ2056 (N-(4’-O-glycosyl)-p-coumaroyl agmatine) and GZ1939 (4-Pyridoxic acid O-hexoside) reached 0.66, possibly due to the similar glycosyl structure. Co-localization analysis showed that they were mapped on the same locus of 24.24-28.42 Mb on chromosome 1, indicating that their glycosylation modification may be regulated by the same site ( Supplementary Tables 3, 5 ).

In addition to the annotated metabolites, we had also detected a series of unknown signals. The correlation between the unknown metabolites and the annotated metabolites, as well as the co-localization analysis of mQTLs can provide an important resource for future efforts in the identification of unknown metabolites and pathways (Shin et al., 2014; Do et al., 2015). For instance, the unknown metabolite GZ2531 had a high correlation with multiple flavonoid metabolites in the Eriodictyol-Apigenin pathway with the average correlation of 0.82 ( Supplementary Table 14 ). They were mapped on the same loucs at 6.30-8.89 Mb on chromosome 9. Moreover, GZ2531 and GZ2536 (Chrysoeriol 8-C-hexoside) had the same Q1, thus they were likely to be isomers of each other ( Supplementary Table 5 ).

We used the high-density genetic linkage map of millet to identify 1,049 mQTLs of 922 metabolites. Consistent with the mQTL research on rice, wheat, Arabidopsis, maize and apple, we found that mQTLs were mainly distributed in the form of hot spots on the genome (Khan et al., 2012; Gong et al., 2013; Wen et al., 2015; Knoch et al., 2017; Shi et al., 2020). For example, 1,005 mQTLs identified using metabolites in mature wheat grains were distributed on 68 hotspots, such as chromosomes 1B, 4B and 7A (Shi et al., 2020); mQTLs for rice grains and flag leaves were also distributed on 4 hotspots and 2 hotspots, respectively (Gong et al., 2013). In particular, there was a major hotspot in the 40-44Mb interval on chromosome 5 in millet where most of the mQTLs were located, which was consistent with the research in apple (Khan et al., 2012). Different from other plants, the mQTLs of the same class metabolites in millet were not evenly distributed on chromosomes. For example, the mQTLs of flavonoids were concentrated on two hotspots of chromosome 5 and 9, while the nucleotide mQTLs were only concentrated on the hotspots of the chromosome 5.

We compared the present results with previous mGWAS study based on the genetic mechanism of metabolites in natural variation populations of millet, and found that there were differences in the genetic sites identified for the same metabolite (Wei et al., 2021). There were 49 metabolites with significant association interval, including 59 mQTLs and 350 significant SNPs. From 59 mQTLs, we screened 34 metabolite-related 36 important mQTLs (LOD >5, PVE >10%), and only 12 metabolites could detect similar genetic loci by both study, including 9 unknown metabolites (GZ0797, GZ0814, GZ0894, GZ1055, GZ1201, GZ1272, GZ2255, GZ2256, GZ2300) and 3 annotated metabolites (GZ0755 pentoside Caffeate, GZ2305 (+)-Gallocatechin -hexoside, and GZ2800 Tricin O-vanilloylhexoside). Unexpectedly, we found that different metabolites of the same pathway may be detected at the same genetic locus in different populations. Therefore, mQTL and mGWAS can not only verify the accuracy of metabolite genetic locus mining, but also help improve the metabolite pathway.

A total of 132 pQTLs were identified by single-environment QTL analysis using the phenotypic data of three years. By analyzing the distribution of all the pQTLs, we found 12 hotspots (all of which contained at least 3 pQTLs) on six chromosomes including 81 pQTLs. Among them, hotspots_9 contained the most pQTLs and covered with 15 traits, was located at 42-45 Mb on chromosome 5. In hostspots_9, the range of LOD was 3.63 to 28.24, and the range of PVE was 0.39 to 53.20% ( Supplementary Table 15, 16 ). Based on the physical coordinates of the QTL confidence intervals, we compared the present results with previous reports. Four pQTLs (qPw5-2, qPh1-1, qPh5-1, qTgw5-1) were overlapped with the genomic regions of qpw5, qph1, qph5 and qtgw5 that were isolated from 439 RIL populations in foxtail millet (Zhang et al., 2017). The physical position of qPw6-1 was overlapped with that of qpw6.2 for PW detected in a backbone line Ai88 × Liaogu1 F2 population (Zhi et al., 2021). Above all, the QTL analysis in this study was reliable to further explore the genetic relationship between metabolites and phenotypes.

Metabolites are the closest phenotypic link in the process from heredity to phenotype (Luo, 2015). Metabolome has striking advantages in analyzing less observable phenotypes, which can clearly reflect small changes in phenotypes from metabolite content levels. Through correlation and co-localization analysis, we identified strong associations between metabolites and agronomic traits in millet. For example, PC and two anthocyanin metabolites GZ2284 (Cyanidin 3-O-malonylhexoside) and GZ2423 (Pelargonidin 3O-malonyl-malonylhexoside) showed a significant positive correlation. GZ2284 and GZ2423 could effectively regulate plant color, which was consistent with the research in Arabidopsis and apple (Jiang et al., 2019; Kim et al., 2019). Studies in rice, sorghum and wheat have shown that the downstream metabolic pathways of phenylpropane, tryptophan and tyrosine of shikimic acid can regulate the yield of cereal crops, showing a significant positive correlation with yield (Chen et al., 2020b). However, we found the yield-traits of millet were negatively correlated to GZ0329 (L-Phenylalanine), GZ0200 (L-(-)-Tyrosine) and flavonoid metabolites (based on Tricin-related metabolites), and positively correlated to unsaturated fatty acids, succinic acid. This may because of the specificity of metabolites accumulation pattern in different species or tissues.

The genetic basis of metabolites and phenotypes is mainly divided into three types. First, functional genes that regulate metabolites are distributed near phenotype-related genes, there may be close linkages between these genes. Selection for agronomic traits during the breeding process will also select genes that regulate metabolites through the free-riding effect (Zhu et al., 2018). Secondly, transcription factors that regulate both phenotypes and metabolic pathways can simultaneously affect specific agronomic traits and metabolites (Pillet et al., 2015). Finally, genes can influence the phenotypes by regulating the content of metabolites that affect the phenotypes (Jiang et al., 2019). At present, the identification technology of trait genetic locus has been gradually mature, but the process of mining the candidate genes remains a huge challenge. Many studies only identified QTL or SNP markers and failed to dig out important candidate genes within the locus. According to further genetic mechanism research, we found that metabolites and phenotypes with strong relevance were located in the same locus, which was consistent with the research in wheat and rice (Gong et al., 2013; Shi et al., 2020). In this study, combined analysis of mQTL and pQTL can effectively help narrow the candidate interval and validate key phenotype-related genes.

With the progress of various omics technologies, multi-omics data have been used to predict complex agronomic traits. Such application enriches the methods of molecular marker-assisted breeding and brings breakthroughs for the genetic improvement of crops (Xu et al., 2014; Wang et al., 2017; Xu et al., 2017). Due to the strong association between seedling metabolites and important agronomic traits throughout the growth period, we attempted to predict agronomic traits using multi-omics data in this study. Based on least absolute shrinkage and selection operator (LASSO) and BLUP models, genomic (2,202 bins of 33,579 SNPs integrated), metabolomic (3,452 metabolic signals) and multi-omics data (genomic and metabolomic data integration) were used to predict the 63 agronomic traits, respectively. Due to the sparse solution characteristics of the LASSO model, the phenotypes were not completely predicted. The average predictability of different data were 0.60, 0.57, and 0.60, respectively. The BLUP model can complete the prediction of all phenotypes with an average predictability of 0.68, 0.66 and 0.79. The predictability of all phenotypes using the BLUP model and multi-omics data has reached more than 0.5, indicating that this method can be used to predict the phenotype of millet and achieve satisfactory results ( Supplementary Table 17 ). Among the phenotypes, PN_2017 showed the highest predictability, with an average of 0.85 under all methods ( Supplementary Table 18 ).

Comparing the prediction results of different omics data by two models, the prediction using BLUP model and multi-omics data was significantly better than other methods. ( Figure 7A and Supplementary Table 17 ). These results were consistent with the reports in rice and maize (Riedelsheimer et al., 2012; Xu et al., 2016). Based on the correlation analysis of the relationship between the predictability and heritability, we found that the average correlation between predictability of genome participation and heritability was 0.71, 0.69, 0.58, 0.71, respectively, while the correlation between the predictability of phenotype by the metabolome and heritability was smaller, showed 0.36, 0.44, respectively ( Supplementary Table 17 ). It showed that the use of genomes for phenotype prediction comparisons relies on the heritability, while metabolomes have advantages in predicting low-heritability phenotypes. For example, the heritability of TPN_2016 was only 0.09, and its predictability based on the BLUP model using genomic and metabolomic data reached 0.48 and 0.79, respectively ( Figure 7B ). Similar results were obtained in rice, metabolomic prediction for YIELD with low heritability was almost twice as efficient as genomic prediction (Xu et al., 2016), and it may be one of the reasons for the better predictability of multi-omics data.