Transcriptomic and Metabolomic Analyses Reveal the Response to Short-Term Drought Stress in Bread Wheat ( Triticum aestivum L.)

: Drought stress, a major abiotic stress, significantly affects wheat ( Triticum aestivum L.) production globally. To identify genes and metabolic pathways crucial for responding to short-term drought stress, we conducted transcriptomic and metabolomic analyses of winter wheat cultivar Jimai 418 at four developmental stages: jointing (GS31), booting (GS45), anthesis (GS65), and 8 days after anthesis (DAA8). Transcriptomic analysis identified 14,232 differentially expressed genes (DEGs) under drought stress compared to the control. Specifically, 1387, 4573, 7380, and 892 DEGs were identified at the four developmental stages, respectively. Enriched pathways associated with these DEGs included plant hormone signal transduction, mitogen-activated protein kinase (MAPK) signaling, galactose metabolism, and starch and sucrose metabolism. Totals of 222, 633, 358, and 38 differentially accumulated metabolites (DAMs) were identified at the four stages, respectively. Correlation analysis of both datasets revealed DEGs and DAMs associated with plant hormone signal transduction, arginine and proline metabolism, ABC transporters, and amino acid biosynthesis. These findings offer significant insights into Jimai 418’s molecular response to short-term drought stress. The identified DEGs, DAMs, and enriched pathways contribute to our understanding of wheat drought tolerance. This research will facilitate further investigations into drought tolerance mechanisms and guide the breeding of wheat varieties with enhanced drought resistance.


Introduction
Wheat (Triticum aestivum L.), a crucial staple crop with global significance, plays a vital role in sustaining the world's population and driving economic growth.With a rapidly growing population, wheat production needs to significantly increase to meet rising demand [1].This challenge is further amplified by the consistent threat of drought [2].In China, the major wheat-growing regions are concentrated in the Northwest and North China Plain, with the Huang-Huai-Hai Plain particularly renowned for its wheat cultivation.The precipitation timing in this region is often inconsistent with the most critical stage of wheat growth that requires water.Therefore, insufficient rainfall is considered the primary constraint on wheat production, and precipitation shortages during wheat growth in the region exacerbate the growing food demand [3,4].
Climate change is an indisputable reality and ensuring food security for the 21st century stands as the main challenge for humanity.Global warming directly reduces crop yields [5,6].Abiotic stresses cause severe yield losses in crop production [7].Among the abiotic factors, drought significantly affects crop productivity and food safety.It can cause various detrimental effects and physiological responses, such as wilting, stomatal closure, Agronomy 2024, 14, 704 2 of 26 metabolic changes, reduced CO 2 assimilation during transpiration and photosynthesis, inhibited growth, and increased production of antioxidants.In severe cases, drought can even lead to plant death [8][9][10].Water stress exerts various effects on wheat's reproductive growth, causing significant declines in growth parameters, photosynthetic pigments, gas exchange attributes, relative water content (RWC), and yield [11].Additionally, drought stress impacts wheat variously at different growth stages, with significant yield losses observed during flowering and grain filling [2,12].
Winter wheat is the predominant cultivar in the Huang-Huai-Hai Plain, and drought has become severe in the region at the heading stage due to rising temperatures and declining precipitation.Therefore, drought is the main factor limiting crop yield increases in this region [13][14][15].Short-term drought stress is particularly critical, as winter wheat is especially susceptible to drought during critical growth stages, including sowing, jointingheading, and milky ripening/physiological maturity.The adverse impacts of drought on winter wheat during these crucial growth stages can result in yield loss [2,16].Plants adopt a variety of strategies to deal with biotic or abiotic stresses to ensure survival, and short-term stress responses are crucial for their growth and development [17].Previous research has documented yield reductions (2.03-64.39%) in winter wheat exposed to shortterm drought, particularly at different periods after the jointing period, compared to plants with sufficient water supply [12].Moreover, drought stress can impede wheat growth and reduce yield by affecting the anthesis and grain-filling processes, as well as decreasing leaf water potential, stomatal conductance, and photosynthesis [18].
Some previous studies have examined the responses of different cultivars and species to drought stress [19,20], which are mediated by a complex network of cellular and molecular processes [21].Phytohormones, such as abscisic acid (ABA) and auxin, play essential roles in abiotic stress responses [22].Several phytohormones, including ABA, auxins, cytokinins (CTKs), jasmonic acid (JA), salicylic acid, and brassinosteroids can enhance drought tolerance [23].In response to drought stress, plants activate their antioxidant defense system, upregulating enzymatic antioxidants like superoxide dismutase (SOD), catalase (CAT), and peroxidase (POD), to counteract reactive oxygen species (ROS) generation [24,25].Polyethylene glycol (PEG)-induced drought stress in wheat cultivars has shown increased activity of antioxidant enzymes, along with elevated contents of free proline (Pro), hydrogen peroxide (H 2 O 2 ), malondialdehyde (MDA), and glycine betaine (GB) [26].Proline production is regulated by pyrroline-5-carboxylate, which contributes to enhancing plant drought tolerance [27,28].Drought stress can induce complex signal transduction pathways linked to hormone metabolism.Changes in phytohormone content can lead to physiological and morphological responses in plants under various abiotic stresses.ABA and auxin are critical in response to drought stress.Under drought conditions, ABA content in plant leaves increases and stomatal conductance is inhibited, which results in a decreased transpiration rate and prevents excessive water loss [29].
Transcriptomic analyses offer valuable insights into the molecular level effects of environmental stress [4,[30][31][32].RNA sequencing (RNA-seq) and microarray technology have enabled the identification of numerous genes and pathways involved in wheat's drought response [33,34].A previous study performed transcriptomic analysis and found that DEGs associated with rehydration were enriched in pathways related to hormone metabolism and de-peroxidative stress in two different drought-resistant winter wheat cultivars [35].The high sensitivity, throughput, wide coverage, and fast separation capabilities of metabolomics analyses have led to their widespread use in plant research.Metabolomic analyses have identified key metabolites in tomato [36], cabbage [37], and crops in the Triticeae family [38][39][40].
In recent years, several studies have combined metabolomics and transcriptomics analyses to clarify the mechanisms underlying metabolite biosynthesis in plants [41][42][43][44][45][46].One study analyzed three wheat lines with varying drought tolerance using transcriptomic and metabolomic approaches under moderate and severe drought conditions and revealed that increased antioxidant activity and flavonoid content led to enhanced amino acid biosynthesis and ROS scavenging under drought stress [47].In this study, we performed transcriptomic and metabolomic analyses of droughttolerant Jimai 418 wheat at four developmental stages to identify the molecular pathways and gene expression patterns associated with its response to short-term drought stress.Our findings offer significant insights into the molecular mechanisms of wheat response to short-term drought and hold promise for breeding drought-tolerant wheat cultivars.

Plants and Experimental Site
In this study, we selected a wheat variety Jimai 418, known for its drought characteristics and yield stability [47].Field pot experiments were conducted at the experimental station of Hebei Agricultural University (Gaocheng County, Shijiazhuang, China), which is a major wheat-producing area in the Huang-Huai-Hai Plain (38.03 • N, 114.84 • E, 56 m above sea level).Each pot was initially filled with loam soil and had an equivalent weight (13.0 kg).The soil field capacity (FC, mass water content) was determined using the weighing method.To protect the plants from precipitation, a portable rain shelter was constructed.
On 12 October 2019, 30 full-sized, uniform, and healthy seeds of Jimai 418 were sown into each of the 30 metal buckets (0.9 kg; 40 cm high, bottom diameter 20 cm, top diameter 30 cm).The loam soil in the buckets was supplemented with 225 kg/ha pure nitrogen, 90 kg/ha P 2 O 5 , and 150 kg/ha K 2 O as fertilizer.

Treatments and Sampling
We selected four main developmental stages for the 5-day short-term water deficit treatments, namely GS31 (jointing, >80% visibility of the first internodes of plants), GS45 (booting, flag leaf extension, ears were not extracted), GS65 (anthesis), and DAA8 (8 days after anthesis).Each stage consisted of both a water deficit treatment and a normal control group, with 3 buckets in each group, totaling 24 buckets (Figure S1).The experiment adopted a completely randomized design.
Before each treatment period at all four stages, the soil moisture content was maintained at 60-65% FC.Five days prior to each sampling, the deficit treatment group's soil water content was reduced and maintained at 45-50% FC, while the control group was increased and maintained at 75-80% FC.
Total weight measurements were taken at each stage.Plant weight measurements were ignored; the bucket weight (0.9 kg) and dry soil weight were subtracted from the total weight, and the amount of water required for watering or dispersion was calculated according to the soil FC of the treatment (Table 1; Figure S2).The growth stages were documented using the Zadoks scale [48].The upper expanded leaves from five individual plants per bucket were collected between 10 and 11 a.m., and these samples were pooled as one replicate (Figure S3).The samples were frozen in liquid nitrogen and preserved at −80 • C. For transcriptome analysis, three biological replications were used, while six biological replications for metabolome analysis.Soil moisture content of all groups was maintained at 75-80% FC after sampling.Leaf areas (length × width × 0.75) were measured 15 days post-flowering.The mature grains were harvested on 31 May 2020.At the fully mature stage, the above-ground dry matters of 20 plants per pot were collected, then the total stems were counted to calculate the spike per plant, and air-dried to constant weight; the total weight was the biomass.Grains were then separated by hand, counted, and 1000-grain weight was calculated.

Transcriptome Sequencing and Functional Annotation
Total RNA samples were collected utilizing TRIzol reagent (Invitrogen, Carlsbad, CA, USA).The RNA concentration and integrity were measured using a Qubit fluoromete and polyacrylamide gel electrophoresis, respectively.RNA-seq was conducted on the Illumina NovaSeq platform with paired-end reads of 150 bp.Low-quality reads (i.e., reads with Q < 30) were filtered and excluded.Clean reads were then aligned to Chinese Spring wheat reference genome (v2.0) with a maximum of two alignments allowed.HTSeq (v0.6.1) was employed to determine the fragment per kilobase million (FPKM) values for each gene [53].The DESeq package was employed to identify DEGs with the criteria padj < 0.01 and |log 2 foldchange| > 2 [28].Enrichment analysis of GO and KEGG pathways was conducted utilizing R package 'clusterProfiler' (3.4.4) with a statistical significance threshold of less than 0.05.

Metabolome Analysis
Metabolites were extracted using a previously described protocol [54].The identification of untargeted metabolites was conducted using UHPLC-MS/MS (ultra-performance liquid chromatography and tandem mass spectrometry) on a Vanquish UHPLC coupled with an Orbitrap Q Exactive TM HF mass spectrometer (Thermo Fisher, Dreieich, Germany) at Beijing Novogene Biotechnology Co., Ltd.(Beijing, China).The analytical conditions were consistent with those described in a previous study [55].
For the metabolites, KEGG pathway analysis was performed.Statistical significance of differences was determined using orthogonal partial least squares discriminant analysis (OPLS-DA), with p-value below 0.05 considered statistical significance.Differentially accumulated metabolites (DAMs) were identified using p < 0.05, VIP (variable importance in projection) > 1, and FC ≥ 2 or ≤0.5 as the criteria.

Quantitative Real-Time PCR (qRT-PCR) Validation
Total RNA was extracted from all 24 samples and quantified using a DS-11 Spectrophotometer (DeNovix Inc., Wilmington, DE, USA).EasyScript All-in-One First-Strand cDNA Synthesis SuperMix for qRT-PCR (TransGen Biotech, Beijing, China) was employed to reverse-transcribe 1 mg of total RNA into first-strand cDNA with oligo dT primers.The resulting cDNA templates were preserved at −20 • C.
qRT-PCR was conducted following established protocols [35], with three replicates per sample.The reaction involved an initial denaturation step (95 • C, 3 min), followed by 40 amplification cycles (95 • C for 10 s and 60 • C for 30 s).The 2 −△△Ct method was utilized to determine relative gene expression levels [35], with wheat β-actin gene as the internal control.Primer sequences are listed in Table S1.

Statistical Analysis
All measurements and analyses were conducted on at least 3 replicates.Statistical analyses and figure generation were conducted using Microsoft Excel 2010 and SPSS v18.0.Data were analyzed using ANOVA tests between control and drought treatments at different growth stages.The mapping data were expressed as mean ± standard deviation (SD), with a significance level of p < 0.05.Correlations between DEGs and DAMs were evaluated utilizing Pearson's correlation coefficient (p < 0.05).Principal component analysis (PCA) and OPLS-DA analysis were conducted utilizing the OmicShare tool (https://www.omicshare.com/tools,accessed on 8 June 2022).

Transcriptome Sequences and Differential Expression Analysis
RNA sequencing of the 24 leaf samples (4 stages × 2 treatments × 3 replications) using Illumina paired-end sequencing generated over 201 million clean reads, with an approximate GC content of 50%.Each library yielded a total of 78-89 million reads.Mapping of these clean reads to the wheat reference genome sequences revealed a success rate of 88.66% to 90.05% (Table S2).PCA was conducted on the data of plants at different stages and under different treatments (Figure 2a).The first and second principal component (PC1 and PC2) accounted for 28.80% and 26.34% of total variance, respectively.Pearson's correlation coefficients of the gene expression data were calculated, and R 2 values were greater than 0.89 (Figure S4).Thus, the gene expression level for each sample was determined by averaging FPKM values from its three replicates.

Transcriptome Sequences and Differential Expression Analysis
RNA sequencing of the 24 leaf samples (4 stages × 2 treatments × 3 replications) using Illumina paired-end sequencing generated over 201 million clean reads, with an approximate GC content of 50%.Each library yielded a total of 78-89 million reads.Mapping of these clean reads to the wheat reference genome sequences revealed a success rate of 88.66% to 90.05% (Table S2).PCA was conducted on the data of plants at different stages and under different treatments (Figure 2a).The first and second principal component (PC1 and PC2) accounted for 28.80% and 26.34% of total variance, respectively.Pearson's correlation coefficients of the gene expression data were calculated, and R 2 values were greater than 0.89 (Figure S4).Thus, the gene expression level for each sample was determined by averaging FPKM values from its three replicates.
DEGs were identified by comparing the short-term drought stress samples with their respective controls at each developmental stage (Figure 2b).In the GS31_D (short-term drought stress) vs. GS31_C (the control, non-stress) group, 1387 DEGs were identified, with 838 DEGs showing upregulation and 549 DEGs showing downregulation.In the GS45_D vs. GS45_C group, there were 2934 upregulated and 1639 downregulated DEGs.For the GS65_D vs. GS65_C group, 3515 DEGs were upregulated and 3865 were DEGs were identified by comparing the short-term drought stress samples with their respective controls at each developmental stage (Figure 2b).In the GS31_D (short-term drought stress) vs. GS31_C (the control, non-stress) group, 1387 DEGs were identified, with 838 DEGs showing upregulation and 549 DEGs showing downregulation.In the GS45_D vs. GS45_C group, there were 2934 upregulated and 1639 downregulated DEGs.For the GS65_D vs. GS65_C group, 3515 DEGs were upregulated and 3865 were downregulated.For the DAA8_D vs. DAA8_C group, there were 609 upregulated and 283 downregulated DEGs.Venn diagram analysis (Figure 2c) revealed only 49 DEGs were shared among the four groups, suggesting substantial regulatory differences in gene expression across the various developmental stages.
Excluding the common DEGs, 10,317 DEGs were identified in the short-term drought stress vs. control comparison group at the four developmental stages (Figure 3; Table S3).Hierarchical clustering analysis revealed significant changes in unigene expression across all samples, and this revealed variation in the transcriptome among developmental stages.Drought had the most substantial impact on the expression of DEGs in GS45 and GS65, which was the period of endosperm cell development.The expression of DEGs at GS31 differed from those at GS45 and GS65, and the expression of DEGs at DAA8 was similar under short-term drought stress and control conditions (Figure 3).Hierarchical clustering analysis revealed significant changes in unigene expression across all samples, and this revealed variation in the transcriptome among developmental stages.Drought had the most substantial impact on the expression of DEGs in GS45 and GS65, which was the period of endosperm cell development.The expression of DEGs at GS31 differed from those at GS45 and GS65, and the expression of DEGs at DAA8 was similar under short-term drought stress and control conditions (Figure 3).

GO and KEGG Enrichment Analysis of DEGs in Response to Short-Term Drought at Four Developmental Stages
GO analysis was performed to investigate the functional roles of DEGs in response to short-term drought stress (Table S4).For the GS31_D vs. GS31_C comparison, the most enriched GO terms were associated with embryo development, phosphoric ester hydrolase activity, oxidoreductase activity, protein disulfide oxidoreductase activity, and stress-responsive activities that response to water, response to abiotic stimulus, and response to acid chemical.For GS45_D vs. GS45_C, significantly enriched GO terms included drought stress-responsive activities (e.g., response to water, response to acid, and response to abiotic stimulus), sugar metabolism and oxidoreductase activity (e.g., catalase activity, histone acetyltransferase activity, and peptide N-acetyltransferase activity).For GS65_D vs. GS65_C, the top enriched GO terms were the same as those in GS45_D vs. GS45_C, with the exceptions of photosynthesis-antenna proteins, phenylpropanoid biosynthesis, cyanoamino acid metabolism, glycerophospholipid metabolism, and fatty acid elongation, which were significantly enriched in GS65_D vs. GS65_C but not in GS45_D vs. GS45_C.For DAA8_D vs. DAA8_C, the molecular function category exhibited the highest enrichment of GO terms.
To reveal the biological pathways associated with short-term drought stress, we performed KEGG enrichment analysis (Figure 4).In the GS31_D vs. GS31_C group, the most enriched pathways included: plant hormone signal transduction; ribosome; MAPK signaling pathway; galactose metabolism; and starch and sucrose metabolism pathways (Figure 4a).In the GS45_D vs. GS45_C group, significantly enriched pathways included glyoxylate and dicarboxylate metabolism, starch and sucrose metabolism, galactose metabolism, plant hormone signal transduction, and MAPK signaling pathways (Figure 4b).For the GS65_D vs. GS65_C group, we observed significant enrichment in pathways including plant hormone signal transduction, photosynthesis-antenna proteins, glyoxylate and dicarboxylate metabolism, starch and sucrose metabolism, phenylpropanoid biosynthesis, cyanoamino acid metabolism, the MAPK signaling pathway, glycerophospholipid metabolism, galactose metabolism, and fatty acid elongation (Figure 4c).For the DAA8_D vs. DAA8_C group, many DEGs exhibited enrichment in the following pathways: zeatin biosynthesis; cutin, suberin, and wax biosynthesis; glycerophospholipid metabolism; plant hormone signal transduction; and arginine and proline metabolism (Figure 4d).The consistent enrichment of plant hormone signal transduction and the MAPK signaling pathways at stages GS31, GS45, GS65, and DAA8 suggests their potential crucial roles in wheat's response to drought stress.
We generated a heatmap to visualize the expression patterns of the 49 common DEGs among the four comparison groups (Figure 5a; Table S5).Out of these 49 DEGs, 38 (77.6%) exhibited upregulated expression (Table S5).GO analysis revealed these common DEGs were mainly enriched in plant response-related terms, including response to abiotic stimulus, response to acid chemical, and response to water (Figure 5b).KEGG enrichment analysis highlighted the notable enrichment of these DEGs in the plant hormone signal transduction and MAPK signaling pathways (Figure 5c).

Metabolomic Response of Wheat to Short-Term Drought Stress
Metabolites were identified in the two treatments and four developmental stages, with six replicates conducted for each sample.OPLS-DA revealed that PC1 correlated with developmental stage accounted for 14.4% of total variation, while PC2 correlated with drought treatment accounted for 21.6% of total variation.These results indicate that most of the variation in metabolites stemmed from developmental stage and drought treatment, and notable metabolite differences were observed between treatments (Figure 6a).In total, 1251 DAMs were identified in the 48 samples.In the GS31_D vs. GS31_C group, 165 downregulated and 57 upregulated DAMs were detected.In the GS45_D vs. GS45_C group, 363 downregulated and 270 upregulated DAMs were identified.A total of 358 DAMs, including 206 downregulated and 152 upregulated DAMs, were identified in the GS65_D vs. GS65_C group.A total of 38 DAMs, including 24 downregulated and 14 upregulated DAMs, were identified in the DAA8_D vs. DAA8_C group (Figure 6b, Table S6).
The abundances of metabolites differed significantly among the four comparison groups.The Venn diagram revealed 103, 336, 92, and 13 DAMs in the GS31_D vs. GS31_C, GS45_D vs. GS45_C, GS65_D vs. GS65_C, and DAA8_D vs. DAA8_C comparison groups, respectively (Figure 6c).Five metabolites (two increased and three decreased in abundance) were common to all four comparison groups (Table S7).Specifically, the concentration of N-(2,5-dioxo-4-imidazolidinyl) urea (allantoic acid) was 5.75-, 56.54-, and 37.80-fold higher in GS31, GS45, and GS65 samples than in control samples, respectively.The DLarginine concentration was 3.89-and 5.98-fold higher in GS45 and GS65 samples than in control samples, respectively (Table S7).These findings indicate that the metabolic pathways mediating drought tolerance in plants vary among stages.Further analysis through KEGG pathway enrichment revealed that these DAMs were primarily enriched in several metabolic processes, including amino acid metabolism, biosynthesis of other secondary metabolites, metabolism of cofactors and vitamins, carbohydrate metabolism, lipid metabolism, and nucleotide metabolism (Figure 6d).exhibited upregulated expression (Table S5).GO analysis revealed these commo were mainly enriched in plant response-related terms, including response to stimulus, response to acid chemical, and response to water (Figure 5b).KEGG enri analysis highlighted the notable enrichment of these DEGs in the plant hormon transduction and MAPK signaling pathways (Figure 5c).concentration of N-(2,5-dioxo-4-imidazolidinyl) urea (allantoic acid) was 5.75-, 56.54-, and 37.80-fold higher in GS31, GS45, and GS65 samples than in control samples, respectively.The DL-arginine concentration was 3.89-and 5.98-fold higher in GS45 and GS65 samples than in control samples, respectively (Table S7).These findings indicate that the metabolic pathways mediating drought tolerance in plants vary among stages.Further analysis through KEGG pathway enrichment revealed that these DAMs were primarily enriched in several metabolic processes, including amino acid metabolism, biosynthesis of other secondary metabolites, metabolism of cofactors and vitamins, carbohydrate metabolism, lipid metabolism, and nucleotide metabolism (Figure 6d).

Combined Transcriptomic and Metabolomic Analysis
A combined transcriptomic and metabolomic analysis was conducted to clarify the regulation of metabolic network.A total of 29, 42, 49, and 13 co-enriched pathways were identified at GS31, GS45, GS65, and DAA8, respectively.Three co-enriched pathways were detected in all four developmental stages, namely galactose metabolism, pentose and glucuronate interconversions, and beta-alanine metabolism (Figure 7).A total of 35 co-enriched pathways were detected in both GS45 and GS65, including membrane transport, global and overview maps, carbohydrate metabolism, biosynthesis of other secondary metabolites, energy metabolism, lipid metabolism, nucleotide metabolism, amino acid metabolism, metabolism of other amino acids, metabolism of cofactors and vitamins and translation.These shared pathways are likely to play a crucial role in conferring drought tolerance to wheat.Most DEGs were enriched within these co-enriched pathways, exhibiting expression patterns that highly aligned with the corresponding metabolite levels.Additionally, several commonly enriched pathways were observed, such as plant hormone signal transduction, biosynthesis of amino acids, arginine and proline metabolism, and ATP-binding cassette (ABC) transporters.

Combined Transcriptomic and Metabolomic Analysis
A combined transcriptomic and metabolomic analysis was conducted to clarify the regulation of metabolic network.A total of 29, 42, 49, and 13 co-enriched pathways were identified at GS31, GS45, GS65, and DAA8, respectively.Three co-enriched pathways were detected in all four developmental stages, namely galactose metabolism, pentose and glucuronate interconversions, and beta-alanine metabolism (Figure 7).A total of 35 co-enriched pathways were detected in both GS45 and GS65, including membrane transport, global and overview maps, carbohydrate metabolism, biosynthesis of other secondary metabolites, energy metabolism, lipid metabolism, nucleotide metabolism, amino acid metabolism, metabolism of other amino acids, metabolism of cofactors and vitamins and translation.These shared pathways are likely to play a crucial role in conferring drought tolerance to wheat.Most DEGs were enriched within these coenriched pathways, exhibiting expression patterns that highly aligned with the corresponding metabolite levels.Additionally, several commonly enriched pathways were observed, such as plant hormone signal transduction, biosynthesis of amino acids, arginine and proline metabolism, and ATP-binding cassette (ABC) transporters.

DEGs Involved in Plant Hormone Signal Transduction Network in Response to Short-Term Drought
A total of 30 DEGs associated with plant hormone signal transduction were identified.These DEGs mainly participated in the synthesis and decomposition of ABA, auxin, gibberellin, CTK, and brassinosteroid (Figure 8a; Table S8).The expression of PYL gene (TraesCS2B02G105300) encoding ABA receptor PYR/PYL family members was downregulated at all four developmental stages.Genes encoding protein phosphatase 2C (PP2C), including TraesCS4A02G094300, TraesCS5D02G188600, TraesCS1B02G441400, TraesCS3D02G212100, TraesCS3B02G394600, and TraesCS2A02G000400 exhibited upregulated expression in each of the four comparison groups.The expression levels of three sucrose non-fermenting-1-related protein kinase 2 genes (SnRK2, TraesCS2B02G521800, TraesCS1A02G215900, and TraesCS1D02G271000) and two ABAresponsive element binding factor (ABF) genes (TraesCS6B02G364000, and TraesCS3D02G371900) were upregulated in all comparison groups (Table S8; Figure 8a).Metabolomic analysis revealed that ABA contents were 5.74-, 12.51-, and 9.94-fold higher under drought treatment than under control conditions in the GS31_D vs. GS31_C, GS45_D vs. GS45_C, and GS65_D vs. GS65_C comparison groups, respectively (Figure 8b).

DEGs Involved in Plant Hormone Signal Transduction Network in Response to Short-Term Drought
A total of 30 DEGs associated with plant hormone signal transduction were identified.These DEGs mainly participated in the synthesis and decomposition of ABA, auxin, gibberellin, CTK, and brassinosteroid (Figure 8a; Table S8).The expression of PYL gene (TraesCS2B02G105300) encoding ABA receptor PYR/PYL family members was downregulated at all four developmental stages.Genes encoding protein phosphatase 2C (PP2C), including TraesCS4A02G094300, TraesCS5D02G188600, TraesCS1B02G441400, TraesCS3D02G212100, TraesCS3B02G394600, and TraesCS2A02G000400 exhibited upregulated expression in each of the four comparison groups.The expression levels of three sucrose non-fermenting-1-related protein kinase 2 genes (SnRK2, TraesCS2B02G521800, TraesCS1A02G215900, and TraesCS1D02G271000) and two ABA-responsive element binding factor (ABF) genes (TraesCS6B02G364000, and TraesCS3D02G371900) were upregulated in all comparison groups (Table S8; Figure 8a).Metabolomic analysis revealed that ABA contents were 5.74-, 12.51-, and 9.94-fold higher under drought treatment than under control conditions in the GS31_D vs. GS31_C, GS45_D vs. GS45_C, and GS65_D vs. GS65_C comparison groups, respectively (Figure 8b).

Combined Transcriptomic and Metabolomic Analyses of the ABC Transporter Pathway and Arginine and Proline Metabolism Pathway
Correlation analysis revealed a notable enrichment of DEGs and DAMs within the ABC transporter pathway at GS65 (Figures S5 and S6).Six DEGs exhibited significant correlation with nine DAMs, namely, L-cystine, L-lysine, melibiose, inositol, L-phenylalanine, thiamine, biotin, choline, and norfloxacin (Figure 10a and Figure S6; Table S9).A heatmap showed that the abundances of these nine metabolites varied among the different groups (Figure 10b).L-lysine, inositol, and L-phenylalanine exhibited high abundances at GS45.The abundances of inositol and L-phenylalanine were high at GS65, while melibiose exhibited high abundance at DAA8.To confirm gene expression patterns identified through RNA-seq, six genes (novel.13863,TraesCS3A02G246000, TraesCS1A02G343300, TraesCS3D02G371900 TraesCS7D02G204300, and TraesCS1B02G441400) were chosen randomly from DEGs (Figure 9a).qRT-PCR data for these six genes exhibited a high correlation with the corresponding data from RNA-seq (R 2 = 0.8246, Figure 9b).
exhibited high abundance at DAA8.
The abundance of D-proline was 4.43-, 104.46-, and 6.94-fold higher under drought treatment than under control conditions in the GS31_D vs. GS31_C, GS45_D vs. GS45_C, and GS65_D vs. GS65_C comparison groups, respectively.4-Oxoproline abundance was 0.38-and 2.80-fold higher under drought treatment than under control conditions in the GS31_D vs. GS31_C and GS45_D vs. GS45_C comparison groups.5-Aminovaleric acid abundance was 19.22-and 14.56-fold higher under drought treatment than under control conditions in the GS45_D vs. GS45_C and GS65_D vs. GS65_C groups, respectively (Figure 11a).Increased proline synthesis was induced by glutamate 5-kinase (proB, 2.7.2.11), ornithine-oxo-acid transaminase (rocD, OAT, 2.6.1.13),arginase (rocF, ARG, 3.5.3.11), and proline iminopeptidase (pip, 3.5.3.11)mainly via the arginine and proline metabolism pathways (Figure 11a).RNA-seq analysis revealed that P5CS expression (TraesCS1D02G280700, 2.89-fold change) was notably upregulated in the GS65_D vs. GS65_C comparison group.The expression of TraesCS3B02G538100 increased 2.34-and 2.32-fold, and the expression of TraesCS3D02G483400 encoding P5CR increased 2.09-and 2.59-fold in the GS45_D vs. GS45_C and GS65_D vs. GS65_C groups, respectively.The expression of TraesCS5D02G384500 encoding ornithine aminotransferase (δ-OAT) was 2.55-and 2.76-fold higher under short-term drought stress than under control conditions in the GS45_D vs. GS45_C and GS65_D vs. GS65_C groups, respectively.The expression of TraesCS2A02G035700 encoding rocF, also called arg, was 2.30-fold higher under short-term drought stress than under control conditions in the GS65_D vs. GS65_C group.The expression of TraesCS1D02G336300 increased 2.72-, 3.47-, and 2.79-fold under drought stress than under control conditions in the GS45_D vs. GS45_C, GS65_D vs. GS65_C, and DAA8_D vs. DAA8_C groups, respectively.The expression of novel.2587encoding pip increased 2.49-and 2.60-fold under drought stress than under control conditions in the GS45_D vs. GS45_C and GS65_D vs. GS65_C groups, respectively.The expression of ProDH (TraesCS1A02G209100) significantly decreased over the experimental period, and these genes were mainly associated with the proline degradation pathway (Figure 11b).A correlation analysis of the 49 common genes and five common metabolites in the four comparison groups was conducted (Figure S7).Further analyses were conducted on DEGs and metabolites based on a Pearson's correlation coefficient greater than 0.8 or less than −0.8 as the selection criteria.A total of 29 genes were positively correlated with allantoic acid, 10 genes were positively correlated with DL-arginine, TraesCS2B02G245900 was positively correlated with maltopentaose, and TraesCS3A02G055500 exhibited a positive correlation with melibiose (Table S10).TraesCS1D02G237500 and TraesCS1B02G249000 displayed upregulation in the GS31_D vs. GS31_C, GS45_D vs. GS45_C, and GS65_D vs. GS65_C groups, but downregulated in the DAA8_D vs. DAA8_C group.Allantoic acid was not a DAM in the DAA8_D vs. DAA8_C group.

Discussion
Frequent drought events can have deleterious effects on winter wheat yield [56].In the Huang-Huai-Hai Plain, water scarcity is the main factor limiting crop production, and growing demand for food in this region has further exacerbated the effects of water shortages [3].Annual precipitation in this region is low and unevenly distributed.Drought stress can affect wheat growth at various developmental stages, with differing impacts during the sowing, jointing-heading, and milky ripening/physiological maturity stages.Plants undergo various physiological and biochemical changes to adapt to drought stress [57].Drought can have deleterious effects on wheat yield at various developmental stages, for example biomass yield, population, leaf area, grain number, and 1000-grain weight [58].In this study, we investigated short-term drought responses in wheat variety Jimai 418 at four developmental stages.Drought during the reproductive stage leads to spikelet sterility in wheat.The 1000-grain weight was lowest in plants at GS65_D, showing a 5.8% reduction compared to GS65_C (Table 2).Leaf area was also notably smaller under short-term drought stress than under control conditions (Table 2).Previous studies have also reported a negative impact of drought stress during anthesis on grain size, in the most drought-tolerant cultivars such as Halberd and SYN604 [58].Thus, drought at anthesis had a negative effect on yield.The effects of drought on yield were weaker in late stages compared to early stages.These findings align with previous studies of wheat, demonstrating that sufficient water supply can promote tiller growth and panicle formation, drought during GS31 causes floret and whole-spikelet death, and drought during grain filling reduces grain size and weight [58,59].
The response of the Jimai 418 cultivar to drought stress varied among developmental stages, as determined by the physiological parameters measured in this study.POD acts as a detoxifying enzyme and removes toxic reductants during plant responses to abiotic stress [22].Our findings demonstrated that POD activity increased in response to short-term drought treatment at all four developmental stages (Figure 1).Cell membrane integrity and stability are essential for drought tolerance under water stress conditions [57].Increased levels of MDA under stress conditions suggest that drought stress may induce membrane lipid peroxidation via over production of ROS [60,61].For this reason, MDA is frequently used to assess drought tolerance [59].Our study revealed significantly higher MDA content under short-term drought stress compared to control at GS65 and DAA8 (Figure 1).Based on these observations, drought stress resulted in physiological changes at GS45, GS65, and DAA8.Irrigation regimes can be optimized through an understanding of crop development stage responses to water deficit stress [62].In our study, total plant biomass and yield were assessed to evaluate the influence of short-term drought on wheat growth and development.Drought had the greatest impact on yield at GS31.
Our study combined transcriptomic and metabolic analyses to gain valuable insights into how wheat responds to short-term drought stress.Transcriptomic analyses revealed substantial differences in gene expression among the four developmental stages.KEGG pathway analysis revealed the activation of plant hormone signal transduction and MAPK signaling pathways at GS31, GS45, GS65, and DAA8, suggesting their conserved and important roles in wheat's response to short-term drought stress.Plant hormones are known to trigger a series of responses to drought stress [63,64].Plant hormone signal transduction plays a critical role in both roots and leaves of Populus euphratica to short-term drought stress [17].Our study identified 30 DEGs involved in plant hormone signaling pathways, with 12, 12, 4, and 2 DEGs associated with the ABA, auxin, CTK, and JA pathways, respectively.Endogenous ABA promotes drought tolerance by activating stressrelated gene expression [65].Our findings revealed 12 ABA-related genes, including one PYR/PYL, six PP2Cs, three SnRK2s, and two ABFs, were differentially expressed at the four stages, suggesting the pivotal roles of ABA-related genes in Jimai 418's drought response.ABA binds to PYP/PYL and promotes interaction with PP2C, which inhibits PP2C activity and activates the protein kinase SnRK2 [66].Activated SnRK2 phosphorylates downstream TFs, which leads to an ABA signal response [67].Our metabolic profiling revealed ABA content was 5.74-, 12.51-, and 9.94-fold higher at GS31, GS45, and GS65 under short-term drought treatment than under control, respectively; however, no differences in ABA content under short-term drought and control conditions were observed at DAA8.In our study, PP2C expression was upregulated in response to short-term drought.This might stem from the fact that ABA induces PP2C expression, which regulates SnRK2 activity and induces ABA signaling.In a future study, we plan to evaluate the efficacy of applying exogenous ABA application in different drought-tolerant wheat cultivars.
IAA is related to drought tolerance and has positive effects on plants under drought stress [68].The application of IAA can alleviate PEG stress's adverse effects and enhance barley growth [69].IAA can mitigate increases in ABA levels via phosphate-related degradation through the MAPK signaling pathway [70].Our findings are in line with these conclusions, further emphasizing the essential role plant hormones play in Jimai 418's response to short-term drought stress.
We also observed the impact of metabolites, including proline, on wheat's response to short-term drought stress.Through RNA-seq analysis, we identified key droughtresponsive genes that play important roles in mediating the response to short-term drought stress in wheat.Through a targeted gas chromatography-mass spectrometry approach, Bowne et al. examined three cultivars of bread wheat with varying levels of drought tolerance under drought stress, and found that levels of proline, tryptophan, and branchedchain amino acids increased in all cultivars under drought stress [71].In our study, D-proline content was 4.43-, 104.46-, and 6.94-fold higher in GS31, GS45, and GS65 after drought stress treatment, respectively.OAT, a crucial enzyme in the ornithine pathway of proline biosynthesis, plays a significant role in multiple conditions [72].Overexpression of the OsOAT gene in rice notably enhanced both drought tolerance and osmotic stress tolerance [73].The DEGs associated with the arginine and proline metabolism pathway, including key enzymes (P5CS, PDH, OAT, and ARG1) [74], were upregulated in DAA8, GS31, GS45, and GS65, respectively (Figure 11a), and metabolite contents were similar to the gene expression.
ABC transporters are one of the largest protein families, functioning as both exporters and importers [75].They play key roles in the transmembrane transport of a variety of molecules in plants to mediate adaptation to rapidly changing adverse environments [76].Cells require the absorption of various nutrient elements, the discharge of endogenous toxins, and the exchange of signaling molecules to survive changes in abiotic conditions [76].Thus, ABC transporters have various functions, and the regulation of stress responses is complex.Our findings revealed a notable enrichment of DEGs and DAMs in ABC transporters within the leaves of wheat Jimai 418.According to the heatmap analysis results, a total of six DEGs exhibited a strong correlation with nine metabolites clustered into different groups (Figure 10).These ABC transporters may play important roles in wheat's response to short-term drought stress.Overall, the response to short-term drought stress involves multiple pathways, including those related to ABC transporters, amino acid metabolism, and plant hormone signaling, and common interactions enhanced wheat's drought tolerance.

Conclusions
In this study, transcriptome sequencing and metabolite profiling were conducted on the drought-tolerant wheat variety Jimai 418 at four developmental stages.In total, 14,232 DEGs and 1251 metabolites were identified between the four developmental stages.Comprehensive analysis of transcriptomic and metabolomic data uncovered the prominent involvement of pathways such as plant hormone signal transduction, ABC transporter, and arginine and proline metabolism in the response to drought stress.These findings provide a comprehen- drought treatments at the same stages based on a two-way ANOVA test.*, **, and *** indicate significance at the 5% (p < 0.05), 1% (p < 0.01), and 0.1% (p < 0.001) level, respectively.

Figure 1 .
Figure 1.Changes in physiological parameters, including MDA, POD, IAA, and ABA contents, in wheat plants at four growth stages under short-term drought stress.MDA, malondialdehyde; POD, peroxidase; IAA, indole-3-acetic acid; ABA, abscisic acid.Values represent mean ± standard deviation (SD) from three biological replicates.Different lowercase letters within a column indicate significant differences (p < 0.05) between control and drought treatments at the same growth stage according to a one-way ANOVA test.

Figure 1 .
Figure 1.Changes in physiological parameters, including MDA, POD, IAA, and ABA contents, in wheat plants at four growth stages under short-term drought stress.MDA, malondialdehyde; POD, peroxidase; IAA, indole-3-acetic acid; ABA, abscisic acid.Values represent mean ± standard deviation (SD) from three biological replicates.Different lowercase letters within a column indicate significant differences (p < 0.05) between control and drought treatments at the same growth stage according to a one-way ANOVA test.

Figure 2 .
Figure 2. Inter-sample analysis of Jimai 418 leaf transcriptome data.(a) Principal component analysis (PCA).There are relative coordinate points on the principal component after the samples

Figure 2 .
Figure 2. Inter-sample analysis of Jimai 418 leaf transcriptome data.(a) Principal component analysis (PCA).There are relative coordinate points on the principal component after the samples are analyzed by dimension reduction.Closer proximity between points indicates higher similarity between samples.The dots of different colors in the figure represent different stages and treatments.(b) Number of differentially expressed genes (DEGs) that were upregulated and downregulated at the stages GS31, GS45, GS65, and DAA8.(c) Venn diagram showing the overlap of DEGs among the indicated comparison groups.

Figure 3 .
Figure 3. Heatmap of DEGs accumulation patterns in comparison groups of Jimai 418.Log2(FPKM + 1) values were extracted from DEGs expression profile and visualized using hierarchical

Figure 3 .
Figure 3. Heatmap of DEGs accumulation patterns in comparison groups of Jimai 418.Log 2 (FPKM + 1) values were extracted from DEGs expression profile and visualized using hierarchical clustering.The normalized gene expression levels are represented by a color gradient from red (high expression) to white (intermediate expression) to blue (low expression).

Figure 5 .
Figure 5. GO classification and KEGG enrichment analysis of 49 common DEGs across all comparison groups.(a) Heatmap of accumulation patterns of 49 common DEGs in the comparison groups.(b) Significantly enriched GO terms (corrected p value ≤ 0.05).(c) KEGG pathway enrichment scatter diagram.Dot colors represent the correlated q value, while dot sizes represent the number of DEGs.

Agronomy 2024 ,
14, x FOR PEER REVIEW 16 of 29 Metabolome analysis included six independent biological replicates.(c) Venn diagram of DAMs across the four comparisons.(d) KEGG pathway analysis of DAMs identified in the four comparison groups.

Figure 7 .
Figure 7. KEGG pathways of DEGs and DAMs co-enrichment at four comparison groups based on joint analysis of transcriptomics and metabolomics.X-axis (ratio): the proportion of differential genes or metabolites enriched in each pathway relative to the total number of genes or metabolites annotated to that pathway.Y-axis: KEGG pathways enriched by both metabolites and transcriptomics.(a) GS31, (b) GS45, (c) GS65, and (d) DAA8.

Figure 7 .
Figure 7. KEGG pathways of DEGs and DAMs co-enrichment at four comparison groups based on joint analysis of transcriptomics and metabolomics.X-axis (ratio): the proportion of differential genes or metabolites enriched in each pathway relative to the total number of genes or metabolites annotated to that pathway.Y-axis: KEGG pathways enriched by both metabolites and transcriptomics.(a) GS31, (b) GS45, (c) GS65, and (d) DAA8.

Figure 8 .
Figure 8.(a) Heatmap of accumulation patterns of DEGs in plant hormone signal transduction pathway across comparison groups.The normalized gene expression levels are represented by a color gradient from red (high expression) to white (intermediate expression) to blue (low expression).(b) ABA content at four developmental stages.

Figure 8 .
Figure 8.(a) Heatmap of accumulation patterns of DEGs in plant hormone signal transduction pathway across comparison groups.The normalized gene expression levels are represented by a color gradient from red (high expression) to white (intermediate expression) to blue (low expression).(b) ABA content at four developmental stages.

Figure 9 .
Figure 9. (a) Expression levels of six randomly selected DEGs as determined by qRT-PCR and RNAseq data.(b) Correlation analysis on the overall expression levels of six selected DEGs between RNA-seq and qRT-PCR.Log2FC values of RNA-seq data (x-axis) were plotted against log2FC values of qRT-PCR data (y-axis).

Figure 9 .
Figure 9. (a) Expression levels of six randomly selected DEGs as determined by qRT-PCR and RNA-seq data.(b) Correlation analysis on the overall expression levels of six selected DEGs between RNA-seq and qRT-PCR.Log2FC values of RNA-seq data (x-axis) were plotted against log2FC values of qRT-PCR data (y-axis).

Figure 10 .
Figure 10.Heatmap of Key DEGs (a) and DAMs (b) in the ABC transporter pathway.Figure 10.Heatmap of Key DEGs (a) and DAMs (b) in the ABC transporter pathway.

Figure 10 .
Figure 10.Heatmap of Key DEGs (a) and DAMs (b) in the ABC transporter pathway.Figure 10.Heatmap of Key DEGs (a) and DAMs (b) in the ABC transporter pathway.

Table 1 .
The calculation according to the soil field capacity of treatment.

Table 2 .
Impact of short-term drought stress on leaf areas, biomass, yield, grain number, 1000-grain weight, and spike number at different growth stages.