Comparative transcriptomics reveals key genes contributing to the differences in drought tolerance among three cultivars of foxtail millet ( Setaria italica )

Foxtail millet ( Setaria italica L.) is an important food crop with strong drought stress tolerance. The response of foxtail millet to drought stress is a complex regulatory network. It is of great guiding significance for agricultural production to continuously explore candidate genes of foxtail millet drought resistance and reveal the molecular regulatory mechanism and metabolic pathway of foxtail millet drought tolerance. This study investigated three different cultivars of foxtail millet with different drought resistance. Drought stress reduced the water and chlorophyll content, and increased the Malondi-aldehyde (MDA) level and soluble protein of foxtail millet leaves. From strong to weak, the drought resistance order was Jigu39, Jingu21, and Longgu16. The transcriptome analysis of these three cultivars was carried out. 2954, 1531, and 2344 deferentially expressed genes under drought stress were identified in Jigu39, Jingu21, and Longgu16. The GO and KEGG pathway analysis identified DEGs significantly enriched in photosynthesis, chlorophyll metabolism, amino acid metabolism, and carbohydrate metabolism in all three cultivars. In addition, we identified 46 genes whose trends of transcription change were consistent with the drought resistance trends among three cultivars of foxtail millet. Among them, 32 genes were first identified related to drought response in foxtail millet, which was worth further research.


Introduction
Global climate conditions are deteriorating due to increased human activities such as urbanization, rapid industrialization, and overuse of non-biodegradable commodities, thereby exposing plants to harsh and extreme climatic conditions, suffering abiotic stresses such as drought, salinity, and cold, adversely affecting plant morphology, growth and development, and molecular processes (Chaudhry and Sidhu 2022). Drought is considered to be one of the major abiotic stresses limiting crop production worldwide because it reduces the supply of water in the crop and seriously affects the growth rate and biomass of plants, thereby affecting global food security (Delauney and Verma 1993;Mishra and Cherkauer 2010;Anjum et al. 2011).
Foxtail millet (Setaria italica L.) is an important food crop in arid and semi-arid regions (He et al. 2015). As a C4 plant, foxtail millet can maintain a high photosynthetic rate under drought conditions than other traditional C3 crops such as wheat and rice (Pant et al. 2016). Foxtail millet has a small genome (515 Mb), less repeated DNA, and a short life cycle, making it an ideal model system for the C4 plant study (Li and Brutnell 2011;Bennetzen et al. 2012;Zhang et al. 2012;Qin et al. 2020). Foxtail millet also requires less soil moisture, and its transpiration efficiency and water use efficiency (WUE) are much higher than other species, such as sorghum, maize, and wheat (Tang et al. 2017a). Research focusing on foxtail millet will help understand the mechanism of plant stress tolerance.
The adaptation of foxtail millet to drought stress includes biochemical and physiological changes. Foxtail millet weakens or offsets the adverse effects of drought by reducing its Communicated by Xingchun Wang. net photosynthetic rate, regulating stomatal closure, synthesizing osmotic regulators to balance cell osmotic pressure, and synthesizing various anti-oxidants (Lawlor and Cornic 2002). However, the response of foxtail millet to drought stress was a complex regulatory network. Previous studies have shown that the physiological and biochemical levels of foxtail millet responses to drought stress vary among different cultivars, and studies have found that some physiological traits and indicators related to drought resistance, such as the ratio of ascorbate peroxidase to superoxide dismutase, are key factors in alleviating drought stress (Bhatt et al. 2011;Bartwal et al. 2016;Mukami et al. 2019). Therefore, we can screen out foxtail millet cultivars with drought-resistant potential and further explore drought-resistant genes through transcriptomics.
At present, the deep transcriptional studies on the drought tolerance in foxtail millet mainly focus on the different development stages and tissues of one or two foxtail millet cultivars with different drought tolerance (Tang et al. 2017b;Xu et al. 2019;Qin et al. 2020;Zhang et al. 2022). In our previous study, we detected the transcription levels of 52 SNARE (soluble N-ethylmaleimide-sensitive factor attachment protein receptor) genes in Jingu21 and Longgu16 under drought stress ) and found that NPSN11 (Novel Plant SNARE 11, Seita. 4G280200) was only downregulated in Jingu21, which was different from the report of Qi et al. (2013). The above results indicated that transcriptional regulation of some drought-induced genes was different in different foxtail millet cultivars. Therefore, In this study, we further expanded the scope of the comparative transcriptome and compared the drought resistance of three foxtail millet cultivars (Jigu39, Jinu21, and Longgu16). Jigu39, Jingu21, and Longgu16 were selected and bred by the Chinese Academy of Agricultural Sciences in Hebei, Shanxi, and Gansu provinces. The results showed that they had different drought tolerance. Jigu39 had the strongest drought resistance, followed by Jingu21 and Longgu16. Then, we used these three cultivars as materials and performed transcriptome analysis under drought stress by stopping watering. Differentially expressed genes (DEGs) related to drought resistance were isolated and identified. Through gene clustering and related metabolic pathway analysis, the mechanisms of foxtail millet response to drought were systematically expounded from the physiological, biochemical, and transcriptome levels. In total, 46 genes were identified whose trends of transcription change were consistent with the drought resistance trends among three cultivars. Furthermore, 32 genes were first reported as drought-related in foxtail millet. Our research provided useful information and a theoretical basis for further exploring the drought resistance mechanism of foxtail millet.

Plant materials, growth conditions and stress treatment
In this study, three different cultivars of foxtail millet (Longgu16, Jingu21, and Jigu39) were used as experimental materials. All three seeds were collected from Prof. Lizhen Zhang (School of Life Science of Shanxi University) for free and used for research only. Three cultivars of foxtail millet were planted in a mixture (nutrient soil: vermiculite = 1:1) to grow in a greenhouse for 4 days, and growth-consistent seedlings were moved to the plastic pot. Seedlings were grown for 14 days in a greenhouse at 16 h light/8 h darkness, 23-26 °C, 50,000 Lux light, and 60-70% relative humidity, and then treated with drought stress by stopping watering for another 14 days. The leaves were collected from the plant with normal watering or drought stress, immediately frozen in liquid nitrogen and stored at − 80 °C for RNA sequencing and qRT-PCR analysis. There were three biological replicates and three technical replicates per experiment to ensure the reliability of experimental data.

Physiology assay
The water loss rate was measured as a previous report (Aguilar et al. 2000). The detached fresh leaves of 14-day-old seedlings were dried at room temperature for 48 h (15-25% relative humidity), and the weight of the leaves was measured. The weight of detached leaves calculated the water loss rate at each marked time. The weight of 0 h was fresh weight (FW), and the weight of leaves at the corresponding time was recorded as dry weight (DW). The leaf water loss rate at the corresponding time was calculated according to the following formula.

RNA Isolation and Library Preparation
The transcriptome sequencing and preliminary analysis were carried out by OE Biotech (Shanghai, China). Total RNA was extracted by the TRIzol method (Ambion, Thermo Scientific Inc., USA). After removing the DNA in total RNA using DNase, the mRNA was enriched by magnetic beads and interrupted by the reagent. The double-stranded cDNA was gradually synthesized as a template, and the end of the double-stranded cDNA was filled. A tail was added, and Leaf water loss rate (%) = (FW − DW)∕FW × 100.
the adapter was added. Finally, the library was constructed by PCR amplification, and the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) was used for quality inspection.

RNA sequencing and differentially expressed genes analysis
The constructed libraries were sequenced using the Illumina HiSeq X Ten platform, processed to obtain clean readings and mapped via HISAT2 (Kim et al. 2015) to the Setaria italica v2.2 transcripts in the JGI transcriptome reference database (https:// phyto zome-next. jgi. doe. gov/).
After the genes were obtained, their fragments per kilobase of transcript per millions mapped reads (FPKM) values were calculated, readings were taken using HTSeq count, and differential expression analysis was performed using the DESeq R software package to obtain DEGs. P value < 0.05 and foldchange > 2 or foldchange < 0.5 were used as the thresholds for significant difference expression. Hierarchical clustering analysis, GO enrichment and KEGG pathway enrichment analysis were performed on DEGs.

qPCR validation in the illumina data
The accuracy of the RNA-seq data needs to be verified by quantitative real-time PCR (qRT-PCR). The total RNA of foxtail millet leaves was extracted using TransZol™ UP Plus RNA Kit (TransGen Biotech, Beijing, China) and reversely transcribed using EasyScript® One-Step gDNA Removal and cDNA Synthesis SuperMix Kit (TransGen Biotech, Beijing, China). QRT-PCR was performed using the Perfect-Start Green qPCR SuperMix kit (TransGen Biotech, Beijing, China), and the results were analyzed by the 2 −ΔΔCT method (Livak and Schmittgen 2001). SiAct2 (Seita. 8G043100) is used as an internal reference. All primers used in the experiment were synthesized by Sangon Biotech, Shanghai (Supplement Table 1). Three biological and technical replicates were performed for each sample in the qRT-PCR to reduce the impact of biological and technical differences.

Statistical analysis
Each experiment has at least three completely randomized block designs with biological repetitions to ensure the accuracy of the experimental results. Student's t-tests were performed in SPSS Statistics 22.0 to evaluate the data of physiological analysis (SPSS Inc., Chicago, IL, USA). Each data value is expressed as mean ± standard error (SE). P < 0.05 indicated that there was a statistically significant difference.

Longgu16 shows better drought stress tolerance than the other two cultivars
After germination, all cultivars showed similar growth and development pattern during the time range of the experiment, suggesting there was no significant development difference among the three cultivars (Supplement Fig. 1). The phenotype response of Jigu39, Jingu21 and Longgu16 to drought stress are shown in Fig. 1a. After drought treatment, the leaves of all three cultivars became yellow and wrinkled to a certain extent. Longgu16 showed the most serious of wrinkles. After plants are subjected to drought stress, they will reduce water loss by closing stomata. The water loss in detached leaves can reflect the drought resistance of plants (Carmody et al. 2016). Here, the water loss rate of detached leaves within 48 h was detected, and the water loss rate was calculated. The results showed that the leaf water loss rate of Longgu16 was significantly higher than that of Jigu39 and Jingu21 at the same time, and that of Jigu39 was the slowest (Fig. 1b), supporting their phenotypes shown in Fig. 1a.
MDA is an indicator of membrane lipid peroxidation under stress. The change of MDA content can reflect the damage degree of membrane lipid under drought stress and indirectly reflect the drought resistance of crops (Hu et al. 2013). As shown in Fig. 1c, the MDA contents of the three cultivars of seedlings increased, indicating that their cell membrane was damaged under drought stress (Fig. 1c). Chlorophyll content can indicate the utilization rate of light energy in crops. When crops are under drought stress, the enzyme synthesis of chlorophyll is hindered (Impa et al. 2012), resulting in a decrease of chlorophyll content (Meher et al. 2018). This study found that the chlorophyll contents of Jigu39, Jingu21 and Longgu16 leaves decreased significantly after drought stress (Fig. 1d). Again, Jigu39 showed a higher level of chlorophyll content, even after drought stress, which might be related to its better drought tolerance.
Soluble proteins are important osmoregulation substances. The accumulation of soluble proteins improves the water retention ability of cells and protects the vital substances and biofilm of cells, so they are often used as plant stress resistance indicators (Smirnoff 1993;Sharma and Dubey 2005). As shown in Fig. 1e, it was found that after drought stress, the soluble protein contents in Jigu39 leaves were significantly increased, while the soluble protein contents in Longgu16 leaves had no significant change (Fig. 1e), which also supported the results above.
In conclusion, Jigu39 showed the highest drought stress tolerance, while Longgu16 was the weakest among the three cultivars.

RNA-seq analysis
To identify the key genes of foxtail millet in response to drought stress, leaves from three foxtail millet cultivars were taken for RNA-seq analysis with or without stop watering treatment. Each treatment had three biological replicates, and 18 cDNA libraries were constructed (marker as Jigu39 and Jigu39X, Jingu21, and Jingu21X, Longgu16 and Long-gu16X, with X meaning sample under drought stress). A total of 120.99 G of clean data was acquired from the referenced transcriptome sequencing of 18 samples. The effective data volume of each sample was distributed in 6.45-7.09 G, the Q30 base was distributed in 94.92-95.33%, and the average GC content was 53.87%. Use HISAT2 (Kim et al. 2015) to align clean reads with the specified reference genome. The sequencing and alignment results are summarized in Supplement Table 2. The reads were mapped to the reference genome and showed high alignment rates. Alignment rates ranged from 96.85 to 98.05%, and unique mapped reads ranged from 94.33% to 95.81%.

Differential expression analysis of transcripts
Htseq-count software (Anders et al. 2015) was used to obtain the number of reads on protein-coding genes. The cufflinks software (Trapnell et al. 2010) was used to calculate these genes' FPKM value (Roberts et al. 2011). According to the gene expression of all samples, principal component analysis (PCA) and cluster analysis were performed by R software to evaluate the repeatability and correlation of samples. It can be seen from Supplement Fig. 2 that the repeatability of three biological repeats in each treatment is satisfactory, and the difference between treatments is obvious, indicating that the current RNA-Seq data is reliable for subsequent analysis.
The three repeated data of each treatment were summarized and compared as Jigu39X with Jigu39 (Jigu39X/ tein. Asterisk indicates that the gene expression after drought stress treatment has a significant difference compared with the control (**P < 0.01) 1 3 Jigu39), Jingu21X with Jingu21 (Jingu21X/Jingu21), and Longgu16X with Longgu16 (Longgu16X/Longgu16) to further identify differentially expressed genes (DEGs) in response to drought stress, P value < 0.05 and |log 2 Fold-Change|> 2 was set as the threshold for significantly differential expression. The number of DEGs in different cultivars was summarized in Fig. 2a. Before and after drought stress, 24966 DEGs were identified in three cultivars of foxtail millet. After drought stress, there were 4581 up-regulated DEGs and 4896 down-regulated DEGs in Jigu39, 3442 upregulated DEGs and 3160 down-regulated DEGs in Jingu21, and 4068 up-regulated DEGs and 4819 down-regulated DEGs in Longgu16. The overall distribution of differentially expressed genes is shown in the volcano map (Fig. 3).
In order to further understand the regulatory mechanism of response to drought stress in foxtail millet, we compared DEGs among Jigu39X/Jigu39, Jingu21X/Jingu21, and Longgu16X/Longgu16. A total of 2717 (18.69%) overlapping genes were found in the three cultivars. 94.81% of overlapping genes were up-regulated and down-regulated in the same trend among the three cultivars, and these overlapping DEGs were more likely to participate in the metabolic Differentially expressed volcano map. Gray is the gene with a non-significant difference, red and green are the genes with a significant difference, the horizontal axis is log 2 FoldChange, and the vertical axis is − log 10 PValue. (Color figure online) pathways that affected the growth and development of foxtail millet under drought stress. In addition, there were 2954, 1531 and 2344 cultivar-specific DEGs in Jigu39, Jingu21, and Longgu16, respectively. Among these cultivar-specific DEGs, there were 1655, 809, and 1310 up-regulated genes and 1299, 722, and 1034 down-regulated genes in three cultivars (Fig. 2a). The expression specificity of many differential genes in the three cultivars of foxtail millet showed different regulation modes of drought stress, which may be one of the reasons why the tested cultivars were sensitive or tolerant to drought stress.

Gene ontology (GO) analysis of DEGs
After obtaining DEGs, GO enrichment analysis of DEGs was performed, and their functions were described combined with GO annotation results. Based on GO annotation, in Jigu39X/Jigu39 comparison group, drought stress mainly affects two biological processes: cellular process and metabolic process. The cellular component mainly affects the cell, cell part, membrane, complex, membrane part, organelle, and organelle part. For molecular function, drought stress mainly affects binding and catalytic activity, structural molecule activity, and transporter activity (Supplement Fig. 3). The main changes in biological process, cellular component, and molecular function of Jingu21X/Jingu21 and Longgu16X/Longgu16 are similar to that of Jigu39X/ Jigu39 (Supplement Fig. 3). In addition, the results for significantly enriched Top30 GO terms (P adj -value < 0.05) are presented in Supplement Table 3. In addition, to further understand the regulatory network of drought stress response of three foxtail millet cultivars, we also performed GO enrichment analysis on the overlapping and specific DEGs in this study (Supplement Fig. 4).

Kyoto encyclopedia of genes and genomes (KEGG) metabolic pathway analysis of DEGs
Based on KEGG pathway enrichment analysis, 2052 (Jigu39X/Jigu39), 1227 (Jingu21X/Jingu21) and 1838 (Longgu16X/Longgu16) DEGs were assigned to 372, 355 and 377 pathways (Supplement Table 4). These three comparisons were different in the metabolic pathways of DEGs. Jigu39, with strong drought resistance, was significantly enriched in photosynthesis antenna proteins, biosynthesis of ansamycins, ribosome, photosynthesis, ribosome biogenesis in eukaryotes, glyoxylate and dicarboxylate metabolism, and carbon fixation in photosynthetic organism pathways (Fig. 4a). DEGs in Jingu21 are mainly enriched in plant hormone signal transduction, MAPK signalling pathway-plant, type I diabetes mellitus and other pathways (Fig. 4b). DEGs' main enriched metabolic pathways in Longgu16 included biosynthesis of ansamycins, photosynthesis-antenna proteins, glyoxylate and dicarboxylate metabolism, photosynthesis, porphyrin and chlorophyll metabolism, and carbon fixation in the photosynthetic organisms (Fig. 4c).

Drought stress-responsive transcriptional factors
Transcription factors (TFs) play a central role in regulating biological and abiotic stress responses and various biological processes by regulating plant functional genes (Mickelbart et al. 2015;Baillo et al. 2019). This study identified 122 TFs from 26 TF families in the detected DEGs (Supplement Table 5, Supplement Table 6). These TFs identified in the leaves of three cultivars of foxtail millet included many WRKY, ethylene responsive factor (ERF), myeloblastosis (MYB), no apical meristem; Arabidopsis transcription activation factor; and cup-shaped cotyledon (NAC), basic helix-loop-helix (bHLH), basic leucine zipper (bZIP), and heat-shock transcription factor (HSF) TFs that have been confirmed to be involved in the transcriptional regulation of multiple stress responses in plants (Lata et al. 2013;Muthamilarasan et al. 2015;Todaka et al. 2015;Wang et al. 2018). The number of target genes regulated by the DOF family was the largest. In addition, in this study, the number of DEGs identified by DOF and C2H2 families in Jigu39 with the strongest drought resistance was higher than that of the other two cultivars, while the number of DEGs identified by the MIKC-type MADS family was the largest in Longgu16 with the weakest drought resistance (Supplement Table 5). We also found that some TFs had different expression patterns among the three foxtail millet cultivars in response to drought stress (Supplement Table 6). Notably, three TFs, including Seita.1G043800 (belonging to the bZIP family) and Seita.9G432200 (belonging to the MYB family), and Seita.2G206000 (belonging to the WRKY family), expressed opposite patterns in Jigu39X/Jigu39 and Jingu21X/Jingu21. The first two were down-regulated in Jigu39X/Jigu39 but up-regulated in Jingu21X/Jingu21, whereas Seita.2G206000 was the opposite. Of interest, none of the three TFs was expressed significantly in Longgu16X/Longgu16 (Supplement Table 6). In summary, the expression changes of the three cultivars under drought stress indicated that the TF encoding gene played an important role in the drought tolerance of foxtail millet.

Prediction of drought resistance genes in differentially expressed genes
According to the expression levels and differential expression multiples of the differential genes, we identified 46 genes, and the upward or downward trend of their expression levels under drought stress was consistent with the resistance trend of three cultivars (Table 2). These genes may play an In order to verify the results of DEGs in RNA-seq data, we randomly selected nine genes from the above 46 genes whose transcription levels are consistent with the drought resistance trend among three genotypes of foxtail millet and performed the qRT-PCR analysis (Supplement Table 1). The qRT-PCR results of the candidate genes are shown in Supplement Fig. 5. The expression patterns of these genes obtained by qRT-PCR were consistent with the results of RNA-seq (Fig. 9). Pearson correlation coefficient showed a significantly positive correlation between RNA seq and qRT-PCR data. In Jigu39X/Jigu39, r value is 0.814, in Jingu21X/ Jingu21, r value is 0.907, and in Longgu16X/Longgu16, r value is 0.912 (Fig. 9).

Discussion
Plants are challenged by various biotic and abiotic stresses in the whole life process (Boyer 1982). Crops respond to these stresses through biochemical and physiological processes (Shinozaki et al. 2003). Drought is one of the major abiotic stresses affecting agriculture. In recent years, in order to improve crop productivity in regions with limited water resources, the transcriptional regulation mechanism of drought response in many grain crops such as rice (Lenka et al. 2011) and Sorghum (Dugas et al. 2011), wheat (Ergen et al. 2009) and maize (Zheng et al. 2010) has been extensively studied. Foxtail millet, as a member of Gramineae plants, has outstanding tolerance to drought and has the potential to become a new C4 model plant (Peng and Zhang, 2021). Although the research on foxtail millet is increasing daily, we still know little about the complex network of drought response of foxtail millet.
Moreover, most transcriptional research is concentrated on one or two cultivars, and there is a lack of comparison data among multiple cultivars. Therefore, it is still quite challenging to deeply explore the drought-resistant mechanism of foxtail millet. This study performed transcriptome sequencing of three foxtail millet cultivars with different drought resistance to identify the transcriptome characteristics common to the three cultivars, thereby defining the core genes for drought resistance in foxtail millet, eliminating the specificity caused by differences in individual genetic backgrounds.

Three cultivars of foxtail millet had different drought tolerance mechanism
This study measured leaf physiological status changes and transcriptome changes of three foxtail millet cultivars after drought stress. From the physiological data, the drought resistance of Jigu39 was better than that of Jingu21 and Longgu16. Its water loss rate and chlorophyll content decreased mildly. Moreover, Its soluble protein content increased significantly (Fig. 1b, d and e). Under drought stress, the MDA content of three cultivars of seedling leaves increased to varying degrees (Fig. 1c). The above results supported previous reports that O 2 produced in chloroplasts could receive electrons from the photosynthetic electron transport chain and transform into O 2− under drought stress, causing oxidative damage to the photosynthetic pigment and plasma membrane (Gill and Tuteja 2010;Exposito-Rodriguez et al. 2017). In addition, to eliminate the difference in growth trend among millet cultivars, we studied the growth of three cultivars at different development stages during the Fig. 8 Statistical histogram of differentially expressed genes in amino acid metabolic pathways. The horizontal axis is the name of the three comparison groups. The longitudinal axis is the number of differential genes in the comparison group. Up is the number of up-regulated genes with significant differences, and down is the number of downregulated genes with significant differences    Fig. 9 Linear correlation analysis of fold change data between qRT-PCR and RNA-seq using Pearson correlation coefficient (r) experimental time. The results showed that the development states of the three cultivars were similar (Supplement Fig. 1). In summary, the study showed that among the three cultivars of foxtail millet, Jigu39 had the strongest drought resistance, and Longgu16 had the weakest drought resistance (Fig. 1a), which laid the foundation for our further research on drought tolerance mechanism by RNA-Seq analysis. RNA-seq analyzed the transcriptome of these three cultivars under drought stress, and many DEGs were detected in the three cultivars, which revealed the transcriptome reconstruction of foxtail millet under drought stress. In addition, according to GO and KEGG analysis of DEGs in Jigu39X/ Jigu39, Jingu21X/Jingu21, and Longgu16X/Longgu16, there were differences in molecular functions and pathways of DEGs. In order to further understand the differences in the regulatory mechanisms of three cultivars of foxtail millet in response to drought stress, we also conducted GO enrichment analysis on the overlapping and specific DEGs (Supplement Fig. 4). The study found that the DEGs with overlapping of the three genotypes were significantly enriched in the biological process, such as "recognition of pollen", "protein phosphorylation", "protein", and "response to stress" (Supplement Fig. 4a). Interestingly, we found that these DEGs significantly enriched GO terms related to photosynthesis and energy metabolism in cellular and molecular functions, respectively. The overlapping DEGs of the three cultivars might play a conservative role in the drought stress response of foxtail millet. In addition, it should be noted that the specific DEGs of the Jigu39X/Jigu39 were significantly enriched in the GO term of "translation", and the specific DEGs of the Jingu21X/Jingu21 were significantly enriched in the GO terms of "carbohydrate metabolic process" and "metabolic process". The specific DEGs of the Longgu16X/ Longgu16 were significantly enriched in the GO term of "embryo development" (Supplement Fig. 4). These specific genes were enriched into different GO terms, indicating that different cultivars of foxtail millet had different regulatory mechanisms in response to drought stress.
Based on KEGG pathway analysis, we found that the up-regulated DEGs in Jigu39X/Jigu39 were significantly enriched in "ribosome" (ko03010), "ribosome biogenesis in eukaryotes" (ko03008). The ribosome is the molecular factory of protein synthesis, and ribosome biosynthesis is the basis of cell growth ability (Lempiainen and Shore 2009). This result indicated that Jigu39 could ensure the correct assembly of ribosomal proteins under drought stress. Among three cultivars, fewer DEGs were identified in Jin-gu21X/Jingu21. The DEGs of Jingu21X/Jingu21 involved in chlorophyll metabolism, amino acid metabolism, and carbohydrate metabolism pathways were fewer than those in Jingu21X/Jingu21. There were a large number of DEGs enriched in some amino acid metabolic pathways in Long-gu16X/Longgu16, such as "glycine, serine and threonine metabolism" (ko00260), "cysteine and methionine metabolism" (ko00270), "valine, leucine and isoleucine degradation" (ko00280), "arginine and proline metabolism" (ko00330), "alanine, aspartate and glutamate metabolism" (ko00250).
In conclusion, different physiological and transcriptomic changes of the three cultivars suggested that they have different mechanisms of the drought stress response.

Drought response pathways and related genes in foxtail millet
Photosynthesis is important for plants accumulating organic matter (Xiong et al. 2019). Based on the GO enrichment results, it was found that many down-regulated DEGs were enriched in the GO term of "photosynthesis" (GO:0015979) in the three cultivars of foxtail millet, indicating that drought stress would affect the photosynthesis (Supplement Table 3), consistent with previous reports (Galmes et al. 2007;Chaves et al. 2009;Zhang et al. 2016;Shi et al. 2018). We have sorted out the DEGs that are significantly enriched in the pathway of photosynthesis (Fig. 5). Cluster analysis showed that the transcriptional regulation of the photosynthesis pathway in Jigu39 and Longgu16 was similar (Fig. 5). The studies found that the transcription levels of most genes related to photosystem II (PSII), cytochrome b6/f (Cytb6f), photosystem I (PSI) and carbon fixation were down-regulated in Jigu39 and Longgu16 (Fig. 6). However, compared with Longgu16, the regulation mechanism of photosynthesis is more complex in Jigu39. The transcription levels of psbA (Seita.J017900), psbB (Seita.J014000), psbK (Seita.1G126100), psbT (Seita.J014100), psbZ (Seita. J018400), psaA (Seita.J011100), petG (Seita.J013400), and atpB (Seita.J020700) were only down-regulated in Jigu (Fig. 6). Further analysis showed that under drought stress, the transcription level of the PsaO subunit responsible for binding to light-harvesting complex II in the PSI reaction centre was significantly down-regulated in three cultivars of foxtail millet (Gorski et al. 2022). The transcription levels of PNSL3, a subunit of the chloroplast NADH dehydrogenase-like (NDH) complex involved in photosystem I (PSI) cycle electron flow, as well as FdC1 and FdC2, two FD-like protein members encoded by the plant genome with an extension at the C-terminal of the 2Fe-2S cluster, were significantly down-regulated in the drought-sensitive Longgu16. Therefore, we speculated that the down-regulation of transcription levels of PsaO, PNSL3, FdC1 and FdC2 in drought-sensitive Longgu16 was higher than that in Jigu39, which might be one of the reasons why the drought resistance of Longgu16 was weaker than Jigu39 (Fig. 5). It indicated that chloroplast photosynthesis was involved in plant response to drought stress.
Photosynthesis is related to photosynthetic pigments, photosynthetic gene expression, and redox regulation (Li and Yi 2020). Chlorophyll is the main photosynthetic pigment in plants and plays a key role in the light absorption of photosynthesis. A stable supply of chlorophyll is the key factor for the plant's development of photosynthesis. (Xiong et al. 2019). In our study, the down-regulated DEGs in Jigu39X/Jigu39 and Longgu16X/Longgu16 were significantly enriched in the 'Porphyrin and chlorophyll metabolism' (ko00860) pathway, and the number of DEGs involved in this pathway in Jingu21 was the least (Fig. 8, Supplement Table 4). Under drought stress, the expression levels of protein-coding genes involved in chlorophyll metabolism in Jigu39X/Jigu39 and Longgu16X/Longgu16 showed significant synchronous changes (Fig. 8). Except for the up-regulation of some enzymes, such as heme a synthase (COX15) and chlorophyll b reductase (NOL, NYC1), the expression levels of most protein-coding genes were down-regulated in all three cultivars. The pheophorbide an oxygenase (PAO) gene was significantly up-regulated in Longgu16X/Longgu16, which played an important role in chlorophyll degradation. The genes of protochlorophyllide reductase (PORA) and magnesium protoporphyrin IX methyltransferase, chloroplastic (CHLM) were significantly down-regulated, which were related to chlorophyll biosynthesis. This result may be one of the reasons for the relatively weak drought resistance of Longgu16 (Fig. 8).
Carbohydrates are substrates for energy production, osmoregulation, and osmoprotection in response to drought and are associated with growth, development, and carbon status. Among the main enriched carbohydrate metabolic pathways, the numbers of DEGs in Jigu39 and Longgu16 were similar in each pathway, while Jingu21 was relatively small (Table 1).
Amino acids are the synthetic units of proteins and play a role in compound biosynthesis, signalling, and stress responses Ming et al. 2021). In plants, the levels of free amino acids are tightly controlled by metabolic and developmental factors. Their degradation not only contributes to osmoregulation to maintain cellular homeostasis but also to plants' energy state, providing a link between carbon and nitrogen metabolism. Under drought stress, amino acid and carbohydrate metabolism are important for stress adaptation (Diniz et al. 2020). This study also sorted out the DEGs involved in amino acid metabolism under drought stress (Supplement Table 7). The number of DEGs in Longgu16 was the largest, while that in Jingu21 was the smallest. After comparison, only 21 overlapping genes (5.98%) were identified in three foxtail millet cultivars, and all DEGs were down-regulated. This situation shows that three cultivars have different regulatory mechanisms in amino acid metabolism under drought stress, which may be one of the reasons they show different tolerance to drought stress.

Transcription factors play an important role in drought tolerance
Transcription factors control the expression of many target genes in plants under drought stress (Qin et al. 2020). In this study, we found that the number of target genes regulated by the DOF family was the most among three cultivars of foxtail millet, followed by the ERF family (Supplement Table 5), indicating that the transcription factors of these two families played an important role in the drought response of foxtail millet. The DOF family involves many basic processes in higher plants, from response to light and plant hormones to flowering time and seed maturation . Studies have shown that the DOF and MYB families play important regulatory roles in barley seed germination (Seven and Akdemir 2020). However, the relationship between DOF-type TFs and abiotic stress tolerance remains largely unknown (Corrales et al. 2017). The ERF family is a large gene family of transcription factors (Nakano et al. 2006). In contrast to the DOF family, many proteins of the ERF family were identified in our study. They are involved in many different functions in plant growth and development, such as hormone signal transduction and disease resistance pathways (Gutterson and Reuber 2004), responses to biotic and abiotic stresses (Gu et al. 2000;Gutterson and Reuber 2004;Mizoi et al. 2012), and regulation of metabolism (van der Fits and Memelink 2000; Aharoni et al. 2004;Broun et al. 2004;Zhang et al. 2005). We identified 5 DOF and 12 ERF differentially expressed in the leaves of three cultivars. The number of differentially expressed transcription factors (DETFs) encoded by DOF and ERF in the tolerant cultivar (Jigu39) was higher than that in the sensitive cultivars (Jingu21 and Longgu16) (Supplement Table 5). In Jigu39, DETFs of four DOF families were all downregulated. The only existed DOF (Seita.1G290700) in the other two cultivars was up-regulated. Its homology was presumed to be a key regulator involved in carbohydrate metabolism and drought tolerance in Betula platyphylla (Tarelkina et al. 2020). In addition, ERF family proteins are considered to be key regulators in response to environmental stimuli in Arabidopsis and rice (Nakano et al. 2006). It is worth noting that two up-regulated DETFs (Seita.1G342700 and Seita.2G230500) and three down-regulated DETFs (Seita.1G055900, Seita.8G081600 and Seita.9G320700) of this family only exist in drought tolerance cultivar (Jigu39). They were likely to play an important role in the drought tolerance of Jigu39. Previous research has shown that the expression level of EREBP1 (homology of Seita.1G342700) 1 3 in drought-resistant cultivars of Lemon Balm was significantly higher than that in drought-sensitive cultivars (Rahimi et al. 2016).
In addition to the above genes, in previous studies, some genes were identified as important regulatory factors involved in plant abiotic stress responses in some plant species, but they have not been reported in foxtail millet. For Example, ACO1 (homology of Seita.1G335200) is related to drought resistance in Arabidopsis and tomato. Plants can regulate the ethylene synthesis in plants by increasing their expression levels, thereby improving drought resistance of plants (Eglous et al. 2013). FTSH6 (homology of Seita.4G099200), identified as a potential drought tolerance marker in Quercus ilex (Guerrero-Sanchez et al. 2021), was highly induced when the tomato was subjected to water stress (Tapia et al. 2021). The expression of R40C1 (homology of Seita.5G069300) and rpl20 (homology of Seita.J013800) were related to the drought tolerance potential of rice and Caragana cultivars, respectively (Jiang et al. 2018;Sahid et al. 2020). In cotton, GT-2 (homology of Seita.7G184300) was highly up-regulated under drought and salt stress. Overexpression of GT-2 enhanced drought tolerance, decreased MDA and H 2 O 2 contents, and increased reactive oxygen species scavenging enzyme activity in cotton (Magwanga et al. 2019). In Arabidopsis, plants expressing CBP (homology of Seita.8G025700) have a better survival rate and root growth under intermittent drought or high salt stress (Tsou et al. 2012).
Lipoxygenase catalyzes the production of oxytocin and plays a key role in the whole growth period of plants (such as seed germination, plant growth, and stress resistance). A previous study by Zhang et al. (2021) showed that the transcription level of Seita.7G113700 was decreased in foxtail millet under drought stress, which was different from our results. In addition, in our study, HSP18.6 (Seita.1G342400) was involved in the response process of foxtail millet to drought stress, which was also different from the result reported by Coca et al. (1996) in sunflowers. These results indicated that the transcriptional regulation of the same gene was different in different species or even among different cultivars under drought stress. As the drought resistance mechanism of foxtail millet is relatively complex and involves various factors, these problems need to be further studied.
We also found some drought-related genes in foxtail millet ( Table 2) that previous studies have not reported as stressrelated. Among these identified new genes, the functions of some genes' homologous in other species have been characterized, such as DMAS1-D (homology of Seita.7G115200) and ZIFL1 (homology of Seita.7G312500), they maintain iron homeostasis in gramineous plants such as rice and barley, and CYP701A6 (homology of Seita.4G235500) necessary for GA biosynthesis in rice Beasley et al. 2017;Lee et al. 2021). However, there are still some genes, such as Seita.7G230500, Seita.4G155500 and Seita.6G096800, whose functions in plant growth and development need further study. In conclusion, these genes were probably related to the drought resistance mechanism of foxtail millet.

Conclusion
The plant drought stress response is complex, involving survival optimization at multiple levels (Diniz et al. 2020). The physiological and transcriptional changes of three foxtail millet cultivars were studied to understand the drought stress response mechanism. We found that these three cultivars had different drought tolerance. Jigu39 has the strongest drought tolerance, followed by Jingu21, and Longgu16 was the weakest. This study performed transcriptome sequencing on these three cultivars to identify genes whose transcriptional level change trends were consistent with drought resistance trends. From the deep transcriptome investigation of these three cultivars of foxtail millet, we found that foxtail millet initiated a system composed of multiple genes and multiple steps, which regulated various pathways to resist drought stress. GO enrichment showed that differentially expressed genes were mainly enriched in the GO terms related to photosynthesis, amino acid metabolism, and carbohydrate metabolism, indicating that some key DEGs involved in these pathways played an important role in drought resistance foxtail millet. Based on the analysis of the KEGG pathway, we found differences in the pathways of apparent enrichment in the three cultivars of foxtail millet. In addition, we found that two transcription factor families, DOF and ERF, played important roles in the drought response of foxtail millet. We identified 46 genes whose expression levels shared the same trends with the drought resistance of the three foxtail millet cultivars. Among them, 32 genes were firstly reported as drought stress response-related genes in foxtail millet. Our results provide a theoretical basis for further research on drought-resistant mechanisms.