We selected 10 S. italica accessions from a previous study (Jia et al., 2013) for field assessment under normal and water-limited conditions during the jointing stage. The 10 accessions and the decrease of phenotypes are introduced in Table 1, Figure 1A, respectively. We selected Ci134 and Ci603 as drought-sensitive genotypes and Ci328 and Ci409 as drought-tolerant genotypes for this trial.

Drought-sensitive accessions (Ci134 and Ci603) and drought-tolerant accessions (Ci328 and Ci409) were planted in pots (diameter, 30 cm; depth, 40 cm) with five uniform individuals selected from 10 seedlings in each pot after 15 days of seed sowing. All of the pots were placed in the greenhouse at the Institute of Crop Sciences, Chinese Academy of Agricultural Sciences (Beijing, China). After 45 days of growth under normal conditions with weekly watering, we selected uniformly developed and healthy individuals for further analysis. For the drought treatment, we randomly selected five pots from each accession to stop watering for 2 weeks, while we continued to water the control individuals every week until 3 days before sampling (Supplementary Figure S2A).

For the RNA-seq, we collected two upper, fully expanded leaves. We also harvested young, elongating stems from the same plant after removing the sheaths and leaves. We collected all of the roots from the same individual (Supplementary Figure S2B). Finally, three individuals were harvested for each accession as biological replicates. All samples were then immersed in liquid nitrogen immediately after collection and stored at −80°C until total RNA extraction.

We measured the photosynthesis rate and the stomatal conductance using the Li-6400 platform (LI-COR Inc., USA). We used the WP4-T Dewpoint Potential Meter (ICT International Inc., USA) to measure leaf water potential. The leaf soluble sugar content was measured using the anthrone colorimetry method described in a previous study (Wen et al., 2005).

We sampled fully expanded leaves and weighed (Wf) them to measure relative leaf water content (RLWC). After immersing the leaves in water for 4 h, we removed the surface water (Wt) and weighed the samples. The samples dried in an oven until a constant weight (Wd) was reached. The RLWC was calculated using this formula:

We extracted the total RNA using a TRIZOL kit (Ambion, Thermo Fisher Scientific Inc., USA) according to the protocol described by instructions. The cDNA library was constructed and then sequenced using the Illumina HiSeq 2100 platform with 150-bp paired-end mode (Illumina Inc., USA), with three biological replications for each sample. Raw data obtained by this study are available on the EMBL-EBI database (http://www.ebi.ac.uk/) with accession No. PRJEB43702.

After removing the adapter and filtering-quality reads, we mapped clean reads to a reference genome for S. italica V2.2 (http://phytozome.jgi.doe.gov) using Hisat2 tools (Kim et al., 2015). All transcripts were assembled using StingTie (Pertea et al., 2015). The gene expression level was quantified by FPKM (Fragments Per Kilobase of exon model per million mapped fragments) and the differentially expressed genes (DEGs) were identified using DEseq2 (Love et al., 2014). The screening threshold for DEGs was fold-changes bigger than two, and the false discovery rate was < 0.05.

Identification of DEGs was conducted through two different approaches. First, the drought condition was compared with normal condition to identify drought induced/repressed DEGs, which were then used for detection of genotype/tissue-specific DEGs, differentially expressed transcription factors (DETFs), and k-means cluster analysis. The transcription factors (TF) list was obtained from PlantTFDB (http://planttfdb.gao-lab.org/). Second, drought-tolerant genotypes were compared with sensitive genotypes to find out DEGs between genotypes under the same condition. After elimination of the genotype differences under normal condition, DEGs identified under drought condition were defined as drought-induced DEGs between genotypes (di-DBG).

We randomly selected the DEGs for validation through qRT-PCR, using primers listed in Supplementary Table 15. We selected Cullin and EF1a as reference genes in stems and leaves and SDH and EF1a as reference genes in roots. The first strand of cDNA was synthesized using 1 μg of total RNA with a 20-μl reaction mixture and then was diluted to a 200-μl volume. We used 2 μl of cDNA as a template for the qRT-PCR reaction (95°C, 10 min; then 40 cycles of 95°C, 10 s; 60°C, 60 s). Four replications were conducted for each gene and calculated the relative expression levels using the 2^−ΔΔCT method.

We conducted functional annotation of transcripts using Gene Ontology (GO, https://geneontology.org) and the Kyoto Encyclopedia of Genes and Genomes (KEGG, https://www.genome.jp/kegg) database. Gene sets functional enrichment analysis was performed with clusterProfiler (Yu et al., 2012), and enrichment results were visualized using ggplot2 (Wickham, 2016).

A total of 9,237 DEGs detected in roots, 4,276 DEGs detected in stems, and 6,923 DEGs detected in leaves were used as inputs for K-means cluster analysis, respectively, after being conversed by log₁₀ (FPKM + 1), and 6 clusters were finally generated from each DEG set.

We filtered all gene expression data by FPKM > 0.1 in at least eight samples. After normalizing the data by log₂ (FPKM + 1), the weighted gene co-expression network was constructed using the WGCNA package (Langfelder and Horvath, 2012) with the default parameters except that power = 5, maxBlockSize = 5,000. We exported and visualized the network using Cytoscape software (version 3.7.2, https://cytoscape.org/).

We used the foxtail millet core collections reported in our previous study (Jia et al., 2013) to select drought-sensitive and drought-tolerant accessions. Many elite cultivars of S. italica exhibit exceptional drought tolerance abilities under water-limited conditions and rarely decrease grain yield under osmotic conditions, while drought-sensitive accessions incur a heavy grain yield loss for adjusting similar water-deficient conditions. In this trial, we screened out 10 accessions according to the grain-weight decreased ratio (1 – drought/control ×100%), and then we evaluated them in field experiments under normal and drought conditions (Figure 1A). To our surprise, many drought-tolerant genotypes of S. italica that we identified at the jointing stage were not consistent with those identified at the seedling stage. For instance, Ci488 is most vulnerable, and Ci309 is generally tolerant to drought conditions at the seedling stage, according to a previous study (Tang et al., 2017), but the grain weight and growth of Ci488 (Figure 1A) were only slightly affected by drought at the jointing stage in the field. This observation indicates that the drought-tolerant abilities of many foxtail millet varieties might not remain stable during the whole development period. Similar results were also characterized in other grass species like sorghum (Varoquaux et al., 2019). Compared with normal conditions, grain weight and plant height of Ci309 and Ci328 were decreased by about 6% under drought environments. Grain weight, and panicle weight of Ci409 was stable, and plant height of Ci409 was decreased by 22% under osmotic conditions compared with the control. Grain weight, panicle weight, and plant height also decreased in Ci134 (≈50%), Ci603 (≈60%), and Ci351 (≈70%) under drought conditions, respectively (Figure 1A). We also evaluated heading dates of relevant accessions (Ci235 and Ci426 were much longer; Ci035 was shorter than others. Data not shown) for accurate determination of developmental stages for drought treatment and sample harvesting. Finally, we selected Ci134 and Ci603 as drought-sensitive genotypes. We selected Ci328 and Ci409 as drought-tolerant genotypes for further RNA-seq analysis.

After 2 weeks of drought treatment, when the soil-water content fell below 15% (Table 2), we sampled from the drought-tolerant genotypes (Ci328 and Ci409) and the drought-sensitive genotypes (Ci134 and Ci603) for RNA-seq. Relative leaf water content (RLWC) of drought-sensitive and drought-tolerant genotypes was comparable (≈95%) under normal conditions, while we observed a difference in the RLWCs between drought-sensitive and drought-tolerant plants after 2 weeks of drought treatment. The RLWCs of two drought-sensitive genotypes decreased by ≈35% (35.8% for Ci134 and 34.1% for Ci603), but drought-tolerant genotypes only decreased by ≈7.5% (7.3% for Ci328 and 7.6% for Ci409) under osmotic-stress conditions (Figure 1B). Leaf water potential (LWP) of drought-tolerant genotypes was higher than drought-sensitive genotypes under drought conditions, and the change of LWPs in tolerant genotypes toward drought treatment was smaller than that of drought-sensitive genotypes (Figure 1C). Stabilities of soluble sugar content in leaves of drought-tolerant genotypes were also higher than drought-sensitive genotypes under drought treatments, implying osmatic adjustments are an essential propagation strategy for optimizing drought-tolerance of foxtail millet under water-limited conditions (Figure 1D). In addition, the photosynthesis rate (Figure 1E) and the stomatal conductance (Figure 1F) of all four accessions under drought treatments were significantly lower than the controls, with tolerant genotypes more stable than sensitive genotypes.

Figures 1G–J illustrate the morphological changes of leaves in sampled accessions after drought treatments. Fully expanded leaves of drought-sensitive genotypes (Ci134 and Ci603) bleached earlier and more severely than drought-tolerant genotypes under similar drought conditions. Notably, fully expanded leaves of the Ci134 genotype died due to dehydration, indicating Ci134 is sensitive to osmotic stress conditions. Based on these findings, we confirmed that Ci134 and Ci603 are drought-sensitive accessions, while Ci328 and Ci409 are drought-tolerant accessions at jointing developmental stages of foxtail millet.

To achieve a better understanding of drought-tolerant mechanisms of S. italica during the jointing stage, we conducted RNA-seq for roots, stems, and leaves of two drought-sensitive (Ci134 and Ci603) and two drought-tolerant genotypes (Ci328 and Ci409) under normal and drought conditions (Supplementary Figure S2). We constructed 72 pair-end libraries for transcriptomic sequencing and generated over 222.24 GB of raw reads. We obtained 1,682,873,236 clean reads and successfully mapped 89.95% of the reads to a reference genome of foxtail millet (Sitalica_312_v2.2, https://phytozome.jgi.doe.gov/) (Supplementary Table S1). The number of reads mapped to annotated genes ranged from 1 to ≈1.5 million, with a median of 352 (Figure 2A). We detected a total of 32,170 expressed genes in all samples (Supplementary Table S2), and the log₁₀(FPKM) value for all samples was a spread tally with a normal distribution (Figure 2B). We performed principal component analysis using expressed genes shared by all 72 samples (Figure 2C). All samples were divided into three main groups: roots, stems, and leaves. Additionally, each group was then separated into normal and drought treatment subgroups. All three biological replicates were mostly closed based on principal component analysis.

In this trial, we screened out expressed genes with more than 2-fold change (FC) and with a false discovery rate (FDR) of less than 0.05 as differentially expressed genes (DEGs). We separately analyzed three tissues of two drought-tolerant (Ci328 and Ci409) and two drought-sensitive genotypes (Ci134 and Ci603) under both drought and normal conditions at the jointing stage. A total of 13,552 DEGs correlated with drought tolerance among twelve comparison groups were identified (Table 3, Supplementary Table S3). We found the most DEGs in roots, followed by leaves, and then stems of all four genotypes. Overall, the DEGs detected in roots (9,237 DEGs, including 3,497 in Ci134, 6,007 in Ci603, 5,656 in Ci328, and 4,583 in Ci409) were 33% higher than in leaves (6,923 DEGs, including 2,557 in Ci134, 3,549 in Ci603, 3,776 in Ci328, and 3,117 in Ci409) and more than 2-fold higher than in stems (4,276 DEGs, including 2,013 in Ci134, 2,229 in Ci603, 1,454 in Ci328, and 587 in Ci409) (Table 3, Figure 2D). These observations suggest that the roots are more sensitive to drought conditions, consistent with observations reported in rice (Wang D. et al., 2011). Moreover, the number of DEGs in the stems of two drought-sensitive genotypes (2,013 DEGs in Ci134 and 2,229 DEGs in Ci603, respectively) was higher than that of two drought-tolerant genotypes (1,454 DEGs in Ci328 and 587 DEGs in Ci409), suggesting stems might play an important role in determining drought tolerance in foxtail millet (Table 3).

To clarify the biological functions of the DEGs responses to drought in foxtail millet, we performed KEGG and GO enrichment analysis for each DEGs set. The functional enrichment analysis of DEGs detected in roots revealed osmotic adjustment and phytohormone-related pathways, playing important roles in drought-tolerant mechanisms of S. italica. KEGG enrichment analysis revealed “starch and sucrose metabolism,” and “galactose metabolism” pathways were enriched in upregulated DEGs, and “phenylpropanoid and diterpenoid biosynthesis” pathways were enriched in downregulated DEGs identified in roots of foxtail millet. As a whole, “glycolysis and glutathione metabolism” pathways were enriched in all DEGs detected in roots of all four accessions (Supplementary Figure S3A, Supplementary Table S4). Moreover, GO enrichment analysis also revealed that regulatory genes participated in responding pathways of heat, water, and growth hormones were obviously affected by drought stress (Supplementary Figure S3B, Supplementary Table S5). Interestingly, two ABA pathway-related genes of ABF4 (Seita.1G331600) and ABA2 (Seita.9G041900) were identified as downregulated DEGs in drought-sensitive genotypes but upregulated DEGs in drought-tolerant genotypes.

Functional categories of the DEGs identified in the stems of drought-sensitive (Ci134 and Ci603), and drought-tolerant (Ci328 and Ci409) genotypes revealed that young stem elongation was affected in drought-sensitive genotypes, and heat shock proteins (HSPs) might contribute to the growth of young stems of foxtail millet under drought conditions. Regulatory genes involved in plant cell organizations were primarily enriched in the DEGs set of drought-sensitive genotypes, and the number of downregulated DEGs of this pathway was higher in drought-sensitive genotypes than those that are drought tolerant. This observation implies that stem elongation of foxtail millet under drought conditions was seriously affected, which coincides with the plant height of both drought-sensitive and drought-tolerant accessions detected in this trial under normal and drought conditions (Figure 1A). KEGG enrichment analysis revealed “plant hormone signal transduction” and “amino sugar and nucleotide sugar metabolism” pathways were enriched in upregulated DEGs detected in the stems of drought-tolerant genotypes (Supplementary Figure S4A, Supplementary Table S4). Interestingly, heat and reactive oxygen species response GO terms were only enriched in drought-tolerant accessions (Supplementary Figure S4B, Supplementary Table S5), and five HSPs (Seita.3G175700, Seita.5G131400, Seita.9G470600, Seita.9G488400, and Seita.2G241600) were identified as upregulated DEGs in drought-tolerant genotypes under stressed conditions, while no significant change was detected in drought-sensitive genotypes. Therefore, we speculate that drought-sensitive genotypes might reduce chaperone synthesis, which could support elongation of stems under drought conditions through a transcription regulatory mode.

Functional enrichment analysis of drought-response DEGs in leaves of all genotypes revealed a clear signature that photosynthesis-related genes were downregulated, but plant hormone transduction-related genes were upregulated under drought conditions (Supplementary Tables S4, S5). Moreover, the “starch and sucrose metabolism” pathway was enriched for downregulated DEGs in sensitive genotypes. In contrast, the “glycerol phospholipid metabolism” pathway was enriched for upregulated DEGs in drought-tolerant genotypes (Supplementary Figure S5), suggesting photosynthetic characters of functional leaves were seriously affected under water-limited conditions.

In this investigation, we detected most DEGs in a tissue-specific manner. In total, only 243, 415, 296, and 134 DEGs out of 6,481, 9,184, 8,497, and 6,953 DEGs were shared in three tissues of Ci134, Ci603, Ci328, and Ci409, respectively (Table 3, Figures 3A–D). There were 551, 1,520, 1,256, and 1,115 upregulated DEGs, 2,011, 3,224, 3,026, and 2,574 downregulated DEGs from roots of Ci134, Ci603, Ci328, and Ci409. Moreover, there were 851, 1,099, 1,010, and 887 induced DEGs, and 893, 1,171, 1,348, and 1,364 repressed DEGs in the leaves of Ci134, Ci603, Ci328, and Ci409 under drought conditions. In stems, there were 386, 350, 223, and 45 upregulated DEGs, and 799, 736, 365, and 188 DEGs were downregulated in Ci134, Ci603, Ci328, and Ci409 under osmotic stress conditions. In total, few of the DEGs were up or downregulated under drought conditions in all three tissues, including 144, 197, 174, and 94 induced DEGs and 88, 130, 78, and 32 repressed DEGs in all four accessions (Ci134, Ci603, Ci328, and Ci409) (Supplementary Figures S6A–D, S7A–D).

We also identified genotype-specific DEGs in this trial. In roots, there were 163 upregulated DEGs and 198 downregulated DEGs only detected in drought-tolerant accessions, while we identified 55 upregulated DEGs and 154 downregulated DEGs in drought-sensitive accessions (Supplementary Figures S6E, S7E). In stems, there were 26 upregulated DEGs and 10 downregulated DEGs only identified in drought-tolerant genotypes, while we only detected 73 upregulated DEGs and 217 downregulated DEGs in drought-sensitive accessions (Supplementary Figures S6F, S7F). In leaves, there were 237 upregulated DEGs and 381 downregulated DEGs only detected in drought-tolerant accessions, while we only identified 194 upregulated DEGs and 253 downregulated DEGs in drought-sensitive genotypes (Supplementary Figures S6G, S7G). The roots of all four genotypes shared a total of 1,727 DEGs (out of 9,237). The leaves of all four accessions shared a total of 786 DEGs (out of 6,923). Furthermore, the stems of all four genotypes shared a total of 123 DEGs (out of 4,276) (Table 3, Figures 3E–G).

To dissect the function of genotype-shared/specified DEGs, KEGG and GO enrichment analysis was then performed (Supplementary Figure S8, Supplementary Table S6), which revealed that water deprivation-related genes were differentially expressed in all four genotypes under drought conditions, and only few of water-deprivation DEGs were shared by different tissues. Moreover, photosynthesis-related genes were affected by drought in leaves of all four accessions, and plant hormone transduction and MAPK-related DEGs were also enriched in leaves and stems.

We detected 35 upregulated and 3 downregulated genes shared by all three tissues of the four accessions under drought treatments (Table 4). Previous studies have confirmed that many functional genes of this DEGs set are correlated with the response of plants to multiple abiotic stress conditions, such as salt, drought, and extreme temperature stresses (Wu et al., 2005; Ali et al., 2013; Nawkar et al., 2017). For instance, 5 LEA proteins, 5 HSPs (HSP70, HSP70T-2, HSA32, and HSP17.6C), 2 PP2Cs (HAI3 and AHG1), and 1 dehydrin family protein (RAB18) were upregulated under drought conditions. Previous study has confirmed that these proteins act as abiotic stress-responding factors (Table 4).

Notably, we detected 21 DEGs (including 4, 1, and 16 DEGs identified in roots, stems, and leaves, respectively), responding to drought with the opposite pattern between sensitive and tolerant genotypes in this study (Figure 4, Supplementary Table S7). Two heat shock protein-encoding genes [Seita.3G216900 and Seita.1G326600 (HSP23.5)] were upregulated in drought-tolerant genotypes but were downregulated in drought-sensitive genotypes under drought conditions. Heat shock protein HSP23.5 has also been confirmed to be induced by salinity and heat conditions in Arabidopsis (Sewelam et al., 2019). Moreover, one heat shock transcription factor (HSFA6B, Seita.9G095900)-has been characterized, and knock-out of HSFA6B in Arabidopsis has elevated ABA sensitivity and stress tolerance of transgenic plants (Wang et al., 2020). Seita.5G031200 (CIPK1) encodes a CBL-interacting protein kinase involved in an ABA-dependent regulatory network in the plant (D'Angelo et al., 2006). In addition, BAG5 (Seita.8G128400) is a known leaf-senescence-related regulator in Arabidopsis (Fu et al., 2017). We also found BAG5 was downregulated in the roots of sensitive genotypes and upregulated in tolerant genotypes in this trial, implying that BAG5 plays a vital role in drought adjustment of the roots in S. italica. Moreover, one chloroplast thylakoid-localized RbcX protein (RBCX1, Seita.6G151500), a chaperone in the folding of Rubisco, and two pentatricopeptide repeats (PPR like), protein-encoding genes (LPE1, Seita.3G308800; Seita.7G212200) were also identified as DEGs in the leaves of foxtail millet under drought conditions. These oppositely regulated DEGs between sensitive and tolerant accessions may play key roles in modulations of drought-tolerant abilities of S. italica during the jointing stage.

Transcription factors play central roles in dealing with abiotic stress conditions through regulating functional genes in plants (Baillo et al., 2019). We identified a total of 915 TFs from the DEGs detected in this trial, including many WRKY, ERF, MYB, NAC, bHLH, bZIP, and HSF transcription factors that have been confirmed to be important for the tolerance of biotic and abiotic stresses in rice (Todaka et al., 2015). These differentially expressed transcription factors (DETFs) identified in foxtail millet during the jointing stage (Figure 5) were similar to observations reported in previous studies of drought tolerance in foxtail millet at the seedling stage (Tang et al., 2017). We detected 200 DETFs shared by all four accessions (Supplementary Figure S9), including 131 downregulated DETFs and 69 upregulated DETFs in this trial. We identified a total of 35 ERF, 35 WRKY, 32 bHLH, 31 MYB, 30 NAC, 23 C2H2, and 18 bZIP DETFs differentially expressed in roots, in at least two genotypes, and 10 MYBs, 7 ERFs, 6 HD-ZIPs, 5 bHLHs, and 5 bZIPs, and we detected 4 NAC family DETFs differentially expressed in stems in more than two accessions, and we identified 18 ERF, 15 MYB, 13 bZIP, 12 NAC, and 12 HD-ZIP family DETFs differentially expressed in leaves in this trial.

Additionally, we identified 127 DETFs (including 51 in roots, 26 in stems, and 50 in leaves) with different expression patterns in response to drought conditions between drought-tolerant and drought-sensitive accessions (Supplementary Table S8). For instance, Seita.6G055700 (LATE ELONGATED HYPOCOTYL) encodes an MYB-like transcription factor, which represses the biosynthesis of ABA in plants (Adams et al., 2018); another MYB transcription factor-encoding gene (Seita.2G199900) identified as a drought-responding regulator in this study has also been confirmed as a negative regulator of ABA signal transduction; overexpressing of which could enhance salt tolerance of Arabidopsis individuals (Cui et al., 2013). Furthermore, in this study, the number of ERF-encoded DETFs identified in tolerant genotypes was higher than that of sensitive ones. For instance, an ERF-encoding gene (Seita.9G500100, ERF-IXc5) has been inferred as a crucial regulator of plant growth and drought tolerance of Hevea brasiliensis (Lestari et al., 2018). Moreover, ERF family proteins have been considered as crucial regulators responding to environmental stimuli in Arabidopsis and rice (Nakano et al., 2006). Notably, two DETFs, including Seita.1G228800 (CO-like family) and Seita.9G095900 (HSF family), were downregulated in the roots of sensitive genotypes but upregulated in tolerant genotypes. Overall, DETFs identified in this trial also suggested that TFs have played essential roles in the drought-tolerant process in foxtail millet during the jointing stage.

To further dissect biological processes and pathways correlated with drought tolerance in each tissue, a k-means clustering analysis of DEGs detected in leaves, stems, and roots were carried out, respectively, and 6 clusters were generated for each DEG set (Figure 6, Supplementary Table S9). Functional enrichment analysis of each cluster (Supplementary Table S10) revealed that different biological pathways were responsible for drought tolerance in the three different tissues.

In this trial, the biological processes of “water deprivation response” (drought) and “response to heat” (heat) were enriched in many clusters identified in three tissues, which implied that drought and heat are both severe environmental stresses often occurring in the field together. Many clusters, including Cluster B in roots (46 DEGs involved in drought and 33 DEGs involved in heat), Cluster A in stems (30 DEGs involved in drought and 22 DEGs involved in heat), and Cluster B in leaves (37 DEGs involved in drought and 27 DEGs involved in heat), were involved in these two pathways (Figure 6, Supplementary Table S10), which were consistent with upregulated DEGs identified under drought treatments. However, only 11 shared DEGs involved in the “response to water deprivation” pathway and 7 shared DEGs involved in the “response to heat” pathway were identified in all three tissues, implying that different tissues might take different actions to achieve the similar goal of drought tolerance in foxtail millet at the jointing stage.

In leaves, “MAPK signal transduction” and “plant hormone signal transduction” were mainly enriched in Clusters A and B, with DEGs induced or repressed, respectively, by drought treatments. Moreover, “nitrogen metabolism” and “photosynthesis” related pathways were enriched in Cluster D with downregulated DEGs influenced by drought conditions, implying that enhanced signal transduction and decreased efficiency of photosynthesis and energy metabolism in leaves are essential for drought tolerance of foxtail millet under water-limited conditions.

In stems, stress-related pathways, such as “response to water deprivation” and “heat, salt, and hydrogen peroxide,” were enriched in Clusters A and F. Pathways of “secondary cell wall biogenesis” and “xylan biosynthetic process” were also enriched in Cluster F with downregulated DEGs identified under drought conditions, indicating suppression of stem cell growth in foxtail millet under drought conditions during the jointing stage. Furthermore, the “Auxin-activated process” and “secondary cell wall biogenesis” pathways were enriched in Cluster D with downregulated DEGs identified under drought conditions, which was consistent with the fact that plant height decreased under drought conditions through suppressing stem elongation at the jointing stage in foxtail millet. Photosynthesis-related biological pathways, such as “light-harvesting process” and “antenna proteins,” were enriched in Cluster C with downregulated DEGs, and “light reaction” and “carbon fixation” categories were enriched in Cluster E with upregulated DEGs identified under drought conditions, implying that stems were also a photosynthesis-related organ during the jointing stage in foxtail millet.

In roots, “plant-pathogen” categories were enriched in Cluster A with downregulated DEGs under drought conditions, which was consistent with conclusions that the mass of arbuscular mycorrhizal symbiosis with roots decreased under drought conditions in sorghum (Varoquaux et al., 2019). Biological pathways of “plant secondary cell wall biogenesis” and “plant hormone signal transduction” were enriched in Cluster C with DEGs oppositely expressed in the roots and stems, implying that expansion of the root system is essential for water exploration under drought conditions and growth of stems might be restricted by water-deficient environments. Previous research has reported similar findings that the growth of crown roots in S. italica and Maize is sensitive to drought conditions (Sebastian et al., 2016).

To achieve a better understanding of the genotypic effect on drought tolerance for all four accessions, we conducted four group comparisons between drought-tolerant and drought-sensitive accessions (Ci328 vs. Ci134, Ci328 vs. Ci603, Ci409 vs. Ci134 and Ci409 vs. Ci603) in all three tissues under both normal and drought conditions, respectively. After eliminating the DEGs that were shared by both normal and drought conditions, a total of 406 (115 in leaves, 79 in stems, and 212 in roots) drought-specific DEGs were screened out, which was defined as drought-induced DEGs between genotypes (di-DBG) (Figure 7A).

The “flavonoid biosynthesis pathway” was enriched in di-DBG identified in roots, and the “nitrogen metabolism pathway” was enriched in di-DBG detected in leaves (Figure 7B). Moreover, the water transport biological process was enriched in di-DBG identified in roots and leaves, photosynthesis, and heat response pathways, which were also enriched in di-DBG detected in leaves (Figure 7C). Furthermore, there were 12 di-DBGs overlapped with DEGs oppositely regulated between tolerant and sensitive genotypes under limited water conditions, which were identified (marked as di-DBG in Supplementary Table S7), such as HSP70 (Seita.3G216900), HSP23.5 (Seita.1G326600), and SnRK3.16 (Seita.5G031200), vital regulators of drought tolerance in plants (Cho et al., 2018; Sewelam et al., 2019; Leonetti et al., 2021).

The weighted gene co-expression network analysis (WGCNA) method constructs relationships of expressed genes identified from multi-samples through RNA-sequencing approaches. Totally, 28,401 expressed genes identified in 72 samples were divided into 34 modules (Supplementary Figure S10, Supplementary Table S11). We defined transcripts with connectivity intro-modular (KIM) higher than 0.8 as hub genes (Supplementary Table S11). There was a positive correlation of purple, yellow, midnight blue, and pink modules with the drought treatments of all accessions. There was a negative correlation of orange, red, and green modules with the drought treatment, and a positive correlation of the purple module with soluble sugar content in all accessions was identified (Supplementary Figure S11).

Among all drought-responding-related modules, we further analyzed the purple and orange modules. There was a total of 506 genes grouped into the purple module under drought conditions. We identified 17 out of 78 hub genes that are associated with drought tolerance in plant species (Figures 8A,B), such as CIPK12 (Yasuda et al., 2017), CBF1 (Wei et al., 2016), SnRK2.10 (Maszkowska et al., 2019), and RD22 (Matus et al., 2014), suggesting these gene sets might also be vital to drought tolerance of S. italica. The biological processes of “response to water” and “response to abiotic stimulus” were enriched in the purple module (Supplementary Table S12), which implied that the purple module might play an important role in foxtail millet for adjusting drought conditions.

A total of 119 genes grouped into the orange module were downregulated in leaves under drought conditions, and “monobactam biosynthesis” and “lysine biosynthesis” pathways were enriched in this module. Notably, previous studies have shown that many genes grouped into the orange module are key regulators in the stress tolerance of plants (Figures 8C,D). For instance, AREB3 (Seita.5G356200) has been identified as a downstream target of the ABA-signaling pathway and could regulate stomata openness in Arabidopsis (Qian et al., 2019). The characterization of OST1 (Seita.9G318200) in rice revealed the functions involved in regulating stomatal closure in response to water-deficient conditions through ABA-signaling modulation (Yoshida et al., 2002). Additionally, GSTF3 (Seita.5G171000) is a glutathione transferase induced by oxidative stress environments in Arabidopsis (Lee et al., 2014). Therefore, the orange module might contain vital regulators responsible for drought tolerance of foxtail millet.

Moreover, pink (734 genes, R = 0.79, P = 3e-16) and yellow (1,966 genes, R = 0.77, P = 2e-15) modules were also positively correlated with drought treatment in this study. Gene sets grouped into these two modules were mainly involved in RNA processing, transport, and the surveillance process (Supplementary Table S13), implying that post-transcriptional modulation of functional regulators might play important roles in drought tolerance of S. italica.

The light yellow module (196 genes, R = 0.57, P = 2e-7) was also positively correlated with the drought treatment in leaves (Figures 8E,F). The biological processes of “response to heat” and “very-long-chain fatty acid metabolism” were enriched in this module (Supplementary Table S12). The light yellow module contains many vital regulators identified in previous studies, such as P5C1 (Seita.5G450800), which has been confirmed to be associated with the catalysis of the final step of proline biosynthesis and involved in stress-induced proline accumulation in leaves of plants (Giberti et al., 2014). TINY (Seita.5G141600) is also essential to stress tolerance, as it binds dehydration-responsive and ethylene-responsive Cis-elements of viral regulators in Arabidopsis (Sun et al., 2008). CIPK20 (Seita.5G145900) encodes a CBL-interacting serine/threonine-protein kinase involved in the ABA-signaling pathway (Wang Z. Y. et al., 2011), and DGAT1 (Seita.3G112300) catalyzes the final step of triacylglycerol assembling, which could contribute to tolerance of freezing in Arabidopsis (Arisz et al., 2018).

Interestingly, the turquoise module was negatively correlated with the drought-treated leaves and positively correlated with drought-treated stems and roots (Supplementary Figure S11). The turquoise module was also detected as a downregulated DEG in leaves of all genotypes under drought conditions, and “photosynthesis-related biological processes” were mainly enriched in this module. This result from the turquoise module yielded more than 20 DEGs functions as components of PSI, PSII, and the light-harvesting complex (Figures 8G,H), which is consistent with decreased photosynthetic rates observed in all four accessions under drought conditions (Figure 1E). Moreover, most photosynthesis-related genes have also been downregulated in the leaves of foxtail millet under drought conditions at the jointing stage (Supplementary Figure S12), suggesting that photosynthesis-related traits change dramatically under water-limited conditions in foxtail millet.

Given that previous study about drought tolerance of S. italica mainly focused on the aboveground part of plant individuals, we have compared DEGs detected in leaves with the previously reported dataset (including 3 seedling stages and 1 germination stage) (Qi et al., 2013; Tang et al., 2017; Qin et al., 2020; Yu et al., 2020). We found 10 shared genes in response to osmotic stress during germination, seedling, and the jointing stage in foxtail millet, including HAB1 (Seita.3G164700, Lim and Lee, 2020), AREB1 (Seita.4G082000, Yoshida et al., 2010), GBF3 (Seita.5G251200, Ramegowda et al., 2017), AHG1 (Seita.6G005300, Nishimura et al., 2007), and HAI2 (Seita.9G460200, Lim et al., 2012) that participate in the ABA pathway in plants, and LEA3-4 (Seita.5G287800) induced by multi abiotic stresses (Zhao et al., 2011) and 4 more genes with unknown functions. A total of 193 jointing stage-specific DEGs were identified, including 14 DEGs involved in plant hormone regulation, 14 DEGs in response to abiotic stress, and 5 DEGs participating in photosynthesis, implying different mechanisms of drought tolerance between lateral and early growth stages in foxtail millet (Supplementary Table S14). However, due to the different genotypes tested in previous studies, and genotype-based differences could not be eliminated in the comparison analysis, further researches are still needed in the future for clarifying the unique drought-tolerant mechanisms in foxtail millet among different developmental stages.

The randomly selected DEGs in response to drought treatments from different tissues were used for qRT-PCR assays to verify the reliability of RNA-seq data. We randomly selected nine DEGs from leaves, four DEGs from stems, and six DEGs from roots, including 1 DEG (Seita.9G095900) with the opposite pattern in response to drought for further analysis. The results of qRT-PCR were highly consistent with RNA-seq data with R = 0.787, indicating that the RNA-seq results obtained in this trial are reliable (Supplementary Figure S12, Supplementary Table S15).