Abstract To elucidate the mechanisms underlying photoperiodic responses, we investigated the genomic and metabolomic responses of two broomcorn millet (Panicum miliaceum L.) genotypes. For this purpose, light-insensitive (D32) and light-sensitive (M51) genotypes were exposed to a 16 h photoperiod (long-day (LD) conditions) and an 8 h photoperiod (short-day (SD) conditions), and various transcriptomic and metabolomic changes were investigated. A total of 1664, 2564, 13,017, and 15548 DEGs were identified in the SD-D, LD-D, LD-M, and SD-M groups, respectively. Furthermore, 112 common DEGs were identified as well. Interestingly, most DEGs in the different groups were associated with photosynthesis and phenylpropanoid and carotenoid biosynthesis. In addition, 822 metabolites were identified under different treatments. The main metabolites, including L-malic and fumaric acids, were identified in the negative mode, whereas brucine and loperamide were identified in the positive mode. KEGG analysis revealed that the metabolites in the different groups were enriched in the same metabolic pathway of the TCA cycle. Furthermore, in negative mode, the metabolites of M51 were mainly D-glucose, whereas those of D32 were mainly L-malic and fumaric acids. One photoperiod candidate gene (C2845_PM11G01290), annotated as ATP6B, significantly increased the levels of L-malic and fumaric acids. In conclusion, our study provides a theoretical basis for understanding the molecular mechanisms of photoperiodic response regulation and can be used as a reference for marker development and resource identification in Panicum miliaceum L.. Keywords: Broomcorn millet, Panicum miliaceum L., Photoperiodism, Metabolomics, Transcriptomics Subject terms: Genetics, Molecular biology Introduction Broomcorn millet (Panicum miliaceum L.) is a short-day cereal crop that is both drought-resistant and barren tolerant^[34]1,[35]2. In contrast, it is highly sensitive to the photoperiod, which affects its adaptability to different geographical regions. Broad adaptability is a major factor that increases the crop acreage. During various growth and developmental stages of different crops, the photoperiod is involved in many complex physiological and biochemical processes, and plants can respond to changes in a timely manner to maintain their growth and development at a relatively optimum state^[36]3,[37]4. One key photoperiod-dependent developmental stage is flowering, which involves complex regulatory mechanisms that ensure the transition from the vegetative to the reproductive stage. Six primary pathways are involved in plant floral regulation: vernalization, photoperiod, gibberellin, autonomous, temperature-sensitive, and age pathways^[38]5. However, the photoperiod pathway is the most conservative flowering response pathway in plants and is the key to regulating flowering. So far, photoperiod-induced flowering mechanisms of Oryza sativa L. and Arabidopsis thaliana have been studied. The OsGI(GI)-Hd1(CO)-Hd3a/RFT1(FT) pathway is a common pathway that promotes flowering in Oryza sativa L. and Arabidopsis thaliana, and at least 64 related genes have been identified^[39]6. More than 600 Quantitative Trait Locis (QTLs) related to the heading date have been identified in Oryza sativa L., and dozens of genes have been screened. Among the cloned genes, approximately 13 genes were linked to the photoperiod pathway^[40]7, including Hd1/Se-1^[41]8, Hd2/Ghd7.1/DTH7/OsPRR37^[42]9,[43]10, Hd3a^[44]11, RFT1^[45]12, Hd4/Ghd7^[46]13, Hd5/Ghd8/DTH8/LHD1^[47]14,[48]15, Hd6/CK2a^[49]16, Ehd1^[50]17, Hd16/EL1^[51]18, Hd17/Ef7^[52]19, Ehd2/RID1^[53]20, Ehd3^[54]21, Ehd4^[55]22, and DTH2^[56]23. Homology analyses of photoperiod genes in different crops revealed that some photoperiod genes are conserved. For instance, Oryza sativa L. genes Hd3a^[57]24, Hd1^[58]25, OsGI^[59]26, Osphy^[60]27, and OsELF3^[61]28 are homologous to FT, GI, PHY, and ELF3 of Arabidopsis thaliana, respectively. Triticum aestivum L. Ppd-1^[62]29, Hordeum vulgare L. Ppd-H1^[63]30, Sorghum bicolor L. SbPRR37^[64]31, and Setaria italica SiPRR37^[65]32 are homologous to AtPRR3 in Arabidopsis thaliana and OsPRR37 in Oryza sativa L.. Additionally, the regulatory pathways of different crop photoperiods on the heading date, such as Ghd7-Ehd1-Hh3a/RFT1^[66]13,[67]33 and PHYA-E1/E1L-FT^[68]34–[69]36, are unique to Oryza sativa L. and Glycine max, respectively. In addition, different crops have specific photoperiod genes, such as Oryza sativa L. Ghd7^[70]13, Ehd1^[71]17, Ehd2^[72]20, Ehd3^[73]21, Ehd4^[74]22, and Glycine max E1^[75]34,[76]35. However, the photoperiod pathway-related genes and their mechanisms remain unknown in Panicum miliaceum L.. Light directly regulates plant physiology and development and influences the response of many plants to their environment^[77]37. Light can regulate plant responses through two mechanisms: photosynthesis and photoreceptor signaling^[78]38. Photosynthesis provides reducing equivalents and energy for the production of metabolites, and plant photoreceptors sense the surrounding red/far red, blue/UV-A, and UV-B light^[79]38. Different photoperiods affect plant metabolite production. For example, blue/UV-A induced higher levels of anthocyanin accumulation in Lactuca sativa L., Fragaria ananassa, and Vitis vinifera L.^[80]39 and increased the contents of chlorogenic, gallic, and ferulic acids in Pisum sativum Linn seedlings^[81]40. Red/far red could promote stem growth by primarily regulating gibberellin biosynthesis in "Manicure fingers" Vitis vinifera L. plantlets^[82]41, and exposure to red and blue light for three months stimulates carotenoid biosynthesis, which is a sign of a slower decay of product compounds after Solanum lycopersicum harvest^[83]42. Moreover, a previous study used metabolomic approaches to show that alkaline stress can alter the biosynthesis of phenylpropanes, flavonoids, flavonols, valine, leucine, and isoleucine, as well as the metabolism of arginine, proline, tryptophan, and ascorbic acid in Panicum miliaceum L.^[84]43. However, the effects of photoperiod on the changes in metabolites in Panicum miliaceum L. remain unclear. Therefore, this study aimed to investigate the physiobiochemical basis of the photoperiodic responses of Panicum miliaceum L. and the key genes involved. In this study, two types of Panicum miliaceum L. genotypes (photoinsensitive: D32 and photosensitive: M51) were subjected to short-day (8 h) and long-day (16 h) conditions, and transcriptomics and metabolomics were employed to explore the core differentially expressed genes (DEGs) and differential metabolites that regulate photoperiodism in Panicum miliaceum L.. Materials and methods Experimental details Two Panicum miliaceum L. genotypes (M51, photoperiod-sensitive and D32, photoperiod-insensitive) were obtained from the germplasm bank of the Center for Agricultural Genetic Resources Research, Shanxi Agricultural University (37.4244° N, 112.5784° E), China. The experiment was conducted under controlled environmental conditions in a growth incubator (temperature 25 °C; humidity 60%; 16 h as a long-day conditions and 8 h as a short-day photoperiod; ARC-41L2-AR (Percival Scientific, Inc.), the phenotypes of M51 and D32 were significantly different under short-day (8 h) and long-day (16 h) photoperiod treatments (Supplementary Figure [85]S1). Small plastic pots (10 mm in diameter) were filled with the substrate, and ten seeds were sown in each pot. After germination, only five seedlings were retained per pot and were exposed to different photoperiods. Four treatments were established: LD-D (D32 Panicum miliaceum L. genotype under long-day treatment), LD-M (M51 Panicum miliaceum L. genotype under long-day treatment), SD-D (D32 Panicum miliaceum L. genotype under short-day treatment), and SD-M (M51 Panicum miliaceum L. genotype under short-day treatment). Each treatment was repeated five times in a total of 20 pots. After that, samples were taken at the beginning of the tricuspid phase, as well as were rapidly cut with scissors, transferred to centrifuge tubes, and rapidly transferred to liquid nitrogen, stored at − 80 °C, and used for the subsequent transcriptomic and metabolomic analyses. Transcriptome sequencing Total RNA was extracted from the seedlings and RNA integrity was detected using an Agilent Bioanalyzer 2100 system. RNA sequencing (RNA-seq) was performed by Shanghai Personal Biotechnology Co., Ltd. (Shanghai, China). High-quality total RNA samples were reverse transcribed into cDNA and used to construct cDNA libraries. In addition, total RNA samples were subjected to next-generation sequencing (NGS), and the quality of the reads was checked using FastQC Fastp to remove adapters and low-quality reads. HISAT2 ([86]https://daehwankimlab.github.io/hisat2/) software was used to align the filtered reads with non-redundant protein (NR), gene ontology (GO), uniport knowledgebase (UniProt), kyoto encyclopedia of genes and genomes (KEGG), and evolutionary genealogy of genes: non-supervised orthologous groups (eggNOG) databases. Gene expression levels were normalized to fragments per kilobase of transcript per million fragments (FPKM), and genes with FPKM > 1 were expressed. Differential analysis of gene expression was performed using the DESeq package in R software, and |log[2] fold change (FC)| > 1 and P < 0.05 were used as thresholds for screening DEGs^[87]44–[88]46. Furthermore, gene ontology (GO)and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed using the Cluster Profiler package in R software. Extraction of metabolites and metabolomics analysis of the seedlings in the different groups Seedlings in the different groups (n = 3) were ground with liquid nitrogen, and 100 mg samples were added to 1 mL of pre-cooled methanol acetonitrile aqueous solution (v: v: v = 2:2:1). After vortex oscillated and low temperature ultrasound for 30 min, the samples were placed at − 20 °C for 60 min, and then centrifuged at 14000 g for 20 min at 4 °C. The supernatant was concentrated to dry in vacuum, and the dried powder was redissolved in 100 μL acetonitrile aqueous solution (v: v = 1:1). After vortexing and centrifugation at 14000 × g for 15 min at 4 °C, the supernatant was used for ultra-high-performance liquid chromatography-mass spectrometry (UHPLC-MS) detection. Analyses were performed using an UHPLC (1290 Infinity LC, Agilent Technologies) coupled to a quadrupole time-of-flight (AB Sciex TripleTOF 6600). The temperature of sample injector and column was 4 °C and 25 °C, and the flow rate and the injection volume were 0.4 mL/min and 2 μL, respectively. The mobile phase A consisted of 25 mM ammonium acetate and 0.5% formic acid in water, mobile phase B was methanol. The gradient elution procedure was as follows: 0–0.5 min, 5% B; then B changed to 100% linearly from 0.5 to 10min; 10–12. 0 min, B was maintained at 100%; From 12.0 to 12.1 min, B changed linearly from 100 to 5%; 12.1–16 min, B was maintained at 5%. Positive and negative ion modes of electrospray ionization (ESI) were used for MS detection. The ESI source parameters were as follows: gas1 60, gas2 60, curtain gas 30 psi, ion source temperature 600 °C, and spray voltage ± 5500 V (positive and negative modes). The detection ranges for the primary mass charge ratio and secondary mass charge ratio were 60-1000 Da, and 25-1000 Da, respectively. The accumulation time of the primary and secondary MS was 0.20 and 0.05 s/spectra, respectively. The secondary MS was obtained using data-dependent acquisition mode, and the peak intensity value screening mode was used with the cluster-removing voltage of ± 60 V (positive and negative modes), and the collision energy of 35 ± 15 eV^[89]47. The raw MS data (wiff.scan files) were converted to MzXML files using ProteoWizard MSConvert before importing into freely available XCMS software. For peak picking, the following parameters were used: centWave m/z = 10 ppm, peakwidth = c (10, 60), prefilter = c (10, 100). For peak grouping, bw = 5, mzwid = 0.025, minfrac = 0.5 were used. CAMERA (Collection of Algorithms of MEtabolite pRofile Annotation) was sued for annotation of isotopes and adducts. In the extracted ion features, only the variables having more than 50% of the nonzero measurement values in at least one group were kept. Compound identification of metabolites was performed by comparing of accuracy m/z value (<10 ppm), and MS/MS spectra with an in-house database established with available authentic standards. After normalized to total peak intensity, the processed data were analyzed by R package, where it was subjected to multivariate data analysis, including Pareto-scaled principal component analysis (PCA) and orthogonal partial least-squares discriminant analysis (OPLS-DA). The variable importance in the projection (VIP) value of each variable in the OPLS-DA model was calculated to indicate its contribution to the classification. Subsequently, differential metabolites between the two groups were screened based on the thresholds of Variable Importance for the Projection (VIP) > 1 and P < 0.05^[90]48, and the identified differential metabolites were subjected to KEGG pathway enrichment analysis. Combined analysis of transcriptomic and metabolomic data For the combined analysis of transcriptomic and metabolomic data, the identified DEGs and differential metabolites were first extracted based on the results of quantitative detection analysis of the metabolome and transcriptome. Afterwards, according to the metabolite information in the KEGG database, the corresponding DEGs of related enzymes were extracted, and then the metabolites and related DEGs were mapped to related metabolic pathways^[91]49,[92]50. All the DEGs and metabolites with P < 0.05 enriched in the common pathway were submitted for correlation analysis and Pearson correlation coefficient calculation. Interaction pairs with a Pearson correlation coefficient > 0.96 were retained and used to construct a network using Cytoscape. Real-time fluorescent quantitative PCR (qRT-PCR) To validate the accuracy of the photoperiod transcriptome results of Panicum miliaceum L., 10 genes were randomly selected for qRT-PCR, and a candidate gene (C2845_PM11G01290) was chosen for 24 h rhythmic expression level determination using qRT-PCR. Based on the coding gene sequences, the primers of qRT-PCR were designed using primer premier 5.0 software, and all primer sequences were synthesized by Sangon Biotech (Shanghai, China) (Supplementary Table [93]S1). Briefly, total RNA was used to generate cDNA using M-MLV Reverse Transcriptase (TaKaRa Bio, Beijing, China) according to the manufacturer’s instructions. RT-qPCR was performed using an ABI 7300 Real-Time PCR System (Applied Biosystems). The reaction conditions were initiated at 95 °C for 30 s, and a total of 40 cycles at 95 °C for 5 s, 55 °C for 31 s, and 72 °C for 31 s. Actin (No. LOC112893813) was used as a reference gene, and the relative gene expressions of the selected genes were calculated using the 2^−ΔΔCt method^[94]51. The experiments were done in triplicates, and the data were expressed as mean ± standard error (SE). Data analyses One-way analysis of variance (ANOVA) was performed using statistical product and service solutions (SPSS) (version 22.0) software ([95]https://www.ibm.com/spss), and Origin 2018 (version 9.0) software ([96]https://www.originlab.com/) was used to plot the datasets. Both principal component analysis (PCA) and redundancy analysis (RDA) were performed using R (version 3.5.1) software ([97]https://www.r-project.org/); and the KEGG pathway was plotted with the “ggplot2” package ([98]https://ggplot2.tidyverse.org) in R software. Results Transcriptome analysis under different photoperiodic treatments The results showed that in the 12 samples, more than 95.28% of the clean reads were successfully aligned to the reference genome (Shanghai Center for Plant Stress Biology and Center of Excellence in Molecular Plant Sciences, Chinese Academy of Sciences). The alignment efficiencies of these reads to the reference genome were approximately 87.86%. Furthermore, most of the mapped reads were concentrated in the exon region of the gene, indicating that both the sequencing data and the chosen reference genome are of high quality and suitable for analysis (Supplementary Table [99]S2). Expression fold difference |log[2]FC| > 1 and P < 0.05 were used as criteria to screen DEGs for each comparison group. The results showed that a total of 1664 DEGs were identified in different Panicum miliaceum L. genotypes under the short-day treatment (SD), including 920 upregulated and 744 downregulated genes. Similarly, 2564 DEGs were identified under the long-day treatment (LD), including 1449 upregulated and 1115 downregulated genes. A total of 13017 DEGs were identified in the M51 Panicum miliaceum L. genotype under different photoperiod treatments, including 5915 up-regulated and 7102 down-regulated genes, and 15548 DEGs were identified in the D32 Panicum miliaceum L. genotype under different photoperiod treatments, consisting of 7116 upregulated and 8432 downregulated genes. The number of DEGs in the different Panicum miliaceum L. genotypes under the same light treatment was significantly lower than that in the same Panicum miliaceum L. genotype under different light treatments. Therefore, genetic mechanisms may be more responsive to changes in environmental factors (such as light) than to the materials themselves. Analysis using Venn diagrams (Fig. [100]1) showed that there were 544 DEGs specific to LD-D vs. LD-M, 3822 DEGs specific to LD-D vs. SD-D, 1794 DEGs specific to LD-M vs. SD-M, 244 DEGs specific to SD-D vs. SD-M, and 112 common DEGs between these comparison treatment groups, indicating that these identified DEGs may be involved in the response to the photoperiod. Fig. 1. Fig. 1 [101]Open in a new tab Venn diagrams of differentially expressed genes (DEGs) among different comparisons. There were 544 DEGs specific in LD-D vs LD-M, 3822 DEGs in LD-D_vs_SD-D, 1794 DEGs in LD-M_vs_SD-M, 244 DEGs in SD-D_vs_SD-M, and as well as 112 common DEGs between these different comparison treatment groups. LD-D: D32 Panicum miliaceum L. genotype under long-day treatment; LD-M: M51 Panicum miliaceum L. genotype under long-day treatment; SD-D: D32 Panicum miliaceum L. genotype under short-day treatment; SD-M: M51 Panicum miliaceum L. genotype under short-day treatment. Functional enrichment analysis of the identified DEGs Similarly, a 95% confidence interval (P < 0.05) was used as the standard to screen for significant GO terms and KEGG enrichment pathways. The identified DEGs of the LD-D vs. LD-M, SD-D vs. SD-M, LD-D vs. SD-D, and LD-M vs. SD-M groups were significantly enriched in 8, 21, 29, and 21 biological processes of GO terms, respectively. Furthermore, the DEGs between the M51 and D32 Panicum miliaceum L. genotypes involved in biological processes under the same light conditions were significantly different (Fig. [102]2). The DEGs between LD-D and LD-M were mainly enriched in cysteine-type peptidase activity, ubiquitin-like protein-specific protease activity, and UFM1 hydrolase activity under long sunlight, whereas the DEGs between SD-D and SD-M were primarily enriched in organonitrogen compound and protein metabolic processes and in an anchored component of the membrane. The enrichment pathways of the DEGs between the M51 and D32 Panicum miliaceum L. genotypes were similar to those under different light treatments (Table [103]1), and mainly included protein phosphorylation, phosphorus and phosphorus-containing compound metabolic processes, and phosphorylation. NEK (serine/threonine protein kinase) was involved in the phosphorylation of Panicum miliaceum L. protein, and the expression trends of NEK in photosensitive and photo-insensitive Panicum miliaceum L. were different, suggesting that the differential expression levels of NEK3 and NEK4 could distinguish the photoperiod sensitivity of the two types of Panicum miliaceum L.. Fig. 2. [104]Fig. 2 [105]Open in a new tab Gene ontology terms enrichment analysis of the identified DEGs between the LD-D and LD-M, and between the SD-D and SD-M. The horizontal axis represents the rich factor of DEGs, and the vertical axis represents the name of the item. The size of the dot represents the number of DEGs, and the color of the dot represents the enrichment significance. LD-D: D32 Panicum miliaceum L. genotype under long-day treatment; LD-M: M51 Panicum miliaceum L. genotype under long-day treatment; SD-D: D32 Panicum miliaceum L. genotype under short-day treatment; SD-M: M51 Panicum miliaceum L. genotype under short-day treatment. Table 1. The gene ontology terms of the DEGs between the LD-D and SD-D, and between the LD-M and SD-M. GO Term Gene Gene name Regulation LD-D vs SD-D LD-M vs SD-M Protein phosphorylation; C2845_PM03G21270 C2845_PM04G03900 NEK3 Up No Phosphorus metabolic process; C2845_PM07G29640 C2845_PM13G14270 NEK5 Up Up Phosphate-containing compound metabolic process; C2845_PM11G09590 NEK6 Up Up C2845_PM02G07000 NEK2 Up Up Phosphorylation; C2845_PM05G02360 C2845_PM06G33420 NEK2 Down Down C2845_PM05G36560 NEK4 Down Down C2845_PM06G20280 NEK4 No Down [106]Open in a new tab LD-D: D32 Panicum miliaceum L. genotype under long-day treatment; LD-M: M51 Panicum miliaceum L. genotype under long-day treatment; SD-D: D32 Panicum miliaceum L. genotype under short-day treatment; SD-M: M51 Panicum miliaceum L. genotype under short-day treatment. In addition, the DEGs between the LD-D and LD-M, SD-D and SD-M, LD-D and SD-M, and LD-M and SD-M groups were significantly enriched in 17, 9, 35, and 34 KEGG pathways, respectively. The metabolic pathways of the DEGs between the M51 and D32 Panicum miliaceum L. genotypes were significantly different under the same light treatments (Fig. [107]3). Under the long-day treatment, the common pathway was mainly phenylalanine biosynthesis, while the main metabolic pathways were photosynthesis, carotenoid biosynthesis, and glycerophospholipid metabolism; psbM, psal, and psaJ were involved in photosystem I and II reactions. The photoperiod gene PRR95 was significantly expressed under the long but not under the short photoperiod. The main metabolic pathways under short-day conditions were the circadian rhythm plant, diterpenoid, and monobactam biosynthesis. The photoperiod genes HD3A, FT, and PRR37 were significantly enriched in the circadian-plant metabolic pathway under short-day conditions but not under long-day conditions. The metabolic pathways involved in the DEGs between M51 and D32 were similar under different light treatments, including photosynthesis, photosynthesis-antenna protein, and valine, leucine, and isoleucine degradation. The expression of LHCA6 (the harvesting complex gene 6) was downregulated, and the cytochrome psbE gene was upregulated in M51. The expression levels of the amino acid transferase At3g08860, acetaldehyde dehydrogenase ALDH2B4, psbA, psaI, and psaC were significantly upregulated on D32, while the expression levels of IBR3 and ALDH6B2 dehydrogenase were significantly downregulated. Fig. 3. [108]Fig. 3 [109]Open in a new tab Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment pathways of the DEGs between the LD-D and LD-M, SD-D and SD-M, LD-D and SD-M, as well as LD-M and SD-M. The horizontal axis represents the rich factor of DEGs, and the vertical axis represents the name of the item. The size of the dot represents the number of DEGs, and the color of the dot represents the enrichment significance. LD-D: D32 Panicum miliaceum L. genotype under long-day treatment; LD-M: M51 Panicum miliaceum L. genotype under long-day treatment; SD-D: D32 Panicum miliaceum L. genotype under short-day treatment; SD-M: M51 Panicum miliaceum L. genotype under short-day treatment. [110]www.kegg.jp/kegg/kegg1.html. Annotation of metabolites A total of 822 metabolites were identified under different treatments (Fig. [111]4), including 538 and 284 metabolites in the positive and negative ion modes, respectively. Among these, 29.6% were lipids and lipid-like molecules, 16.8% were undefined, 8.9% were organic acids and derivatives, 8% were organoheterocyclic compounds, 7.3% were benzenoids, and 3.2% were nucleosides, nucleotides, and analogs. Fig. 4. [112]Fig. 4 [113]Open in a new tab The classification of the annotated metabolites under the different treatment groups. 822 metabolites were identified, which mainly belonging to lipids and lipid-like molecules, organic acids and derivatives, organoheterocyclic compounds, benzenoids, and nucleosides, nucleotides, and analogues. Statistical analysis of differential metabolites Differential metabolites in each comparison group were screened using the criteria VIP > 1 and P < 0.05. It was found between the LD-D and LD-M groups, there were 8 differential metabolites in the negative ion mode (such as L-malic and fumaric acids and protodioscin, Supplementary Table [114]S3), which were produced by synthesis and decomposition of organic acids and lipids; and 23 differential metabolites (such as isoreserpin, olivil 4'-O-glucoside, and ellipticine) were screened in the positive ion mode (Supplementary Table S3), generated by the synthesis and decomposition of ketones, alkaloids, and lipids compounds. Between the SD-M and SD-D groups, in the negative ion mode(Supplementary Table S3), 10 differential metabolites (such as protodioscin, L-malic acid, and flavone base + 4O, C-Hex-dHex) were found, which were involved in the synthesis and decomposition of lipids, organic acids, and flavonoids, while in the positive ion mode (Supplementary Table S3), 20 differential metabolites (such as cephaeline, cholestenone, and betaine) were involved in the biosynthesis and decomposition of alkaloids and lipids. Between the SD-D and LD-D groups, 21 differential metabolites (such as L-malic, fumaric, and glyceric acids) were detected in the negative ion mode (Supplementary Table [115]S3), which were generated by the synthesis and decomposition of organic acids, lipids, and oxides, and 28 differential metabolites (such as euphol, loperamide, and brucine) were detected in the positive ion mode (Supplementary Table [116]S3), which were mainly produced through the synthesis and decomposition of lipids, benzene ring compounds, and alkaloids. Furthermore, between the SD-M and LD-M groups, there were 14 differential metabolites (including L-malic and fumaric acids and D-glucose) in the negative ion mode (Supplementary Table [117]S3), which were generated through the synthesis and decomposition of organic acids and oxides, and 15 different metabolites (including loperamide, violanthin, and euphol) in the positive ion mode (Supplementary Table [118]S3), which were produced by the synthesis and decomposition of benzenoid compounds, ketones, and lipids. Taken together, the differential metabolites identified in negative mode (L-malic acid and fumaric acid) were mainly produced by the synthesis and decomposition of organic acids, whereas those identified in positive mode (brucine and loperamide) were primarily generated through the synthesis and decomposition of lipids. As shown in the Venn diagram (Fig. [119]5), there were 15, 6, 15, and 1 unique differential metabolites between the LD-M and LD-D groups, between the SD-M and SD-D groups, between the SD-D and LD-D groups, and between the SD-M and LD-M in the positive ion mode, respectively. In negative ion mode, there were three unique differential metabolites between the LD-D and LD-M groups, three between the SD-M and LD-M groups, 12 between the SD-D and LD-D groups, four between the SD-M and LD-M groups, and one common differential metabolite between each comparison group, L-malic acid. Fig. 5. [120]Fig. 5 [121]Open in a new tab Venn diagrams analysis of the differential metabolites between the different comparison groups in positive and negative ion modes. In the negative ion mode, one common differential metabolite was observed among each comparison groups, L-malic acid. LD-D: D32 Panicum miliaceum L. genotype under long-day treatment; LD-M: M51 Panicum miliaceum L. genotype under long-day treatment; SD-D: D32 Panicum miliaceum L. genotype under short-day treatment; SD-M: M51 Panicum miliaceum L. genotype under short-day treatment. KEGG analysis of the identified differential metabolites Differential metabolites identified between the different groups were subjected to KEGG enrichment analysis (Fig. [122]6). The main metabolic pathways of the differential metabolites in the LD-M and LD-D groups were associated with the citrate cycle (TCA cycle), pyruvate metabolism, and carbon fixation in prokaryotes. Differential metabolites between the SD-D and LD-D groups were primarily enriched in the TCA cycle, pyruvate metabolism, and carbon fixation in prokaryotes. Additionally, the differential metabolites in the SD-M and LD-M groups were closely related to carbon metabolism, TCA cycle, and pyruvate metabolism. Moreover, carbon metabolism, TCA cycle, and carbon fixation in photosynthetic organisms were observed in the differential metabolites between SD-M and SD-D. Interestingly, there was a common metabolic pathway among the different groups: the TCA cycle. Fig. 6. [123]Fig. 6 [124]Open in a new tab KEGG enrichment analysis of the differential metabolites between the LD-M and LD-D, SD-D and LD-D, SD-M and LD-M, as well as SD-M and SD-D. The horizontal axis represents the rich factor of DEGs, and the vertical axis represents the name of the item. The size of the dot represents the number of DEGs, and the color of the dot represents the enrichment significance. LD-D: D32 Panicum miliaceum L. genotype under long-day treatment; LD-M: M51 Panicum miliaceum L. genotype under long-day treatment; SD-D: D32 Panicum miliaceum L. genotype under short-day treatment; SD-M: M51 Panicum miliaceum L. genotype under short-day treatment. [125]www.kegg.jp/kegg/kegg1.html. Combined analysis of transcriptomic and metabolomic data A combined analysis was conducted based on differential metabolites and DEGs. The differential metabolites of the different photosensitive genotypes (M51 and D32) of Panicum miliaceum L. were significantly different in the negative mode (Fig. [126]7). The differential metabolites of the M51 Panicum miliaceum L. genotype were mainly d-glucose and those of the D32 Panicum miliaceum L. genotype were primarily L-malic and fumaric acids, which could assist in identifying the photoperiod sensitivity of Panicum miliaceum L.. Fig. 7. [127]Fig. 7 [128]Open in a new tab Differential expression analysis of the differential metabolites and DEGs of the different photosensitive Panicum miliaceum L. in the negative mode. The red columns represent differential metabolites, and the green columns represent mRNA. The horizontal axis means metabolite-mRNA, and the vertical axis means log2 fold change. LD-D: D32 Panicum miliaceum L. genotype under long-day treatment; LD-M: M51 Panicum miliaceum L. genotype under long-day treatment; SD-D: D32 Panicum miliaceum L. genotype under short-day treatment; SD-M: M51 Panicum miliaceum L. genotype under short-day treatment. The number of common KEGG pathways for both differential metabolites and DEGs in the comparison groups of LD-D vs SD-D and LD-M vs SD-M were 16 (Table [129]2), with two common metabolic pathways in the LD-D vs. SD-D comparison group, which included carbon fixation in photosynthetic organisms as well as alanine, aspartate, and glutamate metabolism. There were five common metabolic pathways between the LD-M and SD-M groups: alanine, aspartate, and glutamate metabolism; pyruvate metabolism; glycolysis/gluconeogenesis; the pentose phosphate pathway; and butanoate metabolism. In summary, the DEGs and differential metabolites of the photosensitive M51 Panicum miliaceum L. and photoinsensitive D32 Panicum miliaceum L. genotypes were both significantly enriched in alanine, aspartate, and glutamate metabolism; and the DEGs and differential metabolites in the photosensitive (M51) Panicum miliaceum L. were enriched in four unique pathways, containing pyruvate metabolism, glycolysis/gluconeogenesis, pentose phosphate pathway, and butanoate metabolism; as well as the DEGs and differential metabolites in the photoinsensitive (D32) Panicum miliaceum L. were associated with the unique pathway of carbon fixation in photosynthetic organisms. Table 2. KEGG pathway analysis of differential metabolites and differentially expressed genes (DEGs) in different comparison groups. Treatment Map_ID Map_Name p.value LD-D vs SD-D meta ko00020 Citrate cycle (TCA cycle) 0.0003 ko00620 Pyruvate metabolism 0.0007 ko00190 Oxidative phosphorylation 0.0215 ko00220 Arginine biosynthesis 0.0308 ko00710 Carbon fixation in photosynthetic organisms 0.0308 ko00250 Alanine, aspartate and glutamate metabolism 0.0374 ko00480 Glutathione metabolism 0.0504 ko00650 Butanoate metabolism 0.0556 ko00944 Flavone and flavonol biosynthesis 0.0671 ko00760 Nicotinate and nicotinamide metabolism 0.0722 ko00360 Phenylalanine metabolism 0.0786 ko00630 Glyoxylate and dicarboxylate metabolism 0.0811 ko00240 Pyrimidine metabolism 0.0849 ko00941 Flavonoid biosynthesis 0.0961 ko00350 Tyrosine metabolism 0.1011 ko02010 ABC transporters 0.1714 LD-D vs SD-D trans ko00250 Alanine, aspartate and glutamate metabolism 0.0009 ko00240 Pyrimidine metabolism 0.0021 ko00941 Flavonoid biosynthesis 0.0022 ko00630 Glyoxylate and dicarboxylate metabolism 0.0036 ko00360 Phenylalanine metabolism 0.0055 ko00710 Carbon fixation in photosynthetic organisms 0.0070 ko00650 Butanoate metabolism 0.0134 ko00350 Tyrosine metabolism 0.0774 ko00944 Flavone and flavonol biosynthesis 0.1175 ko00620 Pyruvate metabolism 0.2548 ko00760 Nicotinate and nicotinamide metabolism 0.2640 ko02010 ABC transporters 0.6442 ko00480 Glutathione metabolism 0.7723 ko00220 Arginine biosynthesis 0.8139 ko00020 Citrate cycle (TCA cycle) 0.9666 ko00190 Oxidative phosphorylation 1.0000 LD-M vs SD-M meta ko00020 Citrate cycle (TCA cycle) 0.0002 ko00620 Pyruvate metabolism 0.0004 ko00190 Oxidative phosphorylation 0.0162 ko00220 Arginine biosynthesis 0.0232 ko00710 Carbon fixation in photosynthetic organisms 0.0232 ko00250 Alanine, aspartate and glutamate metabolism 0.0282 ko00010 Glycolysis / Gluconeogenesis 0.0311 ko00030 Pentose phosphate pathway 0.0351 ko00650 Butanoate metabolism 0.0420 ko00944 Flavone and flavonol biosynthesis 0.0508 ko00564 Glycerophospholipid metabolism 0.0518 ko00760 Nicotinate and nicotinamide metabolism 0.0547 ko00360 Phenylalanine metabolism 0.0595 ko00630 Glyoxylate and dicarboxylate metabolism 0.0614 ko00941 Flavonoid biosynthesis 0.0730 ko00350 Tyrosine metabolism 0.0768 LD-M vs SD-M trans ko00250 Alanine, aspartate and glutamate metabolism 0.0028 ko00941 Flavonoid biosynthesis 0.0046 ko00620 Pyruvate metabolism 0.0238 ko00030 Pentose phosphate pathway 0.0279 ko00010 Glycolysis / Gluconeogenesis 0.0329 ko00650 Butanoate metabolism 0.0415 ko00360 Phenylalanine metabolism 0.0452 ko00630 Glyoxylate and dicarboxylate metabolism 0.0826 ko00710 Carbon fixation in photosynthetic organisms 0.0947 ko00944 Flavone and flavonol biosynthesis 0.2323 ko00564 Glycerophospholipid metabolism 0.2802 ko00350 Tyrosine metabolism 0.2909 ko00760 Nicotinate and nicotinamide metabolism 0.3961 ko00020 Citrate cycle (TCA cycle) 0.8584 ko00220 Arginine biosynthesis 0.8895 ko00190 Oxidative phosphorylation 1.0000 [130]Open in a new tab LD-D: D32 Panicum miliaceum L. genotype under long-day treatment; LD-M: M51 Panicum miliaceum L. genotype under long-day treatment; SD-D: D32 Panicum miliaceum L. genotype under short-day treatment; SD-M: M51 Panicum miliaceum L. genotype under short-day treatment. [131]www.kegg.jp/kegg/kegg1.html. Correlation analysis of the DEGs and differential metabolites in the pathway enriched by photosensitive M51 and photoinsensitive D32 showed that these metabolites were involved in the synthesis and decomposition of organic acids in the negative ion mode, and the main metabolites were L-malic and fumaric acids. In the positive mode, the main metabolites were brucine and loperamide, which are involved in lipid synthesis and decomposition. Additionally, LHCA6 (light-harvesting complex I chlorophyll a/b binding protein) was down-regulated in the M51 Panicum miliaceum L. genotype, while the cytochrome psbE gene was up-regulated; as well as in the D32 Panicum miliaceum L. genotype, the expression levels of amino acid transferase At3g08860, acetaldehyde dehydrogenase ALDH2B4, psbA, psaI, and psaC genes were significantly up-regulated, whereas the expression levels of IBR3 and ALDH6B2 dehydrogenase were significantly down-regulated. To screen the differentially upregulated and downregulated genes in the LD-D vs. SD-D and LD-M vs. SD-M comparison groups, the DEGs and differential metabolites that could accurately identify different photoperiod-sensitive Panicum miliaceum L. were mined and combined with the metabolomic data. A total of 18 genes were found to be significantly involved in the differential expression between the photoinsensitive D32 and photosensitive M51 Panicum miliaceum L. genotypes (Fig. [132]8), and candidate gene36259 (name: C2845_PM11G01290), which was annotated as V-type H+-transporting ATPase subunit B (ATP6B), significantly increased the levels of L-malic and fumaric acids. Fig. 8. [133]Fig. 8 [134]Open in a new tab Correlation analysis of DEGs and differential metabolites in different comparison groups. The main metabolites in the negative ion mode were L-malic and fumaric acids; while the main metabolites in the positive mode were brucine and loperamide. A total of 18 genes were significantly involved in the differential expression between D32 and M51 Panicum miliaceum L. genotypes, and the candidate gene36259 (name: C2845_PM11G01290) was found. LD-D: D32 Panicum miliaceum L. genotype under long-day treatment; LD-M: M51 Panicum miliaceum L. genotype under long-day treatment; SD-D: D32 Panicum miliaceum L. genotype under short-day treatment; SD-M: M51 Panicum miliaceum L. genotype under short-day treatment. qRT-PCR validation To verify the reliability of the transcriptomic data, 10 genes involved in photosynthesis, carbon fixation in photosynthetic organisms, and glycerophospholipid, pyruvate, alanine, aspartate, and glutamate metabolism were randomly selected for qRT-PCR analysis. The expression trends of the 10 genes in the different combinations were consistent with those of the transcriptomic data (Fig. [135]9A), with a consistency rate of 100%, indicating the reliability of the transcriptome sequencing results. Fig. 9. [136]Fig. 9 [137]Open in a new tab qRT-PCR verification of the related genes. (A) The expression levels of the 10 genes involved in photosynthesis, carbon fixation in photosynthetic organisms, and glycerophospholipid, pyruvate, alanine, aspartate, and glutamate metabolism in the different groups verified by qRT-PCR. (B) The expression of the candidate gene (C2845_PM11G01290) in 24 h time-period determined in two kinds of Panicum miliaceum L. (M51 and D32). LD-D: D32 Panicum miliaceum L. genotype under long-day treatment; LD-M: M51 Panicum miliaceum L. genotype under long-day treatment; SD-D: D32 Panicum miliaceum L. genotype under short-day treatment; SD-M: M51 Panicum miliaceum L. genotype under short-day treatment. In addition, the expression of the candidate gene (C2845_PM11G01290) over a 24 h time-period was determined in two types of Panicum miliaceum L.. It was found that with an increase in light time, the expression of C2845_PM11G01290 first increased, decreased, and then increased again, and the degree of photosensitivity in Panicum miliaceum L. (M51) was more significant than that of photo-insensitive Panicum miliaceum L. (D32) (P < 0.05, Fig. [138]9B). In the 24 h time-period expression analysis, the peak expression levels primarily occur between 18 and 21 h. Furthermore, the gene expression level of D32 is significantly lower than that of M51, suggesting that the impact on plants is mainly concentrated within the 18 h-21 h time frame. Although the ATP6B gene expression level of D32 (light-insensitive) is higher than that of M51(light-sensitive) at 0 h, 3 h, and 12 h, its influence on the growth and development of different materials is relatively minimal. Throughout the entire 24 h candidate gene expression trend, the magnitude of change in D32 is significantly lower than that of M51, which is consistent with its phenotypic changes (Table [139]3). Table 3. Photoperiodic relative sensitivity (PRS) of D32 and M51 agronomic traits. Agronomic traits Aboveground fresh weight (g) Plant height (cm) Flag leaf area (cm^2) Panicle length (cm) 16 h 8 h PRS 16 h 8 h PRS 16 h 8 h PRS 16 h 8 h PRS D32 20.28 23.95 -0.181 94.8 81.2 0.143 33.67 26.67 0.208 26.7 20.4 0.236 M51 70.12 8.10 0.885 120.7 34.8 0.712 47.88 9.15 0.809 27.2 10.6 0.610 [140]Open in a new tab PRS = (LD-SD) / LD, LD represents the value of a certain trait under long day of 16 h, SD represents the value of a certain trait under short day of 8 h. Discussion Although the mechanisms of photoperiod regulation differ among crops, the key genes involved in photoperiod regulation in different crops may be conserved, providing a theoretical basis for the study of Panicum miliaceum L. photoperiodic responses. Likewise, heading date is not only the most intuitive phenotypic trait to elucidate photoperiod but is also a complex quantitative trait that is regulated by many major and minor genes. It is worth mentioning that the photoperiod sensitivity of Panicum miliaceum L. was higher than that of other gramineous crops. In this study, we used a photoperiod sensitivity evaluation system for Panicum miliaceum L.^[141]52, photoperiod sensitive (M51), and photoperiod-insensitive (D32) genotypes. We identified three photoperiod candidate genes (HD3A, FT, and PRR37) in the circadian-plant metabolic pathway by transcriptome sequencing, which were significantly expressed under short-day conditions but not in response to long-day conditions. These results are consistent with previous reports showing that the metabolic pathways of photoperiods via DEGs mainly focus on circadian rhythms, photosynthesis, and antenna proteins^[142]53,[143]54. Komiya et al.^[144]12 showed that Hd3a and RFT1 could act as flower activators under short-day conditions, and that RFT1 expression is partly regulated by chromatin modification, thus promoting flower transition. Another study indicated that the circadian pathway was the most significantly enriched pathway in Glycine max (Huaxia No. 3) under short-day conditions and that the FT homolog gene may be involved in Glycine max flowering under short-day conditions^[145]55. The Hd1-DTH8-Ghd7-PRR37 module is a photoperiodic flowering core module that improves latitude adaptation by mediating the multiday length perception process in Oryza sativa L.. PRR37, ELF3-1 and Ehd1 mutants can alter the day-length sensing process and exhibit longer main spike lengths and more grains per spike, which implies that their sites have excellent latitude adaptation and yield improvement potential in Oryza sativa L. breeding^[146]56. Taken together, we speculate that HD3A, FT, and PRR37 may be crucial photoperiod-related markers in Panicum miliaceum L.. However, their specific roles in Panicum miliaceum L. under different photoperiods need to be further elucidated. In addition, metabolites are the ultimate carriers of physiological and biochemical activities as well as the final products of gene expression. They adapt to their environment through complex intracellular regulatory mechanisms, ultimately affecting phenotypic traits. Our results showed that the different comparison groups exhibited a common TCA cycle metabolic pathway: L-malic and fumaric acids in the negative ion mode and brucine and loperamide in the positive ion mode. The TCA cycle is essential for cellular energy metabolism and the carbon skeleton supply, particularly for sugar, fat, and protein catabolism^[147]57. In addition to providing energy, the TCA cycle participates in many cellular metabolic pathways, including photosynthesis, photorespiration, abiotic stress, biological clocks, hormone signaling, and glycolysis^[148]58. A previous study demonstrated that the levels of most identified metabolites (primarily amino acids, sugars, and TCA cycle intermediates) increased after 2 h and peaked after 73 h, indicating that the TCA cycle is correlated with metabolic changes in Arabidopsis thaliana plants growing under different light conditions^[149]59. L-malic acid, a TCA cycle intermediate secreted from the roots of Arabidopsis, selectively signals and recruits the beneficial rhizobacterium Bacillus subtilis FB17 in a dose-dependent manner^[150]60. Fumaric acid is a component of the tricarboxylic acid cycle and its concentration increases with plant age and light intensity in Arabidopsis leaves^[151]61. Brucine and loperamide are natural plant alkaloids that have been reported to have cytotoxic and anti-proliferative activities^[152]62,[153]63, however, their effects on plant growth have not been determined. Combined with our results, it can be inferred that L-malic and fumaric acid, brucine, and loperamide may be pivotal metabolites for the photosensitive (M51) and photo-insensitive (D32) genotypes of Panicum miliaceum L. under different light conditions and that the TCA cycle may be the key pathway for Panicum miliaceum L. under different photoperiods. However, the specific roles and mechanisms of the key metabolites and pathways identified in Panicum miliaceum L. under different photoperiods require further investigation. In addition, combined transcriptomic and metabolomic analyses were performed to provide a better understanding of the key regulatory metabolic pathways, genes, and metabolites, thus revealing the molecular response mechanisms of Panicum miliaceum L. to different photoperiods. The identified DEGs and metabolites of photosensitive M51 and photo-insensitive D32 were significantly enriched in the alanine, aspartate, and glutamate metabolism pathways. Zhou^[154]64 pointed out that the coverage of KEGG enrichment analysis results between conserved rhythm transcripts and metabolites was limited in Oryza sativa L., and only in amino acid and nucleotide metabolic pathways, such as alanine and glycine. Consistent with this, we also found that these were mainly concentrated in non-essential amino acid pathways that were strongly affected by the photoperiod. Some pathways did not intersect at the transcriptome and metabolome levels and were only enriched at the transcriptome or metabolome level; however, this does not mean that these pathways lack the regulation of photoperiodic transcriptional genes. Each enrichment pathway was determined based on the significance of the proportion of genes associated with the photoperiod. Not all pathways are significantly enriched, and each transition stage between transcripts, proteins, and metabolites is delayed^[155]65,[156]66. Additionally, a candidate gene (gene36259, name: C2845_PM11G01290), which was annotated to ATP6B, was screened for photoperiod, and its expression level was found to be involved in ATP binding and ATP metabolic processes and affected the levels of L-malic and fumaric acids. We also measured the 24 h rhythmic expression of C2845_PM11G01290 (ATP6B) and found that with an increase in light exposure, its expression first increased, then decreased, and then increased again, and the degree of photosensitivity in Panicum miliaceum L. (M51) was more significant than that in photo-insensitive Panicum miliaceum L. (D32). Plants rely on light energy to produce ATP and reducing agents, as well as a supply of nutrients (inorganic C, N, and S compounds) to successfully produce biomass^[157]67. Extracellular ATP is now recognized as an important signaling factor in plant growth and defense responses to environmental stimuli^[158]68. A previous study showed that Cannabis sativa L. adapted to Pb stress by: accelerating ATP metabolism; enhancing respiration, light absorption and light energy transfer; regulate stomatal development and closure; promoting intercellular transport; and increasing the transmembrane transport of ATP in chloroplasts, which providing important reference protein and gene information for future molecular studies of Pb resistance and accumulation in Cannabis sativa L.^[159]69. These reports, together with our findings, suggest that amino acid and nucleotide metabolic pathways, such as alanine, glycine, and ATPase transport, may participate in different Panicum miliaceum L. under different photoperiods and that C2845_PM11G01290 (ATP6B) may be closely related to the circadian rhythmic pathway of Panicum miliaceum L.. However, the specific mechanisms of the identified pathways and candidate genes (C2845_PM11G01290, ATP6B) in different Panicum miliaceum L. under different photoperiods warrant further exploration. Conclusions The three main components of the photoperiod pathway in higher plants are photoreceptors, circadian clocks, and circadian regulatory genes. We found that the rhythmic genes of Panicum miliaceum L., HD3A, FT, and PRR37 were significantly enriched in the circadian-plant metabolic pathway under short-day conditions. However, the photoperiodic candidate genes and metabolites of different photosensitive Panicum miliaceum L. differed significantly in response to changes in photoperiods. Nonetheless, the main metabolite of the photosensitive genotypes included D-glucose, whereas in the photo-insensitive genotype, L-malic and fumaric acids were potential metabolites. In addition, a significantly different photoperiod candidate gene (C2845_PM11G0129, ATP6B) was identified by combined transcriptomics and metabolomics analyses, and this gene could be used to identify different photosensitive genotypes using molecular characterization. Supplementary Information [160]Supplementary Figure 1.^ (301.2KB, png) [161]Supplementary Table 1.^ (11.2KB, xlsx) [162]Supplementary Table 2.^ (11.8KB, xlsx) [163]Supplementary Table 3.^ (71.5KB, xlsx) Acknowledgements