1. Introduction
Foxtail millet (
Setaria italica L.) is an ancient crop originated in China, renowned for its high drought tolerance [
1]. It has emerged as a prominent cereal in regions with limited water resources, including China, India, and certain African countries [
2]. Foxtail millet grains are abundant in essential nutrients, such as organic acids [
3], vitamin E, various amino acids, and high quality protein [
4]. Additionally, it serves as a valuable source of trace elements like zinc and iron [
5]. Its bran also contains substantial amounts of linoleic and linolenic acids [
6,
7], making it an excellent crude fiber source that aids intestinal digestion and promotes digestive health [
8]. These nutritional attributes have garnered increasing global attention towards this grain. However, its yield per unit area significantly lags behind staple crops like maize, wheat, and rice, and the planting area and total yield have been a declining trend [
9]. Therefore, enhancing the yield potential of foxtail millet becomes imperative.
Cereal crops grain yield is determined by three component traits, the number of panicles, the number of grains per panicle and thousand kernel weight (TKW) [
10]. Kernel weight, the most crucial factor determining yield, is largely influenced by the grain filling duration (GFD) and grain filling rate (GFR) [
11]. Previous studies have shown that GFD is greatly affected by temperature [
12,
13], nitrogenous fertilizer [
14], especially during the grain filling period under water stress [
15]. However, the stability of GFR and its limited susceptibility to environmental factors [
16,
17] suggest a potential genetic determination for this parameter [
18]. More importantly, the GFR plays a crucial role in enhancing grain productivity by facilitating endosperm development and nutrient accumulation [
19]. The performance of cultivars with high filling rates is generally superior under adverse conditions such as abiotic and biotic stresses, as they produce seeds that are fuller and exhibit lower effects compared to cultivars with low filling rates [
20]. Thus, the attainment of crops with high GFR can be accomplished through genetic breeding, and the strategic selection and breeding of genotypes exhibiting high GFR will prove to be a successful approach in ensuring optimal grain yield. Furthermore, exploring and identifying regulatory genes and key loci associated with GFR, as well as analyzing and definiting the regulatory pathways and molecular mechanisms governing GFR, we can facilitate precise breeding of rapid filling characteristics, thereby enhancing both grain yield and quality.
The identification of several pivotal genes involved in GFR, along with the elucidation of their molecular pathways in major crops, has significantly contributed to our understanding of this intricate process. In rice, the enzymes associated with assimilate synthesis and transport in rice, including those involved in starch synthesis and sucrose metabolism, exert a significant impact on GFR.
Grain incomplete filling2 (GIF2), encoding an ADP-Glc pyrophosphorylase large subunit, regulats grain filling and starch biosynthesis during caryopsis development. Inhibition of GIF2 expression level leads to a decrease in GFR and yield [
21]. The sucrose transport-related genes,
Sucrose Transporter1 (
OsSUT1),
OsSUT2,
OsSUT3,
OsSUT4, have been reported to exert significant influence on grain filling due to their involvement in the transportation of sucrose.
OsNF-YB1 acts as an upstream regulator of
OsSUT3 and
OsSUT4 expression, specifically expressed in the aleurone layer. [
22,
23]. Moreover, the SWEET genes, such as
ZmSWEET4c in maize and
OsSWEET4 in rice, regulate the seed filling process by mediating transepithelial hexose transport [
24]. In maize,
Opaque2 activates the expression of sucrose synthase genes
SUS1 and
SUS2 and promotes grain filling [
25]. Furthermore, several studies have indicated that the members of the cell wall invertase (CIN) gene family impact the GFR. In rice, the cell wall convertase genes
OsGIF1/
OsCIN2 have been identified as domesticating selective genes involved in carbon allocation during early growth stage. These genes regulate transport and unloading of sucrose and affect GFR [
26]. And the change in plant hormone levels also strongly influences GFR. Ethylene (ETH) and abscisic acid (ABA) regulate the activity of enzymes involved in starch synthesis, thereby further inhibiting carbohydrates biosynthesis in spikelets, reducing GFR and ultimately affecting yield [
27,
28]. And the low concentrations of auxin (IAA) and cytokinin (CK) in rice lead to a decrease in GFR, resulting in poor grain filling, ultimately leading to reduced grain weight and yield [
29,
30].
DEP1/
qPE9-1 is the rice G protein γ subunit, which regulats GFR by increasing the auxin and CK content of grains [
31]. As well as, microRNAs (miRNAs) play a crucial role during grain filling, and the suppression of miR1432 expression enhances GFR, leading to a significantly increase grain yield in rice [
19].
GRAIN-FILLING RATE1 (GFR1) encodes a membrane-localized protein that enhances the expression of Rubisco genes in the calvin cycle, thereby promoting the synthesis of sucrose and increasing the grain filling rate and yield in rice [
32]. And
QGfr.sicau-7D.1 is a newly discovered QTL for GFR in wheat [
33].
However, the study on molecular mechanisms and regulatory genes related to GFR in foxtail millet is rarely reported compared to major cereal crops. The whole genome sequencing of foxtail millet has facilitated further research in the field of molecular biology and functional genomics for this crop [
34,
35]. High-throughput sequencing technology and multi-omics analysis have made it possible to identify genes associated with grain filling. By comparing the difference in transcriptome [
9,
36,
37], metabolome [
38] and genome-wide DNA methylation in foxtail millet [
39] at different grain filling stages, numerous potential regulators in grain filling and the accumulation of grain metabolites have been identified [
35,
39]. Furthermore, one papers provided preliminary insights into the differences in gene expression and molecular mechanisms of grain filling among foxtail millet cultivars with different panicle types in the North China summer-sowing region [
40]. Nevertheless, those studies solely focus on elucidating the mechanism underlying foxtail millet grain filling process, without delving into the mechanism for differences in the grain filling process of two genotypes (high GFR and low GFR). In particular, the molecular pathways and associated genes governing GFR have not been identified.
In this study, we conducted an analysis on the variations in grain weight and GFR across six different developmental stages, utilizing the grains of high-GFR and low-GFR foxtail millet cultivars. Then, the stages with the most significant disparity in grain filling rate between the two cultivars was chosen for transcriptomic and metabolomic investigations. Our comprehensive analysis of the transcriptome and metabolome revealed potentially crucial metabolites and transcriptional regulators, shedding light on the molecular regulatory mechanism underlying grain-filling rate. Furthermore, this study provides a theoretical foundation for achieving high-yield and high-quality molecular design breeding of foxtail millet.
4. Discussion
Foxtail millet grain filling progress is crucial for seed setting, grain weight, quality, and yield [
39]. This process is highly sophisticated, involving the accumulation of carbohydrates and other nutrients by the developing endosperm to synthesize starch and enhance grain development [
52]. The progress of foxtail millet grain filling is manifested in four stages: the ovule stage before pollination (ovule stage, 1-5 DAF), the milk stage after pollination (milk stage, 6-11 DAF), nutrient storage stage after pollination (dough stage, 12-22 DAF), and the grain maturity stage (maturity stage, 23-30 DAF) [
9]. The interaction of phytohormones, the activity of starch biosynthesis enzymes, the levels of polyamines, and the synthesis and translocation of assimilates all play essential roles in grain filling [
24,
53,
54]. Various factors, such as photosynthetic products, transportation of stored substances, transport tissues, physiological activities of the grains, regulation of plant hormones, and environmental factors, can restrict this progress [
55]. Enhancing GFR, particularly improving the rate at which grains are filled with nutrients and water, has been demonstrated as an effective strategy to enhance crop yield and quality [
56]. The grain filling rate is an important agronomic trait primarily influenced by genetic factors, while research on the molecular mechanisms of grain development in cereal crops has only been carried out in recent years. Some genes and quantitative trait loci (QTLs) associated with the grain size and GFR have been identified and cloned in rice, maize or wheat [
32,
33,
57]. However, the molecular mechanisms related to GFR in foxtail millet have not yet been reported.
CS13 and CN47 are well-suited for cultivation in the mid-late mature spring foxtail millet region of northwest China, exhibiting a growth period of 122 days and a heading date ranging from 48 to 50 days [
58,
59]. The two cultivars had similar agronomic traits, except for the thousand-grain weight. Nevertheless, we found that CN47 had a significantly higher grain filling rate at 14 DAF and 21 DAF stages compared to CS13. Therefore, CN47 exhibited high grain filling rate, whereas CS13 displayed low grain filling rate. The significant disparity in grain filling rate between CN47 and CS13 made it the optimal material for investigating the grain filling rate of spring foxtail millet in the middle and late maturity regions of northwest China. So we sampled and sequenced the grains of foxtail millet in both periods. Therefore, this study has provided experimental materials for exploring genes and pathways related to the regulation of grain filling rate.
The utilization of RAN-seq is indispensable in the investigation of plant development processes, facilitating our comprehension of gene expression levels across various stages or cultivars [
60]. The metabolites, serving as the ultimate outcomes of gene expression, are intricately associated with the development of plant traits. It is a common approach to identify pivotal genes and elucidate the mechanism underlying phenotypic variations through the integrated analysis of transcriptome and metabolome. The transcriptome analysis conducted in foxtail millet at five distinct grain filling stages initially revealed the involvement of pathways related to starch biosynthesis, cell-wall invertases, hormone signal transduction, and polyamine metabolism in the process of grain filling [
36]. Subsequently, the transcript dynamics involved in foxtail millet grain development at four stages were investigated using RNA-seq. The results revealed the intricate mechanisms underlying plant hormone signaling, starch and sugar metabolism, carotenoid metabolism, flavonoid biosynthesis, and folate metabolism pathways during grain development in foxtail millet [
9]. Moreover, a comparative transcriptome analysis also revealed significant variations in the expression levels of genes related to the ABC transporter, photosynthesis, and the photosynthesis-antenna protein pathway. Besides, a large-scale metabolome study indicated dynamic changes in flavonoid, glutathione, linoleic acid, starch , sucrose, valine, leucine and isoleucine during foxtail millet grain filling [
38]. In this study, the transcriptome analysis conducted at 14 DAF and 21 DAF stages revealed inconsistencies with previous reports in terms of significantly enriched pathways. However, metabolome analysis revealed that pathways such as flavonoid biosynthesis, ABC transporters, photosynthesis, and plant hormone signal transduction were consistent with previous research findings. The findings indicated that a limited number of DEGs were implicated in pathways associated with the process of grain filling. The DEGs and DAMs identified between CS13 and CN47 at 14 DAF and 21 DAF stages were found to be associated with pathways such as flavonoid biosynthesis, plant hormone signal transduction, tyrosine metabolism, ABC transporters, beta-Alanine metabolism, as well as stilbenoid, diarylheptanoid and gingerol biosynthesis. The results obtained above demonstrate that flavonoid biosynthesis, ABC transporters, and plant hormone signal transduction pathways play a crucial role in influencing the grain filling process of foxtail millet. Importantly, these pathways have also been identified as potentially significant factors contributing to the disparity in GFR between CS13 and CN47.
Flavonoids are the predominant yellow pigments found in plants and are believed to confer nutritional and pharmacological benefits to plants [
61]. They play essential roles in enhancing the overall quality of plants, encompassing aspects such as color and flavor [
62]. The study revealed that foxtail millets were abundant in flavonoids, such as vitexin, naringenin, and naringenin chalcone. It was widely acknowledged that the flavonoid metabolic pathway significantly influenced the development of a yellow hue in de-husked grains. In essence, the flavonoid metabolic pathway played a crucial role in grain filling and GFR. In the study, the integration of the DEGs and DAMs between CS13 and CN47 at 14 DAF and 21 DAF stages identified co-enrichment of flavonoid biosynthesis. Specifically, the pathway only contained one gene,
SETIT_015292mg, which encoded a methyltransferase involved in the post-modification of flavonoids. The gene clearly encoded flavonoid O-methyltransferase (FOMT), which serves as a key modifying enzyme in the flavonoid metabolic pathway by catalyzing the synthesis of O-methylated derivatives of flavonoids [
63]. Furthermore,
SETIT_015292mg also played a crucial role in the stilbenoid, diarylheptanoid and gingerol biosynthesis pathway (
Table S8). Interestingly, the expression of
SETIT_015292mg was observed exclusively in CN47, but not in CS13, which could potentially serve as a significant contributing factor to the high GFR in CN 47 (
Figure S2).
The ABC transporters are essential for normal plant development and play a pivotal role in diverse biological processes, encompassing auxin transport and seed development in plants [
64]. In cereal crops, numerous genes related to ABC transporters had been found and identified to be linked with grain development and grain filling.
TaABCC3 promoted grain formation in wheat, while its suppression resulted in a reduction in the number of grains [
65]. And wheat
ABCC13 was reported functionally essential for grain development [
66]. Similarly, strong expression of
OsABCI15 and
OsABCI16 was exhibited in the seed, including their involvement in seed development in rice. Overexpression of either gene significantly enhanced the grain yield, while CRISPR/Cas9-mediated loss-of-function mutations in these genes resulted in incomplete filling of developing seeds [
67]. The expression levels of several ABC transporter-related DEGs, including ABCB and ABCC genes, were found to be significantly upregulated during the late stages of grain filling in foxtail millets [
47]. In contrast, our study found that two ABCB1 genes in ABC transporters pathway exerted influences on grain filling at different stages by regulating auxin transport. These genes,
SETIT_032403mg and
SETIT_000076mg, exhibited different expression patterns between CS13 and CN47 at 14 DAF and 21 DAF stages, suggesting diverse roles of ABC transporters in GFR.
Phytohormones play an extremely crucial role in kernel development and grain filling in crops. So far, numerous genes related to plant hormones have been identified. Auxin can induce the expression of growth-regulating genes to regulate the panicle and spikelet development in rice [
68]. Furthermore, a study had reported that auxin-response factors potentially influence early foxtail millet grain development in foxtail millet [
47]. Additionally, some studies showed that CTK and ETH were involved in regulating grain filling [
69,
70]. In our study, we observed that many DEGs were related to plant hormone biosynthesis and/or signaling pathways, these genes included YUCCA, ABCB, AUXIN1, AUX/IAA, SAUR, IPTs, AHP, ARR, EIN3. Notably, Indole-3-acetic acid (IAA) was detected as a DAM, and the majority of genes associated with auxin were found in the plant hormone signal transduction pathway. Genes associated with the auxin signaling pathway may serve as valuable resources for molecular breeding of foxtail millet. Therefore, our primary focus lies on elucidating the mechanisms underlying auxin-mediated grain filling and exploring strategies to regulate the growth factor receptor (GFR).
Auxin, a phytohormone, is known to be crucial for endosperm development and grain filling [
57]. YUCCA flavin monooxygenases are involved in auxin biosynthesis, and the rice auxin biosynthesis gene (
OsYUC11) is essential for grain filling in rice. Mutations in
OsYUC11 have been shown to hinder grain filling and the accumulation of storage proteins [
71,
72]. Blocking auxin transport led to abnormalities of seed [
73]. The intercellular transport of auxin, facilitated by PIN [
74], ABCB [
75], and AUXIN [
76], plays a crucial role in establishing and maintaining optimal concentrations of auxin. It reported that the ABCB proteins were involved in auxin transport, with ABCB1, ABCB4, and ABCB19 being responsible for directional transport in rice [
77]. The auxin signal transduction pathway, encompassing TIR1/AFB co-receptors, Aux/IAA repressors, and ARF transcription factors, is recognized for its regulatory role in grain yield. Nevertheless, its involvement in the grain filling process remains ambiguous [
78,
79]. SAUR genes, which are early auxin-responsive genes, play a role in plant growth, particularly in cell elongation [
80]. In this study, we identified several vital genes related to auxin biosynthesis, transport and signal transduction. Contrary to findings in other plant species, the expression of a YUCCA gene was found to be higher in CS13 with low GFR compared to CN47 with high GFR at 14 DAF. Interestingly, the variation in auxin levels between the two cultivars was observed at 21 DAF instead of at 14 DAF. It is hypothesized that YUCCA may usexert a negative regulatory effect on auxin synthesis, with its impact being delayed. Furthermore, several genes related to auxin transport pathway have been identified, belonging to the ABCB gene family and AUXIN gene family. In rice, the low expression of
OsLAX1/
OsAUX1 was found to increased root length [
81]. Similarly the low expression level of four genes in CN47 may also lead to an increase in grain length by promoting grain filling rate. The expression levels of
SETIT_007373mg and
SETIT_028123mg as auxin response gene exhibited contrasting trends. The expression level of AUX/IAA gene
SETIT_007373mg was found to be higher in CS13, whereas a decrease in expression was observed for the SAUR gene
SETIT_028123mg. The expression profiles of different auxin response genes may vary in their responsiveness to auxin signaling pathways. Based on these findings, it is reasonable to speculate that these genes may play distinct regulatory roles at different stages; further investigation is required to elucidate their involvement in the grain filling process.
Figure 1.
Phenotypic of grain filling rate related in foxtail millets. (a) and (b) Morphological changes in foxtail millet grain of CS13 (a) and CN47 (b) at six stages. Bars=500 μm. (c) Thousand millets weight (de-husked grains) changes of CS13 and CN47 at six stages. (d) Grain-filling rate changes of CS13 and CN47 at six stages. (e) Comparison of the agronomic traits between CS13 and CN47. The data are means ± SD (n=3). Student’s t test was used for statistical analysis (*, P ≤ 0.05 and **, P ≤ 0.01).
Figure 1.
Phenotypic of grain filling rate related in foxtail millets. (a) and (b) Morphological changes in foxtail millet grain of CS13 (a) and CN47 (b) at six stages. Bars=500 μm. (c) Thousand millets weight (de-husked grains) changes of CS13 and CN47 at six stages. (d) Grain-filling rate changes of CS13 and CN47 at six stages. (e) Comparison of the agronomic traits between CS13 and CN47. The data are means ± SD (n=3). Student’s t test was used for statistical analysis (*, P ≤ 0.05 and **, P ≤ 0.01).
Figure 2.
Transcriptomes analysis and the total number of differentially expressed genes of the two foxtail millets at two stages during grain filling. (a) Principal component analysis (PCA) plot analysis of transcriptome; the x-axis and y-axis represent principal component 1 (PC1) and principal component 2 (PC2), respectively. S1–S3 represents the three biological replicates of S1 (14 DAF of CS13), S4–S6 represents the three biological replicates of S4 (21 DAF of CS13), N1–N3 represents the three biological replicates of N1 (14 DAF of CN47), N4–N6 represents the three biological replicates of N4 (21 DAF of CN47). (b) Statistic of expressed genes (DEGs) in the CS13_14 vs CN47_14, CS13_21 vs CN47_21 comparison.
Figure 2.
Transcriptomes analysis and the total number of differentially expressed genes of the two foxtail millets at two stages during grain filling. (a) Principal component analysis (PCA) plot analysis of transcriptome; the x-axis and y-axis represent principal component 1 (PC1) and principal component 2 (PC2), respectively. S1–S3 represents the three biological replicates of S1 (14 DAF of CS13), S4–S6 represents the three biological replicates of S4 (21 DAF of CS13), N1–N3 represents the three biological replicates of N1 (14 DAF of CN47), N4–N6 represents the three biological replicates of N4 (21 DAF of CN47). (b) Statistic of expressed genes (DEGs) in the CS13_14 vs CN47_14, CS13_21 vs CN47_21 comparison.
Figure 3.
(a) and (b) The top 25 Gene Ontology (GO) classification in each category of DEGs in the CS13_14 vs CN47_14, CS13_21 vs CN47_21 comparison. (c) and (d) The top 20 Kyoto Encyclopedia of Genes and Genomes (KEGG) classification in each category of DEGs in the CS13_14 vs CN47_14, CS13_21 vs CN47_21 comparison.
Figure 3.
(a) and (b) The top 25 Gene Ontology (GO) classification in each category of DEGs in the CS13_14 vs CN47_14, CS13_21 vs CN47_21 comparison. (c) and (d) The top 20 Kyoto Encyclopedia of Genes and Genomes (KEGG) classification in each category of DEGs in the CS13_14 vs CN47_14, CS13_21 vs CN47_21 comparison.
Figure 4.
Metabolomic profiling analysis at two stages during grain filling of the two foxtail millets. (a) and (b) Volcano plot of DAMs in the CS13_14 vs CN47_14, CS13_21 vs CN47_21 comparison, respectively. (c) and (d) Heatmap analysis of DAMs in the CS13_14 vs CN47_14, CS13_21 vs CN47_21 comparison, respectively; red and blue means high and low level. (e) and (f) Top 20 KEGG pathway classification of DAMs in the CS13_14 vs CN47_14, CS13_21 vs CN47_21 comparison, respectively.
Figure 4.
Metabolomic profiling analysis at two stages during grain filling of the two foxtail millets. (a) and (b) Volcano plot of DAMs in the CS13_14 vs CN47_14, CS13_21 vs CN47_21 comparison, respectively. (c) and (d) Heatmap analysis of DAMs in the CS13_14 vs CN47_14, CS13_21 vs CN47_21 comparison, respectively; red and blue means high and low level. (e) and (f) Top 20 KEGG pathway classification of DAMs in the CS13_14 vs CN47_14, CS13_21 vs CN47_21 comparison, respectively.
Figure 5.
Co-mapped pathway of DAMs and DEGs at two stages during grain filling of the two foxtail millets. (a) and (b) KEGG pathways analysis of the simultaneous annotation of DAMs and DEGs in the CS13_14 vs CN47_14, CS13_21 vs CN47_21 comparison, respectively. (c) and (d) Venn Diagram of DEGs and DAMs in the CS13_14 vs CN47_14, CS13_21 vs CN47_21 comparison, respectively. (e) The number of shared DAMs and DEGs related to co-mapped pathway.
Figure 5.
Co-mapped pathway of DAMs and DEGs at two stages during grain filling of the two foxtail millets. (a) and (b) KEGG pathways analysis of the simultaneous annotation of DAMs and DEGs in the CS13_14 vs CN47_14, CS13_21 vs CN47_21 comparison, respectively. (c) and (d) Venn Diagram of DEGs and DAMs in the CS13_14 vs CN47_14, CS13_21 vs CN47_21 comparison, respectively. (e) The number of shared DAMs and DEGs related to co-mapped pathway.
Figure 6.
The genes expression of auxin biosynthesis, transport and signal transduction pathway. The rectangles color represent the expression level of genes, blue means low level, red means high level.
Figure 6.
The genes expression of auxin biosynthesis, transport and signal transduction pathway. The rectangles color represent the expression level of genes, blue means low level, red means high level.
Figure 7.
The genes expression of cytokinin (CTK) biosynthesis, transport and signal transduction pathway. The rectangles color represent the expression level of genes, blue means low level, red means high level.
Figure 7.
The genes expression of cytokinin (CTK) biosynthesis, transport and signal transduction pathway. The rectangles color represent the expression level of genes, blue means low level, red means high level.
Figure 8.
Quantitative real-time PCR (qRT-PCR) validation and RNA-seq data of 12 selected DEGs.
Figure 8.
Quantitative real-time PCR (qRT-PCR) validation and RNA-seq data of 12 selected DEGs.