Transcriptome Sequencing and Metabolome Analysis Reveals the Regulatory and Molecular Mechanism of Grain Filling Rate in Foxtail Millet (Setaria italica L.)

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.

2. Materials and Methods

2.1. Plant Materials and Growth Conditions

The foxtail millet cultivars (Setaria italica L.) utilized in this study were Changsheng 13 (CS13, with a low grain filling rate) and Changnong No. 47 (CN47, with a high grain filling rate), which were provided by the Millet Research Institute of Shanxi Agricultural University/Shanxi Academy of Agricultural Sciences (Changzhi, Shanxi, China). Foxtail millets were cultivated at the experimental station of the Millet Research Institute of Shanxi Agricultural University/Shanxi Academy of Agricultural Sciences in Changzhi, Shanxi, China (36°21’N, 113°14’E) from May to October in 2023. The experimental field features a temperate semi-humid continental monsoon climate, characterized by an average annual temperature ranging from 8.6 °C to 10.5 °C, an average annual rainfall exceeding 600 mm, and a frost-free period lasting between 155 and 184 days. A randomized block design is employed for the field experiments with each plot measuring 5 m × 2.4 m and three biological replicates. The soil consisted of brown loam with nutrient contents as follows: organic matter (33.51 g/kg), available nitrogen (46.14 µg/g), available phosphorus (20.14 µg/g), and available potassium (168.76 µg/g). Regular agronomic practices were implemented in the field.

2.2. Investigation of the Grain Development and Grain-Filling Rate

The experimental plots were labeled with sixty plants exhibiting spikelets of the same age and size in the central sections, all flowering synchronously. On each plant, 200 flowers were marked with a black pen on the glume surface. After seven days of marking, fresh grains from five randomly selected plants were harvested and subjected to drying in an oven at 105°C for 30 minutes, followed by further drying in another oven at 80°C until reaching a constant weight. Subsequently, these dried grains were manually hulled and weighed using a balance with a precision of up to 0.001 g. Two hundred foxtail millet grains chosen randomly per plant were weighed and averaged. Similarly, fresh grains collected at 7, 14, 21, 28, 35, and 42 days after marking underwent harvesting, drying, hulling, and weighing procedures. The grain-filling rate was then calculated for each stage using the formula: GFRi = (GWi - GWi-1) /7 , where GFRi represents the grain-filling rate (mg grain-1 d-1), GWi denotes the dried grain weight (mg grain-1 d-1) at stage i (i = 1, 2, 3, 4, 5). A stereo microscope (SZX10, Olympus, Japan) was employed to document changes in developing grains at different stages including 7, 14, 21, 28, 35, and 42 DAF(days after flower) [32].

In addition, grains were sampled at the 14, and 21 DAF from of the same age and size in the central sections of marked plant. Grains from five spikelets within the same field plot were pooled together as one biological replicate. Three biological replicates were selected for each developmental stage from three different field plots. All the samples were immediately frozen in liquid nitrogen and stored at -80 °C for subsequent experiments.

2.3. RNA Sequencing and Transcriptomic Profiling Analysis

The grains of CS13 and CN47 at the 14 DAF stage and 21 DAF stage were collected for RNA extraction and sequencing. RNA-seq was performed by Novogene Co., Ltd. (Beijing, China). Each treatment group consisted of three biological replicates. The integrity of RNA was assessed using the Bioanalyzer 2100 system (Agilent Technologies, CA, USA) with the RNA Nano 6000 Assay Kit. The cDNA library was constructde and sequencing were conducted using the Illumina (HiSeq/MiSeq) system by Genepioneer. Subsequently, clean reads were obtained by filtering out unqualified raw reads, based on Q20, Q30 and GC content criteria. Trimmomatic (version 0.36) [41] was used to trim the clean reads while HISAT2 (version v2.0.5) [42] mapped them to the reference genome of foxtail millet [39]. Differential expression analysis of genes was performed using DESeq package with BaseMean value estimation [43]. Differentially expressed genes (DEGs) were identified based on | log2(fold change) | ≧ 1 and adjusted P-value < 0.05 as screening criteria. Gene Ontology (GO) [44] and Kyoto Encyclopedia of Genes and Genomes (KEGG) [45] pathway enrichment analyses were performed by the Novogene cloud analysis tools ( https://magic.novogene.com/ ).

2.4. Metabolome Profiling Analysis

The metabolome analysis material was identical tp that used for transcriptome analysis, with six biological replicates per treatment group (Metabolite extraction and LC-MS analysis also were also conducted by Novogene Co., Ltd. in Beijing, China). Tissues weighing 100 mg each were individually ground using the liquid nitrogen grinding method, and the resulting homogenate was resuspended in prechilled 80% methanol throuth vigorous vortexing. The samples were then incubated on ice for 5 min before being centrifuged at 15,000 g, 4°C for 20 minutes. A portion of the supernatant was diluted with LC-MS grade water to a final concentration containing 53% methanol. Subsequently, the samples were transferred to a fresh Eppendorf tubes and centrifuged again at 15000 g, 4°C for 20 minutes. Finally, the supernatant was injected into the LC-MS/MS system for analysis [46]. UHPLC-MS/MS analyses were performed using a Vanquish UHPLC system (ThermoFisher, Germany) coupled with an Orbitrap Q ExactiveTM HF mass spectrometer or Orbitrap Q ExactiveTMHF-X mass spectrometer (Thermo Fisher, Germany) at Novogene Co., Ltd. Samples were injected onto a Hypersil Gold column (100×2.1 mm, 1.9μm) using a linear gradient over a duration of 12 minutes at a flow rate of 0.2 mL/min. The column temperature was maintained to 40 °C throughout the progress. For positive polarity mode, the eluents consisted of eluent A (0.1% FA in Water) and eluent B (Methanol). For negative polarity mode, the eluents comprised eluent A (5 mM ammonium acetate, pH 9.0) and eluent B (Methanol). The solvent gradient was set as follows: 2% B, 1.5 min; 2-85% B, 3 min; 85-100% B, 10min; 100-2% B, 10.1 min; 2% B, 12 min. The Q ExactiveTM HF mass spectrometer was operated in positive/negative polarity mode with spray voltage of 3.5 kV, capillary temperature of 320 °C, sheath gas flow rate of 35 psi and aux gas flow rate of 10L/min, S-lens RF level of 60, Aux gas heater temperature of 350°C.

The raw data files generated by UHPLC-MS/MS were processed using Compound Discoverer 3.3 (CD3.3, ThermoFisher) for peak alignment, peak picking, and metabolite quantitation. Metabolites quantification involved integrating the peak area from the mass spectra of all the substances and correcting for variations in mass spectrometry peaks across different samples to ensure the accurate quantification. Metabolites were subjected to principal component analysis (PCA) and partial least squares discriminant analysis [47]. The differences in metabolites between treatments were analyzed using a t-test. The metabolites with VIP > 1 and P-value< 0.05 among with | log2 (fold change) | > 0.5 were considered as differential metabolites (DAMs). These DAMs were annotated using the KEGG database (https://www.genome.jp/kegg/pathway.html.) and mapped onto KEGG metabolic pathways for pathway analysis and enrichment.

2.5. Conjoint Analysis of Transcriptome and Metabolome

By simultaneously mapping DEGs and DAMs to the KEGG Pathway database, we obtained their shared pathway information. Correlation analysis was performed using quantitative values of genes and metabolites in all samples. Spearman’s correlation coefficient (Cor) and p-value were calculated using the Cor function of R language. Finally, a network map of DEGs-DAMs based on Spearman’s correlation coefficient (with |Cor| > 0.8 and p ≤ 0.05) was constructed using Cytoscape.

2.6. Quantitative qRT-PCR

Total RNA extracted by the TaKaRa MiniBEST Plant Total RNA Extraction Kit from all samples. The SynScript®Ⅲ RT SuperMix for qPCR was used for reverse transcription of RNA, and the synthesized first-strand cDNA was diluted 5 times and was used as a template for qRT-PCR. qPCR was carried out using the ArtiCanCEO SYBR qPCR Mix and Bio-Rad CFX96 (Bio-Rad CFX96, BIORAD, USA). The reaction conditions were as follows: pre-denaturation at 95 ℃ for 5 min, 40 cycles of denaturation at 95 ℃ for 15 s, and annealing, extension at 60 ℃ for 20 s, and extension at 72 ℃ for 20 s. The relative expression of genes was calculated using the 2^−∆∆Ct method. Each pair of primers was set up with three biological replicates and three technical replicates. The gene expression was normalized to housekeeping gene Actin 1. The primers were listed in Table S8.

2.7. Statistical Analysis

The statistical analysis for the experiment on phenotype and qRT-PCR as well as content verification of DAMs, was performed using SPSS software. A significantly level of p ≤ 0.05 was considered to indicate a statistically significant difference in the t-test. The final graph was generated using Origin software.

3. Results

3.1. Characterization of Grain Filling Rate-Related Phenotypes in Foxtail Millet

To investigate dynamic changes in foxtail millet grain development after flowering, spikelets were collected from two cultivars (CS13 and CN47) at 7, 14, 21, 28, 35 and 42 days after flowering (DAF) during the grain filling period. The growth morphology changes of de-husked grains of foxtail millet was observed across six stages using a stereo microscope (Figure 1a,b). It was observed that CN47 displayed significant differences in grain size and millet color compared to CS13 at 14 DAF and 21 DAF stages. Additionally, grain-filling rates and thousand millets weight were measured for both cultivars at each time point mentioned (Figure 1c,d). There were distinct variations in GFR and thousand foxtail millet weight between CS13 and CN47, with CN47 consistently exhibiting higher thousand foxtail millet weight from 14 DAF to 42 DAF. The highest grain-filling rates were observed at 21 DAF for both cultivars, with CN47 exhibiting higher rates at 14 DAF and 21 DAF compared to CS13. No significant differences were found among the agronomic traits analyzed, including heading date (HD), plant height (PH), grain width (GW), grain number per panicle, and panicle number per plot. However, significant differences were observed in grain length (GL), grain weight per panicle, thousand-grain weight (TGW), and yield per plot (YPP) between CS13 and CN47 (Figure 1e). Analysis of agronomic traits indicated that the relatively high values of GL, TGW, and YPP in CN47 could be attributed to its elevated grain- filling rate.

3.2. Transcriptomic Difference in the Grains during Grain-Filling Stages between Foxtail Millet Cultivars

Transcription profiles of CS13 and CN47 were examined at the 14 DAF and 21 DAF stages, with three biological replicates used for each stage. Twelve cDNA libraries were sequenced separately by Illumina sequencing. An average of 6.7 Gb clean bases per library were generated after removing the adapter sequences, low quality and contaminated reads, with a mean Q30 level of 94.8%. This indicated that the library contained high-quality raw reads (Table S2). The read mapping ratio to the reference genome per sample ranged from 92% to 94% (Table S3). Gene expression levels were quantified using fragments mapped (FPKM). Principal component analysis (PCA) revealed distinct clustering of the samples, with tight clusters formed by the three biological replicates within each group (Figure 2a). Differential expression analysis was performed to compare CS13_14 DAF vs CN47_14 DAF and CS13_21 DAF vs CN47_21 DAF, with a fold change threshold of > 1 and false discovery rate FDR < 0.05. The CS13_14 DAF vs CN47_14 DAF comparison revealed a total of 765 differentially expressed genes (DEGs), with 439 genes up-regulated and 326 genes down-regulated. Similarly, 908 DEGs were identified in the CS13_21 DAF vs CN47_21 DAF comparison, with 471 up-regulated genes and 439 down-regulated genes (Figure 2b).

The GO enrichment analyses were performed to gain further insights into the functions and associated biological processes in which the DEGs participated, with the results depicted in Figure 3. The GO analysis revealed that the DEGs were classified into cellular components (CC), molecular functions (MF), and biological processes (BP). We further analyzed the top 10 categories for CC, MF, and BP in the GO annotation (Figure 3a,b). In terms of molecular functions, we observed an enrichment of shared GO terms such as ADP binding, nutrient reservoir activity and cysteine-type peptidase activity. Regarding cellular components, both comparison groups showed enrichment of GO terms related to the DNA-directed RNA polymerase II, holoenzyme, transcription factor complex, DNA-directed RNA polymerase complex, nuclear DNA-directed RNA polymerase complex and RNA polymerase complex. Different GO terms were enriched in two comparision groups for biological processes, with only two shared terms identified. These findings indicated that various factors influence the differences in grain-filling rates between CS13 and CN47.

To reveal the biological pathways associated with GFR, the pathway enrichment analysis of the DEGs was also performted by the Kyoto Encyclopedia of Genes and Genomes (KEGG), with the results depicted in Figure 3c,d. The KEGG enrichment analysis revealed that 12 out of the top 20 terms were shared between the two stages, indicating a high consistency in the pathways involved in grain-filling. We observed significant enrichment in pathways including basal transcription factor, sphingolipid metabolism, zeatin biosynthesis, plant-pathogen interaction, motor proteins, ribosome biogenesis in eukaryotes, nucleotide metabolism, fatty acid degradation, biotin metabolism, tyrosine metabolism, fatty acid biosynthesis, fatty acid metabolism. Moreover, the pathways of plant hormone signal transduction and nitrogen metabolism were identified as pivotal for the process of grain-filling. These findings provided further evidence that the regulation of grain-filling rate involved multiple pathways.

3.3. Metabolite Difference in the Grains during Grain-Filling Stages between Foxtail Millet Cultivars

In order to identify differential metabolites between the two cultivars at 14 DAF and 21 DAF stages, a non-targeted metabolic analysis was conducted,which revealed distinct patterns of metabolite accumulation. Metabolites were primarily categorized into two modes: the positive ion mode (POS) and the negative ion mode (NEG). Metabolites exhibiting with | log2 (fold change) | > 0.5 and VIP ≥ 1 were chosen as differentially accumulated metabolites (DAMs). In the POS and NEG modes, a total of 246 DAMs (146 up-regulated, 100 down-regulated) were identified between CS13_14 vs CSN47_14 DAF, while 268 DAMs (133 up-regulated, 135 down-regulated) were found between CS13_21 DAF vs CSN47_21 DAF (Table S4; Figure 4a,b). The gene clustering heatmap analysis demonstrated excellent biological repeatability, with clear distinctions among the samples (Figure 4c,d). At 14 DAF stage, 44 DAMs were annotated to the KEGG pathway and mapped to 24 pathways (Table S5). At 21 DAF stage, 40 DAMs were annotated to the KEGG pathway and mapped to 23 pathways (Table S6). The top 20 pathway were depicted in Figure 4e,f. A large number of identical pathways have been discovered at 14 DAF and 21 DAF stages. These included flavonoid biosynthesis, biosynthesis of secondary metabolites, flavone and flavonol biosynthesis, carotenoid biosynthesis, isoflavonoid biosynthesis, isoquinoline alkaloid biosynthesis, betalain biosynthesis, plant hormone signal transduction, tyrosine metabolism, sesquiterpenoid and triterpenoid biosynthesis, caffeine metabolism, anthocyanin biosynthesis, sulfur relay system, glutathione metabolism, ascorbate and aldarate metabolism, ATP-binding cassette (ABC) transporters, beta-Alanine metabolism, and stilbenoid, diarylheptanoid and gingerol biosynthesis. The result suggested that these metabolic pathways exhibited similarly between the two stages during grain filling.

3.4. Integrated Transcriptome and Metabolome Analysis

Based on the results of DAMs and DEGs, a joint analysis of the metabolome and transcriptome was then conducted to unveil their relationship and regulatory mechanism. Figure 5 displayed the KEGG pathways that were co-enriched in both DAMs and DEGs during 14 DAF and 21 DAF stages. The results showed that flavonoid biosynthesis, plant hormone signal transduction, tyrosine metabolism, ABC transporters, beta-Alanine metabolism, stilbenoid, diarylheptanoid and gingerol biosynthesis were identified as shared pathways in both comparisons of CS13_14 DAF vs CSN47_14 DAF and CS13_21 DAF vs CSN47_21 DAF (Figuure 5a, 5b). However, the DAMs and DEGs mapped to shared pathway were differed at different stages (Table S7).

To identify key genes related to grain-filling rate, a Venn diagram analysis was conducted to compare all DEGs and DAMs at the 14 DAF and 21 DAF stages of two foxtail millet cultivars. As shown in the Figure 5c, there were 512 identical DEGs for the comparisons of CS13 and CSN47, with 274 genes up-regulated and 238 genes down-regulated showing a similar trend of change, accounting for more 50 percent. At the same time, we found 116 identical metabolites, including 66 up-regulated metabolites and 55 down-regulated metabolites (Figure 5d). The change trend of DAMs were also the same. Further analysis reveals that shared DEGs and DAMs were enriched the same pathway as mentioned above (Table S8; Figure 5e). It can be inferred that flavonoid biosynthesis, plant hormone signal transduction, tyrosine metabolism, ABC transporters, beta-alanine metabolism, stilbenoid, diarylheptanoid and gingerol biosynthesis may be related to the high GFR of the CN47.

3.5. Many Hormones Were Key Regulators of the Grain Filling Rate

KEGG pathway classification analysis was performed on shared DEGs and DAMs, revealing a significant enrichment of share DEGs in hormone signaling. Our results have identified 11 genes involved in phytohormones signaling, including auxin, cytokinin (CTK), and ethylene (ETH) signaling pathways. To investigate the expression patterns of phytohormone-related genes, we exmined gene expression profile clustering. Among the phytohormones, auxin played a crucial role in grain filling stages [48]. We found several key genes enriched in the auxin biosynthesis (tryptophan metabolism), transport (ABC transport) and signal transduction pathway at grain filling 14 DAF and 21 DAF stages. And a network was established to display the relationship between auxin and plant growth and development (Figure 6). It is worth noting that the tryptophan-dependent pathway is the primary pathway responsible for the production of indole-3-acetic acid (IAA) [49]. SETIT_026842mg was identified, which encoded a YUCCA protein associated with tryptophan metabolism that acted as a rate-limiting enzyme for auxin biosynthesis. The gene was expressed at low levels in CN47 at 14 DAF and 21 DAF stages, it may negatively regulate GFR since the IAA concentration of CN47 was higher than CS13. And two genes related to the ABCB1 were identified, namely SETIT_032403mg and SETIT_000076mg. They were linked to ABC transport and coded an ABCB protein that interacted with PIN protein to regulate Polar Auxin Transport. In CN47, SETIT_032403mg showed low expression level at the 14 DAF and 21 DAF stages of grain filling progress. Conversely, the expression patterns of SETIT_000076mg were different. Additionally, we identified several key genes involved in the auxin signal transduction pathway, including three AUX1, one AUX/IAA, and one SAUR, SETIT_0218801mg, SETIT_0218802mg, SETIT_035200mg, SETIT_007373mg and SETIT_028123mg. Compared with the 14 DAF, the downregulated genes SETIT_0218801mg, SETIT_0218802mg, SETIT_035200mg and SETIT_007373mg showed high expression differences at 21 DAF in CN47 compared to CS13. The small Auxin-Up RNA (SAUR) families constituted the most important family of auxin-responsive proteins, potentially playing an important role in the regulation of foxtail millet grain development [50]. Furthermore, SETIT_028123mg, which was highly expressed in CN47 at 14 DAF and 21 DAF stages, was identified related to SUAR. It may be have a positive correlation with the GFR of foxtail millet. These results suggested their potential role in facilitating extracellular-to-intracellular auxin transport. Notably, SETIT_007373mg and SETIT_028123mg, as early response genes to auxin, could regulate downstream gene expression. These data suggest that these genes related to auxin signaling are involved in grain filling progress and GFR. Therefore, the potential function of these key genes related to auxin biosynthesis, transport, and signal transduction underscores the significance of IAA during grain filling and GFR.

Cytokinin (CTK), an important phytohormone, plays a significant role in grain filling progress [51]. The CTK biosynthesis (Zeatin metabolism) and signal transduction pathway exhibited enrichment of several genes during grain filling stages at 14 DAF and 21 DAF, as depicted in Figure 7. Furthermore, a network was established to illustrate the correlation between CTK and plant growth and development. It was noteworthy that the Zeatin metabolism pathway played an important role for the production of CTK. Their genes include a prenyltransferase (SETIT_008985mg), and a cytochrome monooxygenase (SETIT_032720mg), which respectively encodes IPTs protein and CYP735A protein associated with zeatin metabolism. These enzymes act as rate-limiting enzymes for CTK biosynthesis. The expression level of SETIT_008985mg was significantly upregulated at 14 DAF and 21 DAF stages in CN47, while the CTK contents were also higher in CN47 compared to CN13 at these developmental stages. This suggested that SETIT_008985mg may be positively correlated with the GFR of foxtail millet. Additionally, we identified three key genes enriched in the CTK signal transduction pathway during grain filling stages at 14 DAF and 21 DAF, These genes were involved in SETIT_023526mg, SETIT_018243mg, and SETIT_011287mg. Specifically, SETIT_023526mg was identified as encoding the AHPs protein which facilitated the transfer of phosphate groups to response regulators (SETIT_018243mg, SETIT_011287mg), thereby inducing gene expression and regulating plant growth and development. This gene exhibited high expression levels in CN47 at 14 DAF and 21 DAF stages, suggesting a potential positive correlation with foxtail millet's GFR. Additionally, SETIT_018243mg and SETIT_011287mg, as members of the type-A (ARR) family of response regulators, displayed distinct expression patterns between CN47 and CS13 during the 14 DAF and 21 DAF stages. These three genes collectively functioned as rate-limiting enzymes for CTK signal transduction. In addition, this study identified SETIT_029213mg and SETIT_009909mg, which were associated with the ethylene-insensitive 3 (EIN3)-like 1 (EIL1) protein and EIN3-like 5 (EIL5) protein respectively (Figure S1). These genes exhibited high expression levels at CN47 during the 21 DAF stage, suggesting may be a positive correlation with GFR.

3.6. qRT-PCR Verification of RNA-seq Data

To ensure the reliability of our transcriptome data, we selected 12 key genes through qRT-PCR validation (Figure 8). These genes were involved in plant hormone signal transduction, flavonoid biosynthesis, ABC transporters, tyrosine metabolism, beta-Alanine metabolism, as well as stilbenoid, diarylheptanoid and gingerol biosynthesis. The results of our qRT-PCR analysis demonstrated that the observed variation in gene expression levels aligned with those obtained from the RNA-seq data. This concurrence strongly suggests the reliability of our analytical findings.

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.

5. Conclusions

In this study, we conducted a comparative analysis of the transcriptome and metabolome at 14 DAF and 21 DAF stages to investigate the molecular mechanism regulating GFR in CS13 (low GFR) and CN47 (high GFR). The results revealed that a total of 765 DEGs and 908 DEGs, as well as 246 DAMs and 268 DAMs, were identified at 14 DAF and 21 DAF stages both in CS13 and CN47 respectively. The functional analysis of DEGs and DAMs revealed stage-specific properties in the GFR of foxtail millet. Several important pathways, such as flavonoid biosynthesis, plant hormone signal transduction, tyrosine metabolism, ABC transporters, beta-Alanine metabolism, as well as Stilbenoid, diarylheptanoid, and gingerol biosynthesis were enriched via a conjoint analysis of DEGs and DAMs. Notably most of the share DEGs were enriched in hormone signaling pathways, including auxin, cytokinin, and so on. Especially, several key genes related to auxin biosynthesis (YUCCA), transport (ABCB1) and signal transduction (AUX1, AUX/IAA, SAUR) pathway were identified in the ABC transporters and plant hormone signal transduction pathway. These findings highlight the significant role of auxin in regulating GFR and grain filling progress, providing valuable molecular insights for investigating differences in GFR among foxtail millet cultivars and theoretical mechanisms for high-yield breeding.