Wheat (Triticum aestivum) production is increasingly challenged by drought and heat stresses, which often coincide under climate change, leading to severe yield losses. Drought and heat are among the most damaging abiotic threats to wheat, each capable of causing drastic yield reductions (on the order of ~ 40–60% in severe cases) and even greater losses when they occur together [1,2]. These stresses frequently occur simultaneously in many wheat-growing regions, as rising global temperatures and shifting rainfall patterns increase the likelihood of concurrent heatwaves and droughts. Indeed, recent genomic studies in wheat underscore that studying heat or drought separately offers limited gains, whereas targeting combined heat–drought resilience can reveal unique tolerance loci [3]. This highlights an urgent need for innovative strategies to develop climate-resilient wheat varieties that can maintain productivity under simultaneous drought and heat stress.

Stress memory is an emerging concept in plant biology that holds promise for enhancing plant tolerance to adverse conditions. Stress memory refers to the ability of plants to “remember” a prior exposure to stress and to mount faster or stronger responses upon subsequent stress encounters. In essence, a mild or initial stress episode can prime the plant’s physiological and molecular defenses, conferring improved tolerance to later stresses of the same or even different type. This phenomenon has been documented in various plant species and stress scenarios. For example, a preliminary drought event can enhance photoprotective capacity and yield in perennial grasses during a later drought, and heat priming at vegetative stages has been shown to improve thermotolerance during grain filling in rice [4]. Similar priming benefits are reported in other crops, such as radish and tobacco, as well as in wheat itself. Primed wheat seeds exhibit higher germination and stress tolerance than non-primed seeds under adverse conditions. Additionally, early-season stress exposure in wheat can improve antioxidant activity and grain yield under subsequent stresses. Collectively, these studies illustrate that plants can undergo physiological and biochemical adjustments during an initial stress that persist and provide an advantage when stress recurs. While many of these responses are transient, mounting evidence indicates that drought–heat signaling (via ABA, ROS, and Ca²⁺ waves) can recruit transcription factors and chromatin remodelers, depositing epigenetic marks (e.g., DNA methylation changes, H3K4me3/H3K27me3, and small RNAs) that encode stress experience into longer-term and occasionally intergenerational memory [5]. This mechanistic link motivates our tests of whether acute physiological shifts align with heritable phenotypic responses across generations [6,7]. The mechanisms underlying stress recall are complex and multi-level, involving changes in gene expression, accumulation of protective metabolites, and adjustments in hormonal and redox signaling [4]. Notably, many of these changes are associated with epigenetic modifications, heritable alterations in gene regulation (such as DNA methylation and chromatin modifications) that do not involve changes in the DNA sequence. These epigenetic marks can serve as a molecular memory of stress, enabling a primed state that can be reactivated upon future stress exposure.

Crucially, plant stress memory is not limited to within a single generation. There is growing evidence that stress memory can be transmitted to progeny, a phenomenon known as transgenerational stress memory. While most stress-induced epigenetic modifications are reset during reproduction, a subset can escape this “reprogramming” and be inherited by the next generation [8]. Intergenerational memory (effects on the immediate offspring) and true transgenerational event (persistence into grand-offspring or beyond) have been documented in plant systems, indicating that the offspring of stressed plants may exhibit enhanced stress tolerance even if they themselves were not directly exposed. This transgenerational inheritance is thought to be mediated by stable epigenetic changes, for instance, DNA cytosine methylation patterns or histone marks that remain through gametogenesis and embryogenesis. Small interfering RNAs and other chromatin-associated factors can also carry stress “information” into the next generation [9,10]. In Arabidopsis, recurrent heat can establish “thermomemory” via HSFA2 centered regulatory loops and chromatin remodeling (e.g., H3K27me3 dynamics), with some heat-induced states detectable for 1–2 generations, although RdDM/siRNA pathways and remodelers like DDM1/MOM1 generally restrict transmission [11]. In rice, multigenerational selection under drought or salinity reshapes DNA-methylation landscapes, producing stable epimutations linked to improved offspring performance [12]. Understanding and harnessing this phenomenon could open new avenues for crop improvement under climate stress conditions [8].

Genomics now enables the dissection of the genetic and epigenetic bases of stress persistent acclimation in wheat. High-throughput mapping (GWAS/QTL) is revealing loci for multi-generational resilience. While GWAS cannot directly isolate epigenetic inheritance, it can detect loci whose effects are consistent across stress-experienced generations, suggesting a heritable component. We therefore interpret these loci as candidates for heritable stress memory, not definitive proof of epigenetic causality, which will require functional and molecular validation. In wheat, [13] contrasted accessions with vs. without ancestral drought, identifying loci affecting spike length, grain number/weight, and antioxidant levels, including a chromosome 2ALRR candidate associated with superior performance across three generations, implicating stress-signaling–mediated priming [13]. In barley, [14] detected 332 SNPs in 14 regions (2H, 3H, 4H, 5H, 7H) linked to improved growth and biochemical traits in progeny of stressed plants, underscoring polygenic control and pathway-level integration. Complementarily, wheat epigenomic profiling (whole-genome DNA methylation; ChIP-seq) is mapping stress-imprinted marks [15], while CRISPR/dCas9 epigenome editing enables causal tests via targeted methylation/histone modulation [8]. GWAS studies have identified QTLs linked to physiological and agronomic traits under combined or heat–drought target environments, including chlorophyll content, canopy temperature, stomatal conductance, grain-filling duration, and yield stability. Multi-environment GWAS in wheat have repeatedly mapped loci for chlorophyll/SPAD, canopy temperature, and yield components under drought/heat panels, underscoring the value of these traits as selection proxies in stress-prone environments [16]. Associations between canopy temperature and stomatal conductance further anchor thermal phenotypes to gas-exchange biology, with genetic regions for cooler canopies co-varying with higher gs and yield under water limitation [17]. Under explicitly combined heat–drought conditions, additional MTAs for yield and related traits have been reported, providing targets for marker-assisted improvement in harsh mega-environments [18]. Grain-filling dynamics are also genetically tractable: QTL/GWAS efforts detect loci affecting post-anthesis heat responses and grain-filling duration that co-influence final yield stability [19]. Together, integrated mapping, epigenomics, and multi-omics across generations are clarifying how a heritable “ recall” of drought/heat is encoded and maintained. Integrating the concept of transgenerational stress signature into wheat improvement holds great promise for developing climate-resilient varieties. By identifying genetic markers and epigenetic signatures associated with inherited stress tolerance, breeders can incorporate these into marker-assisted selection or genomic selection frameworks. For instance, loci discovered through GWAS of stress imprinting (such as the wheat 2ALRR gene or the barley QTLs) provide targets for breeding or biotechnological intervention to enhance multi-generational resilience. Moreover, understanding stress memory mechanisms can inform novel breeding strategies like “stress priming” of breeding populations or parental lines. This could involve exposing plants to controlled stress at certain developmental stages to induce beneficial persistent acclimation, which is then passed to seeds used for cultivation, effectively pre-conditioning future crops for anticipated drought or heat episodes. Therefore, the emerging field of transgenerational stress memory, bolstered by cutting-edge genetic mapping and molecular techniques, provides a new paradigm for crop improvement. Harnessing heritable stress.

imprinting in wheat could accelerate the development of varieties that maintain high yields under recurrent drought and heat stress, thus contributing to sustainable food security in the face of climate change [8,13]. This work aims to dissect and ultimately exploit the heritable architecture of wheat resilience to combined drought and heat (DH) by mapping loci and molecular features underpinning transgenerational stress signature. Specifically, we seek to (i) quantify multi-generational phenotypic responses under control versus stress regimes, (ii) derive generation-aware stress tolerance indices that capture heritable effects, and (iii) perform genome-wide association/QTL analyses with longitudinal and multi-trait models to identify loci whose effects persist or intensify across generations. We further aim to integrate genomic and transcriptomic signatures with mapping results to prioritize causal candidates, test the hypothesis that inherited priming contributes measurably to stress tolerance, and distinguish memory-specific loci from single-stress effects. The ultimate objective is to validate top candidates through targeted expression/priming assays and deliver breeder-ready markers and genomic prediction models that accelerate the development of climate-resilient wheat cultivars capable of maintaining yield under recurrent combined drought and heat episodes.

Plant materials

This study utilized a collection of 111 spring wheat accessions sourced from the IPK-Gatersleben Genebank. The collection comprised 57 accessions from Europe, 31 from Asia, 10 each from North America and South America, one from Australia, and two of unknown geographical origin. Most of these wheat accessions were originally released in the 1970s. The genetic diversity within the allohexaploid and allotetraploid wheat populations was assessed using a genotyping array consisting of 90,000 gene-associated single-nucleotide polymorphisms (SNPs). Further methodological details are available in Muqaddasi et al. (2016).

Field experiments

The present study was conducted on a diverse panel of bread wheat (Triticum aestivum L.) genotypes (Table S1). To ensure genetic integrity, seeds were propagated through controlled selfing before initiation of the experiments at Fayoum experimental station. The Fayoum site has a hot, arid desert climate with rainfall confined to winter and near-zero precipitation in summer. Temperatures are mild in winter and rise sharply by late spring–summer, with daytime maxima frequently ~ 33–38 °C and warm nights. Relative humidity is generally low–moderate, trending ~ 50–60% during hot spells, creating high evaporative demand that intensifies drought (Fig. S1 and Table S2). All generations (T1–T3) were cultivated in pots placed in the open field at the experimental station. Pots (dimensions and soil mix as specified above) were set on a leveled gravel bed to ensure drainage and prevent root escape, spaced ≥ 40 cm within rows and ≥ 60 cm between rows, with buffer rows around the experiment to minimize edge effects. After emergence, plants were thinned to one plant per pot to avoid intra-pot competition. Experiments were arranged in randomized blocks with three biological replicates per genotype and condition to ensure statistical robustness. Plants were evaluated across three successive cycles of exposure (T1–T3), with each generation subjected to either control conditions or combined drought and heat stress. In the control regime, plants were irrigated to maintain soil at field capacity and grown under ambient temperatures of 22–25 °C throughout the growth cycle. In contrast, the dual stressed regime imposed concurrent water limitation, maintaining soil moisture at ~ 35% of field capacity, and heat exposure by sowing on 1 January each year so that the late-vegetative through post-anthesis and grain-filling stages coincided with daytime maxima of ~ 33 °C. Each treatment–generation combination was designated as T1C1, T1S1, T2C2, T2S2, T3C3, and T3S3. Stress imposition was carefully controlled. Soil water content was monitored gravimetrically, and irrigation was withheld until the desired drought threshold was reached and maintained for 14 days.

Agronomic traits

Phenotypic and agronomic traits were recorded at maturity using standardized descriptors. Plant height was measured from the soil surface to the tip of the spike (excluding awns), while the number of grains per spike were recorded manually from harvested spikes, and thousand-kernel weight was calculated by weighing a subsample of 1000 fully developed grains.

Physiological measurements

Physiological and biochemical traits were assessed using fresh flag leaves collected at anthesis. Proline content was determined using the acid ninhydrin method of [20] and expressed as µmol g⁻¹ fresh weight. Chlorophyll was quantified following [21], with pigments extracted in 80% acetone and absorbance recorded at 663 and 645 nm. Soluble protein content was measured according to [22] using the Coomassie dye-binding method.

Antioxidant assays

Enzymatic antioxidant activities were quantified in crude extracts prepared in phosphate buffer containing polyvinylpyrrolidone to minimize phenolic interference. Superoxide dismutase (SOD) activity was measured following [23] by its ability to inhibit nitroblue tetrazolium reduction. Catalase (CAT) activity was assayed as the rate of H₂O₂ decomposition at 240 nm according to [24]. Ascorbate peroxidase (APX) activity was determined by monitoring the oxidation of ascorbate at 290 nm [25], while glutathione reductase (GR) activity was assayed by following NADPH oxidation at 340 nm [26]. Non-enzymatic antioxidants were also profiled: ascorbate (AsA) was extracted in an acid medium and quantified spectrophotometrically at 265 nm, and glutathione (GSH) was measured through its reaction with DTNB (Ellman’s reagent) at 412 nm.

Genome-wide association study (GWAS)

SNP genotyping with the 90 K wheat array

High-density SNP genotyping was performed using the Illumina Infinium iSelect 90 K wheat SNP array [27]. This array interrogates approximately 90,000 single-nucleotide polymorphisms that are predominantly gene-associated and well-distributed across the A, B, and D subgenomes of polyploid wheat. Genotyping assays were carried out according to the manufacturer’s instructions. In brief, genomic DNA from each sample was amplified, fragmented, and hybridized to the 90 K SNP beadchip, after which allele-specific primer extension and staining were performed. The processed arrays were scanned on an Illumina iScan system, and raw fluorescence intensity data for each marker were collected. Genotype calling was conducted using Illumina’s GenomeStudio software (version 2011.1). To ensure accurate allele clustering in the polyploid wheat genome, we applied the calling procedures described by [27]. Specifically, the GenomeStudio genotyping module (diploid setting) was used to cluster signal data for each SNP, as recommended for hexaploid wheat. Any ambiguous genotype clusters were manually inspected, and residual heterozygous calls were conservatively treated as missing data (since the lines are largely homozygous inbred cultivars). This approach follows the best practices established for the wheat 90 K array to minimize false heterozygotes and genotype calling errors.

Genotype data processing and quality control

Following initial genotype calling, missing data imputation and marker quality filtering were performed to prepare the dataset for analysis. First, any missing genotype calls were imputed using FastPHASE version 1.3 (default parameters) to infer likely allele states at uncalled loci. This step helped to reduce the proportion of missing data and increase the statistical power for downstream analyses. Next, stringent criteria were applied to filter out low-quality or uninformative SNP markers. Quality filtering criteria for SNP markers included markers with more than 10% missing genotype calls were removed from the dataset. Markers with a minor allele frequency (MAF) below 5% (MAF < 0.05) were excluded from further analysis. These filtering thresholds ensured that only high-confidence, polymorphic SNPs were retained. After applying the filters, the remaining SNP markers provided robust genome-wide coverage for the diversity analysis. All subsequent analyses (e.g., population structure, association mapping) were therefore conducted on this curated set of SNP genotypes, as described by [27] and others using the 90 K wheat array.

Baseline population-genomic statistics

To document data quality and panel diversity, we reported Transition/transversion (Ts/Tv) ratio computed from biallelic SNPs as a global QC indicator [28]. Polymorphic Information Content (PIC) per SNP following [29]. Expected heterozygosity (He) and Nei’s gene diversity across loci and per chromosome/subgenome [30]. Genome-wide and subgenome-specific LD decay: pairwise r² within sliding 10-Mb windows, LOESS fit of r² vs. distance, and decay distance defined at r² = 0.2 [31,32]. These metrics contextualize mapping resolution and SNP-to-gene distances in subsequent locus annotation.

Population structure and relatedness

We inferred structure and kinship from the 90 K SNP dataset after QC (missing ≤ 10%, MAF ≥ 0.05) and LD-pruning (PLINK). Population structure was summarized by PCA on the pruned set (top 3 PCs retained as covariates) and by ADMIXTURE with K = 1–10, selecting K by minimum cross-validation error; ancestry proportions (Q) were extracted for downstream use [33,34,35]. Pairwise relatedness was captured with a marker-based kinship matrix (VanRaden method) computed from the unpruned QC set [36,37]. For GWAS, we fitted mixed models with PCs (and/or Q) to control structure and K to model relatedness (e.g., FarmCPU with K/PC covariates), following standard wheat-array practice [38].

GWAS models and cross-model validation

Given the modest panel size and wheat’s LD, we prioritized models that balance power and confounding control. We ran two complementary, state-of-the-art multi-locus approaches implemented in GAPIT: FarmCPU and BLINK [39,40][38]. FarmCPU iteratively fits associated markers as fixed covariates while optimizing random effects, thereby reducing confounding and improving power [39]. BLINK replaces the random effect with Bayesian-informed marker selection using LD information to speed computation and reduce overfitting [40]. To strengthen robustness and address concerns about model-specific signals, we report as high-confidence, those SNP-trait associations significant in both FarmCPU and BLINK with concordant effect directions and similar peak intervals. Associations detected by a single model are labeled as model-specific and interpreted cautiously, following recommended practice for reproducibility [39,40]. Genome-wide significance used Benjamini–Hochberg FDR q < 0.05 on the full marker set for each trait; suggestive signals at q < 0.10 are reported as exploratory [41]. To provide distribution-free support, we conducted 1,000 phenotype-label permutations per trait in FarmCPU to derive empirical p-value thresholds, which closely matched FDR cutoffs. Genomic inflation (λ) was monitored from QQ-plots; λ values close to 1 indicated adequate confounding control [42]. Manhattan/QQ plots are provided for every trait/model. Significant SNPs within LD-based windows were clumped into quantitative trait loci (QTL) using physical distance informed by the panel’s LD decay per chromosome/subgenome (e.g., ± LD₀.₂ Mb around the lead SNP). When multiple peaks occurred within a window but were in low LD (r² < 0.2), they were retained as distinct loci [32,43].

Candidate-gene annotation

Lead SNPs and their LD blocks were intersected with high-confidence gene models (± LD_0.2) to list nearby candidates. Functional annotation used Ensembl Plants/RefSeq gene descriptionshttps://plants.ensembl.org, GO terms, and InterPro domains. Where available, we summarized stress-relevant expression from public wheat expression browsers to help prioritize candidates [44].

Sensitivity analyses

To evaluate robustness, we re-ran both models under alternative covariate sets (0, 3, 5 PCs; with/without K), and under stricter QC (MAF ≥ 0.07; SNP missingness ≤ 0.05). Concordance of lead SNPs (Jaccard index) and stability of effect sizes across runs were summarized. We also compared results using BLUEs vs BLUPs for traits with heterogeneous error variance [45].

Software and parameters

All analyses were performed in R (≥ 4.3) using GAPIT [38], GAPIT3 updates, and the FarmCPU/BLINK implementationshttps://zzlab.net/GAPIT [39,40]. Additional utilities included PLINK 1.9/2.0 for pruning/LDhttps://www.cog-genomics.org/plink [43,46], ADMIXTURE v1.3 [35],https://github.com/DReichLab/EIG and Beaglehttp://dalexander.github.io/admixture [47].

Plant material and stress treatment

Two wheat genotypes contrasting for stress response were used: a tolerant line and a susceptible line. Plants were grown in pots containing a peat:perlite mix (2:1, v/v) in a controlled environment ([photoperiod, e.g., 16 h light/8 h dark], [light intensity], [RH]) at 22/18 °C day/night until the [tissue stage; e.g., fully expanded third leaf/flag leaf, Zadoks 37–39].

Combined drought and heat was imposed by withholding irrigation to 35% pot field capacity (verified gravimetrically) while raising air temperature to 36/26 °C day/night. Leaf relative water content (RWC) was monitored on subsamples to verify drought progression. Control plants were kept well-watered at 22/18 °C.

Sampling design

Leaf tissue was sampled at 0, 12, 24, and 36 h from the onset of stress (0 h = pre-stress baseline for each genotype). For each genotype × time, n = 3 biological replicates (independent plants) were harvested between 10:00–12:00 to minimize diurnal effects. Tissue was flash-frozen in liquid N₂ and stored at − 80 °C.

RNA extraction and quality control

Total RNA was isolated from ~ 80 to 100 mg frozen tissue using [kit/protocol; e.g., Spectrum Plant Total RNA Kit (Sigma) or TRIzol] following the manufacturer’s instructions, including on-column DNase I treatment to remove genomic DNA. RNA quantity and purity were assessed by spectrophotometry (A260/280 = 1.9–2.1; A260/230 > 2.0). Integrity was checked by agarose gel electrophoresis and/or Bioanalyzer; samples with RIN ≥ 7.0 were used.

cDNA synthesis

One microgram of RNA per sample was reverse-transcribed using [reverse transcriptase] with a mixed oligo(dT) + random hexamer primer strategy in a 20 µL reaction ([cycle: 25 °C 10 min; 50 °C 30–45 min; 85 °C 5 min]). cDNA was diluted 1:5 with nuclease-free water for qPCR.

Candidate genes and primer design

Candidate genes were: Arginine/serine-rich splicing factor, Trehalose-6-phosphate phosphatase, Vacuolar iron transporter, Pectinesterase, Methyltransferase type-12 domain-containing protein, and Glycosyltransferase. Primers (90–180 bp amplicons) were designed in Primer3 with Tm = 60 ± 1 °C, GC 40–60%, and checked for hairpins/dimers. Where possible, primers spanned exon–exon junctions. Specificity was confirmed by BLAST against the wheat genome [44] and by single-peak melt curves and single bands on 2% agarose gels. Primer efficiency (E) was determined from a 5-point fivefold dilution series: E = 10^(− 1/slope) and accepted within 90–110% (R² ≥ 0.99). Reference gene candidates (e.g., ACTIN, GAPDH, EF-1α) were evaluated across all genotypes and timepoints. Stability was assessed with geNorm / NormFinder; the geometric mean of the top 2–3 reference genes was used as the normalization factor [48]. Relative expression for each target was computed as 2⁻^ΔΔCq using 0 h within each treatment as the calibrator [49]. MIQE guidelines were followed for assay design and reporting [50]. All the primers used for the candidate genes are described in Table S3.

qPCR conditions

qPCR was performed on a [instrument model] using SYBR Green chemistry ([master mix]) in 10 µL reactions: 1 × master mix, 0.3 µM each primer, and 2 µL diluted cDNA. Cycling: 95 °C 2 min, then 40 cycles of 95 °C 15 s, 60 °C 30 s, followed by a melt curve (65–95 °C, 0.5 °C increments). Each biological replicate had two technical replicates; NTCs and no-RT controls were included. Technical replicates differing by > 0.5 Ct were repeated or excluded.

Statistical analysis

Analyses were conducted in Python (scipy, statsmodels) and plotted with matplotlib (300 dpi). For each gene and genotype, Welch’s t-test compared 12, 24, 36 h to the corresponding 0 h baseline (two-tailed). Where multiple time point tests were run, Benjamini–Hochberg FDR correction was applied within each gene × genotype family (q < 0.05 deemed significant). An orthogonal analysis was performed using a two-way model (genotype, time, genotype × time) on ΔCt values, which was fitted when normality/homoscedasticity assumptions were met. Tukey-HSD post-hoc tests were used to examine simple effects. Effect sizes were summarized as fold-change (2^−ΔΔCt) and mean differences in ΔCt (± SE) in Table S4. Barplots show mean 2 − ^ΔΔCt ± SEM (n = 3 biological replicates) for tolerant and susceptible genotypes at 0, 12, 24, 36 h. Asterisks indicate significance versus 0 h within genotype (* p < 0.05;** p < 0.01;*** p < 0.001). A two-factor analysis of variance (ANOVA) was conducted to examine the effects of genotype (different wheat accessions) and treatment (control versus combined cold and water deficit stress) on essential phenotypic traits, including the investigation of possible interactions between these variables. This statistical method allowed for the detailed evaluation of genotype-specific responses under both stressed and optimal conditions [51].

Multigenerational responses to combined drought and heat stress

Exposure to combined drought and heat stress across three successive cycles (T1–T3) induced marked and consistent alterations in both biochemical defense systems and yield-related agronomic traits (Fig. 1). Boxplot comparisons demonstrated highly significant control–vs–stress contrasts in most traits, with patterns of sustained induction or partial attenuation depending on the physiological process. Enzymatic antioxidants, including ascorbate peroxidase (APX), catalase (CAT), superoxide dismutase (SOD), and glutathione reductase (GR), were consistently upregulated under dual stress relative to control conditions across all three generations. The magnitude of separation was greatest in T1 and T2, with elevated activities persisting into T3 (Fig. 1). Similarly, non-enzymatic antioxidants, ascorbate (AsA) and glutathione (GSH) displayed higher pool sizes under stress in every generation. These results indicate a robust and durable reinforcement of ROS-scavenging capacity, consistent with stress-induced redox priming mechanisms reported in cereals. Proline (ProC) accumulated significantly under combined drought and heat stress across T1–T3, mirroring antioxidant dynamics and underscoring its dual roles in osmoprotection and redox buffering. Soluble protein (SP) levels were likewise enhanced by stress, with particularly strong differences by T3 (Fig. 1). This persistent increase suggests proteome remodeling toward stress-responsive proteins and metabolic reallocation, favoring defense over growth. By contrast, chlorophyll content (Chl) was consistently reduced under stressed conditions relative to control plants, reflecting limitations in photosynthetic machinery under combined osmotic and thermal stress. Interestingly, the stress-induced reduction was somewhat less pronounced in T3 compared with earlier cycles, suggesting partial acclimation or retention of pigment stability with repeated exposure.

Boxplots illustrate the phenotypic differences observed across three generations for all evaluated physio-morphological traits. The initial generation (T1) experienced combined drought and heat stress during the 2020/2021 season and was compared to both control (T1C1) and stress-exposed (T1S1) groups. Seeds from T1 were used to cultivate the second generation (T2), which was subjected to either the same combined drought and heat stress (T2S2) or control conditions (T2C2), enabling assessment of intergenerational stress memory effects during the 2021/2022 season. Furthermore, a third generation (T3) was developed to investigate both transgenerational and intergenerational stress memory during the 2022/2023 season, under either control (T3C3) or combined drought and heat stress conditions (T3S3)

Full size image

For morphological and yield traits, dual stress imposed significant penalties on structural and reproductive traits. Plant height (PH), spike length (SL), number of spikelets per spike (NSS), grain number per spike (NGS), and thousand-kernel weight (TKW) were all significantly decreased under stress in every generation (Fig. 1). However, reductions in some yield-related parameters (e.g., PH and TKW) appeared less severe in T3, indicating multigenerational acclimatory adjustments that mitigate, but do not eliminate, yield penalties. These results align with observations that repeated stress exposure can induce a degree of transgenerational resilience, albeit at a cost to absolute productivity.

Stress-tolerant index effects

In general, biochemical/antioxidant traits (CAT, SOD, APX, GR, GSH, AsA), osmolyte and protein pools (ProC, SP) increased under stress, whereas chlorophyll (Chl) and yield/architecture metrics (PH, SL, NGS, NSS, TKW) declined (Fig. 2). Under Stress, CAT, SOD, APX, and GR display clear median up-regulation in all generations, with the separation between Stress and Control widening by T3 for most enzymes (e.g., APX, SOD, CAT). The metabolite pools (GSH, AsA) follow the same direction, and ProC (proline) shows one of the largest and most consistent inductions across T1–T3, reflecting osmotic adjustment. SP (soluble protein) also rises, with a pronounced elevation in T3, consistent with stress-induced synthesis and/or accumulation of protective proteins. These patterns indicate progressive reinforcement of the redox buffer and osmoprotection machinery across generations (Fig. 2). Together, these results support a priming-like trend in biochemical defense. For the photosynthetic pigment Chl, the Level consistently decreases under Stress in every generation. The Chl_STI also trends upward from STI1 to STI3 (with strong significance for several pairwise contrasts), consistent with partial maintenance of photosynthetic capacity in later generations as acclimation develops. All yield-related traits show sharp depressions under stress in T1 (Fig. 2). By T3, the gap is reduced for PH and TKW in particular, indicating partial recovery of stature and grain mass under repeated exposure. The STI trajectories capture these differences: PH_STI improves by T3 (significant for STI2 → STI3 and STI1 → STI3), TKW_STI increases strongly toward STI3 (highly significant contrasts), while SL_STI is largely stable with evidence of an early step change. In contrast, NGS_STI and NSS_STI are flat across generations (non-significant pairwise tests), indicating that grain set/spikelet number remains a bottleneck under combined stress even as other components acclimate (Fig. 2). This divergence, mass and height recovering more than grain number, suggests that reproductive processes are more refractory to transgenerational acclimation than vegetative growth or grain filling.

Boxplots demonstrate stress tolerance indices (STI) across three generations for all evaluated physio-morphological traits. Statistical significance levels are denoted as follows:^nsp, 0.1;*p, 0.05;**p, 0.01;****p, 0.0001

Full size image

Correlation analysis

This integrated Pearson heat-map/Mantel network shows how agronomic and physiological indices co-vary across the three generations (STI1–STI3) as described in Fig. 3. A dense block of strong, highly significant positive correlations is evident among yield-related traits (e.g., PH, NSS, TKW) across generations, indicating stable, coordinated variation. Physiological indices also cluster within their class (e.g., chlorophyll, soluble protein, and antioxidant enzymes), but cross-class links to yield traits are fewer and generally weaker. Notably significant positive associations for chlorophyll and soluble protein with kernel- and spike-related indices, whereas some antioxidant measures exhibit only modest or occasional negative ties. The Mantel network corroborates this pattern: the thick, mostly solid edges concentrate on spike/yield and soluble-protein traits, implying they dominate the overall sample structure, with enzyme activities contributing smaller, trait-specific effects.

Pearson’s correlation analysis was conducted among all stress tolerance indices (STI) for the evaluated physio-morphological traits across the three generations

Full size image

Genome-wide marker–trait associations

We identified 97 markers–trait associations (MTAs) spanning 24 stress-tolerance indices and 56 unique SNPs across 16 chromosomes, with enrichment on 1B and 7B and secondary peaks on 2A/2B/6A using FarmCPU model (Table S5 and Fig. S2 and 3). The strongest signal overall was for plant height under stress (PH_STI2) at wsnp_Ex_rep_c66465_64708628 (2A:686,860,574), where the T allele exerted a large negative effect on stature (β≈ − 8.35; p reported as 0), consistent with stress-adaptive dwarfing. A distal 1B haplotype around ~ 581.20–581.21 Mb emerged as a pleiotropic hotspot: Tdurum_contig11896_550/Kukri_c9752_793/Tdurum_contig69792_167 collectively increased chlorophyll indices (Chl_STI1/2; β≈ + 1.2) and thousand-kernel weight (TKW_STI2; β≈ + 2.05) but decreased soluble-protein under stress (SP_STI3; β≈ − 5.18), indicating coordinated yet trait-dependent effects within the same LD block.

Additional multi-trait loci reinforced this architecture. On 3A (725.74 Mb), wsnp_BE426418A_Ta_2_1 is associated with four traits, elevating SP_STI2/STI3 and PH_STI2 while reducing Chl_STI1, suggesting linked regulation of protein allocation, plant stature, and photochemistry. Antioxidant capacity partitioned cleanly: RAC875_c35438_474 on 2B increased both GR_STI2 and GSH_STI3 (large positive β), whereas Tdurum_contig51640_847 on 6B consistently elevated ascorbate indices (AsA_STI2/3). Yield components were further supported by group-6/7 signals, including IAAV5595 (6B; TKW_STI2 and NSS_STI3) and Kukri_c38390_218 (7A; NGS_STI1/STI3). Overall, a compact set of reproducible intervals—most prominently 2A:686.86 Mb and the distal 1B cluster jointly modulate photochemical stability, antioxidant defense (2B glutathione; 6B ascorbate), and yield formation, with consistent allelic directions within pathways and predictable trade-offs across traits.

Genomic regions associated with stress-tolerance indices

We delineated seven linkage-disequilibrium (LD)–defined genomic regions harboring 19 markers–trait associations (MTAs) across nine stress-tolerance indices (Chl_STI1/2, SP_STI3, AsA_STI2/3, GSH_STI3, GR_STI2, PH_STI2, TKW_STI2) as shown in Table S6. Signals were strongly enriched on chromosome 1B (10/19 MTAs), followed by 2A (3/19), indicating group-1 and group-2 hotspots in this panel (Table S6). The most significant association overall was for plant height under stress (PH_STI2) on chromosome 2A, while a distal 1B cluster (~ 581.20–581.21 Mb) showed reproducible, pleiotropic effects spanning chlorophyll, soluble protein, and kernel weight.

A major 2A locus (Region 5; 2A:686,860,574; wsnp_Ex_rep_c66465_64708628) exerted a dominant effect on stature under stress: the T allele markedly reduced PH_STI2 (β =  − 8.35,p = 3.96 × 10⁻¹⁴). Notably, the same allele improved Chl_STI1 (β = 1.45,p = 2.16 × 10⁻⁶) and TKW_STI2 (β = 2.63,p = 9.63 × 10⁻⁶), suggesting a stress-adaptive architecture coupling reduced height with enhanced photochemistry and kernel mass (Fig. 4). The distal 1B cluster encompassed three tightly linked tags that together mapped a robust pleiotropic interval. At Region 2 (1B:581,201,755; Tdurum_contig11896_550), the C allele increased Chl_STI1 (β = 1.23,p = 9.96 × 10⁻⁸), Chl_STI2 (β = 1.21,p = 8.66 × 10⁻⁶), and TKW_STI2 (β = 2.05,p = 1.19 × 10⁻⁶), but decreased SP_STI3 (β =  − 5.18,p = 6.29 × 10⁻⁶). Concordant profiles were observed at Regions 3 and 4, tagged by Kukri_c9752_793 (1B:581,205,619; G allele: Chl_STI2 β = 1.21,p = 8.66 × 10⁻⁶; TKW_STI2 β = 2.05,p = 1.19 × 10⁻⁶; SP_STI3 β =  − 5.18,p = 6.29 × 10⁻⁶) and Tdurum_contig69792_167 (1B:581,205,619; C allele with identical effect sizes and p values) as shown in Table S6. Together, these tags indicate a shared haplotype spanning ~ 581.20–581.21 Mb that promotes chlorophyll maintenance and kernel weight while modulating soluble-protein accumulation in the opposite direction under stress.

a The locus zoom plot highlights marker wsnp_Ex_rep_c66465_64708628 on chromosome 2A, which is associated with essential traits related to stress memory. The x-axis shows SNP order along chromosomes, while the y-axis indicates the − log10 (p value) for each SNP marker.b Locus zoom of highly significant SNP located on chromosome 2A.cThe diagram illustrates the structure of the candidate geneTraesCS2A02G434000, annotated as a vacuolar iron transporter. With the linkage disequilibrium (LD) pattern.d The quantile–quantile (Q-Q) plots compare the observed SNP association values with the expected distribution for the evaluated traits

Full size image

Two additional regions aligned with antioxidant pathways. Region 6 (2B:106,939,395; RAC875_c35438_474) supported glutathione-linked responses: the A allele increased GR_STI2 (β = 47.26,p = 6.90 × 10⁻⁶) and GSH_STI3 (β = 64.80,p = 8.10 × 10⁻⁶), indicating coherent, positively aligned effects on glutathione metabolism (Fig. 5). Region 7 (6B:674,103,633; Tdurum_contig51640_847) consistently underpinned ascorbate capacity, with the G allele elevating AsA_STI3 (β = 40.66,p = 6.16 × 10⁻⁶) and AsA_STI2 (β = 34.75,p = 1.36 × 10⁻⁵), nominating this 6B interval as a principal determinant of low-molecular antioxidant pools in stress conditions. Finally, Region 1 on 1A (1A:31,779,263; BobWhite_c4499_153) illustrated pleiotropy with opposing directions across traits: the T allele increased Chl_STI1 (β = 1.16,p = 1.88 × 10⁻⁸) but decreased PH_STI2 (β =  − 3.34,p = 5.03 × 10⁻⁶) as shown in Fig. 6. This trade-off enhanced photochemical maintenance alongside reduced stature mirrors the pattern at the 2A locus and reinforces a recurring theme of stress-adaptive reallocation. In sum, a compact set of LD-defined intervals, most prominently 2A:686.86 Mb and the distal 1B cluster around ~ 581.20–581.21 Mb coordinate photochemical stability, antioxidant defense (partitioned between 2B for glutathione and 6B for ascorbate), and yield components under combined stress. The consistency of allelic effects within pathways (2B, 6B) and the trait-dependent directions at 1B and 2A together outline a genetically tractable framework for improving stress resilience while managing trade-offs among physiology and yield.

a The locus zoom plot highlights marker RAC875_c35438_474 on chromosome 2B, which is associated with essential traits related to stress memory. The x-axis shows SNP order along chromosomes, while the y-axis indicates the − log10 (p value) for each SNP marker.b Locus zoom of highly significant SNP located on chromosome 2B.c The diagram illustrates the structure of the candidate geneTraesCS2B02G140200, annotated as a methyltransferase type 12 domain-containing protein with the linkage disequilibrium (LD) pattern.d The quantile–quantile (Q-Q) plots compare the observed SNP association values with the expected distribution for the evaluated traits

Full size image

a The locus zoom plot highlights marker BobWhite_c4499_153 on chromosome 1A, which is associated with essential traits related to stress memory. The x-axis shows SNP order along chromosomes, while the y-axis indicates the − log10 (p value) for each SNP marker.b Locus zoom of highly significant SNP located on chromosome 1A.c The diagram illustrates the structure of the candidate geneTraesCS1A02G049700, annotated as an arginine/serine-rich splicing factor with the linkage disequilibrium (LD) pattern.d The quantile–quantile (Q-Q) plots compare the observed SNP association values with the expected distribution for the evaluated traits

Full size image

Candidate genes and Gene expression using RT-qPCR

Candidate genes near the recurrent peaks reinforce the mechanism: 1B (≈581.20–581.21 Mb) includes Trehalose-6-phosphate phosphatase (TraesCS1B02G350600), Rab7 (TraesCS1B02G350200), and γ-glutamyl hydrolase, linking sugar signaling, vesicle trafficking, and folate metabolism to Chl/SP/TKW (Fig. 7). Chromosome 2A (686.86 Mb) features vacuolar iron transporters (TraesCS2A02G434100/200), XTH (TraesCS2A02G434300), a VQ regulator, actin-bundling protein, and plastid-lipid-associated protein, consistent with growth, wall remodeling, Fe/chloroplast homeostasis (PH/Chl/TKW). 2B (106.94 Mb) carries pectinesterase, methyltransferase-type 12, and a CCT-domain factor, matching wall/redox and transcriptional timing (GR/GSH). 6B (674.10 Mb) harbors F-box, Ser/Thr phosphatase, glycosyltransferase, and shugoshin candidates, pointing to proteostasis/post-translational control of ascorbate responses. Overall, a compact set of pleiotropic loci (1B/2A/2B/6B) underpins both stress physiology and yield variation. Across all six candidates, the tolerant genotype mounted a stronger and earlier transcriptional response than the susceptible genotype, with the largest expression levels consistently observed at 36 h. Baseline (0 h) expression was similar between genotypes, indicating that stress exposure, rather than constitutive differences, drives the observed divergence. Trehalose 6-phosphate phosphatase and pectinesterase showed significant induction by 12 h in the tolerant genotype, pointing to rapid mobilization of osmoprotective carbohydrate metabolism and cell-wall remodeling as first-line responses. The arginine/serine-rich splicing factor also began rising by 12 h, suggesting early adjustments in pre-mRNA processing to reprogram stress-responsive isoforms. By 24 h, induction broadened: vacuolar iron transporter, methyltransferase type 12 domain protein, and glycosyltransferase were significantly elevated in the tolerant genotype. These patterns are consistent with the onset of redox/metal homeostasis buffering (via vacuolar iron sequestration), transcriptional/epigenetic fine-tuning (methyltransferase-linked processes), and metabolite glycosylation to stabilize reactive intermediates. The susceptible genotype generally displayed smaller effect sizes and fewer significant changes at this stage. All genes peaked or remained high at 36 h in the tolerant genotype, with glycosyltransferase and trehalose 6-phosphate phosphatase among the strongest responders. This late amplification is indicative of consolidated metabolic reprogramming, osmolyte maintenance, structural reinforcement, detoxification, and sustained regulatory control, supporting stress endurance under prolonged combined drought and heat. In contrast, the susceptible genotype exhibited only moderate late increases, implying limited ability to sustain these protective pathways. The consistent separation between blue (tolerant) and red (susceptible) bars at 24–36 h across multiple genes highlights a coordinated, polygenic transcriptional program in the tolerant background. Functionally, this program spans RNA processing (splicing factor), carbohydrate/trehalose metabolism (T6PP), metal/iron homeostasis (VIT), cell-wall dynamics (pectinesterase), epigenetic regulation (methyltransferase-type 12), and glycosylation/secondary metabolism (glycosyltransferase), a mechanistic ensemble plausibly underlying superior physiological resilience under combined stress.

The relative gene expression for the genes encodinga Arginine/serine-rich splicing factor,b Trehalose 6-phosphate phosphatase,c Vacuolar iron transporter,d Pectinesterase,e Methyltransferase type 12 domain-containing protein, and (f) Glycosyltransferase under combined drought + heat. Bars show mean 2-ΔΔCt ± SEM (n = 3) for tolerant and susceptible wheat genotypes at 0, 12, 24, and 36 h. Asterisks denote significance versus 0 h within the same genotype (Welch’s t-test:* p < 0.05,** p < 0.01,*** p < 0.001)

Full size image

Transgenerational responses to combined drought and heat stress in wheat

Across three successive cycles (T1–T3), combined drought and heat elicited a durable reinforcement of wheat’s redox-defence and osmotic machinery. Antioxidant enzymes (SOD, CAT, APX, GR) were consistently higher under combined drought and heat stress than controls in all generations, indicating sustained ROS-scavenging capacity [51,52]. Non-enzymatic pools (AsA, GSH) likewise increased, consistent with redox priming whereby prior stress “trains” stronger subsequent responses [53,54]. Proline accumulated across T1–T3, supporting osmoprotection and redox buffering [55,56] and aligning with reports under single/combined stresses in cereals [57]. Elevated SP, most pronounced by T3, points to proteome remodeling toward protective functions. Together, these shifts reflect multigenerational priming of defence pathways in wheat [58,59,60]. By contrast, Chl was consistently depressed under dua; stress, in line with combined-stress damage to the photosynthetic apparatus [61,62]. Notably, the Chl penalty narrowed by T3 and Chl-STI rose from STI1 → STI3, suggesting partial transgenerational acclimation potentially via enhanced chloroplastic protection or heat-shock/antioxidant buffering [59,63]. Even so, combined stress still imposed a clear photosynthetic cost each generation. Growth and yield were strongly curtailed under stressed condition, including PH, SL, NSS, NGS, and TKW all fell relative to controls [61,64]. Yet by T3, PH and TKW losses were less severe, consistent with multigenerational adjustments that aid stature maintenance and grain filling [59,65]. In contrast, NGS and NSS remained persistently low across generations, indicating reproductive development (floret fertility, pollen/grain set) is more refractory to priming than vegetative growth or kernel filling [66]. Defence/compatible-solute traits (APX, CAT, SOD, GR, ProC, SP; AsA/GSH) showed rising STI toward T3 (significant pairwise gains), evidencing a strengthening priming effect [59]. Traits that decline under stress (Chl, PH, TKW) displayed improving STI, smaller proportional penalties by T3, whereas NSS/NGS STIs were flat, reinforcing a resilience gap between vegetative/physiological buffering and reproductive success [67]. Overall, repeated drought and heat exposure drive an adaptive re-balancing: defence and osmotic adjustment are up-shifted while photosynthesis and yield are down-shifted, with partial mitigation by T3 [68,69]. These patterns underscore the potential of transgenerational priming to bolster resilience, particularly height and kernel mass, while highlighting persistent constraints at grain set that breeders must address [59,63,64].

Genome-wide association highlights key chromosomal regions

To decipher the genetic basis of these multigenerational stress responses, a GWAS was conducted for various stress tolerance indices and trait values across generations. The GWAS identified numerous significant marker–trait associations distributed across the wheat genome, with certain chromosomal hotspots emerging where multiple traits’ associations co-localized. Three genomic regions stood out for their pleiotropic influence on stress-memory traits. For instance, a cluster of SNPs on 1B was associated with several traits, including Chl_STI2, SP_STI3, and TKW_STI2. All these traits showed their top GWAS signals at virtually the same locus on 1B (around 581.20–581.21 Mb). The fact that one genomic region influences a photosynthetic trait (chlorophyll), a biochemical trait (protein content), and a yield trait (kernel weight) under stress is intriguing and suggests a major pleiotropic QTL for combined stress tolerance. It implies that there may be one or more regulatory genes in this 1B interval that orchestrate broad stress responses affecting metabolism and yield. This co-localization is consistent with earlier studies where major effect loci were found to simultaneously affect multiple drought-response traits for instance, a QTL on wheat 1B was reported to influence both grain weight and leaf senescence under heat stress in another panel [70]. In our study, the lead SNP in this region (e.g.Kukri_c9752_793) had a notable effect size (β) and was consistently significant for different trait indices, highlighting its robustness. Another prominent locus was on chromosome 2A at approximately 686.86 Mb, which emerged as the single strongest association in the entire dataset for PH_STI2. In fact, a SNP here (wsnp_Ex_rep_c66465_64708628) showed an exceptionally low p-value (~ 10⁻¹⁴) for the height index in the second generation, indicating a highly significant effect. Intriguingly, the same 2A locus also showed up in associations for Chl_STI1 and TKW_STI2. This suggests that the 2A region contains gene(s) that impact both plant stature and grain development under stress perhaps master regulators of growth or resource allocation during stress [13]. Recently reported a similar finding: they identified a candidate gene (a Leucine-Rich Repeat protein,TraesCS2A02G432800) on wheat 2A that was linked to all major agronomic traits in plants that experienced drought over three generations. LRR proteins are often involved in stress signaling and development, so it’s possible that our 2A locus could harbor such a gene conferring broad stress-memory effects. The recurrence of a 2A hotspot in different studies underscores the importance of this region for stress resilience. Breeders might focus on this segment of 2A to pyramid tolerance to drought and heat. A third key region was identified on chromosome 2B (~ 106.94 Mb), which was associated with the redox-related traits specifically, GR_STI2 and GSH_STI3 under stress. Both traits mapped to the same lead SNP (RAC875_c35438_474) on 2B, suggesting a locus that governs aspects of the cellular redox state under combined stress. This co-localization makes biological sense: GR is the enzyme that regenerates GSH, so a genetic locus affecting GR might also influence GSH levels. The 2B region could contain a gene involved in antioxidant regulation or glutathione metabolism. Indeed, among the candidate genes in that region, we found a putative pectinesterase and a methyltransferase. Pectinesterases can modulate cell wall properties and also influence the production of methanol and pectic fragments that can act as signals during stress, whereas certain methyltransferases may be involved in epigenetic or metabolic regulation under stress. It’s conceivable that genetic variation in such genes could affect how robustly a plant’s antioxidant system responds (e.g., via signaling or transcriptional control of ROS-scavenging pathways). The identification of this 2B locus dovetails with the idea that controlling oxidative damage is central to surviving combined drought and heat, a concept supported by many studies linking QTLs for drought tolerance with genes in the antioxidant network [71].

Beyond these, other noteworthy associations included signals on group-7 chromosomes related to yield components (for instance, SNPs on 7A were linked with spikelet and grain number indices in some generations) and on group-5 and 6 for certain enzyme indices (e.g., a 5A locus for CAT activity in T3, and a 6B locus for ascorbate content in T2/T3). Many of these associations clustered in LD-defined regions, meaning that within each region, multiple SNPs in linkage disequilibrium all showed signals for the trait, bolstering confidence that a true QTL resides there. In total, the GWAS hits could be summarized into about seven distinct genomic regions that collectively explained much of the multi-trait variance. The fact that multiple traits’ QTL overlapped in the same regions suggests the presence of core stress-regulatory loci. These might correspond to genes that integrate various stress signals and allocate resources (for example, trade-offs between growth and defense, as we observed phenotypically). Our findings are reminiscent of results in barley, where GWAS for combined drought–heat tolerance also found clusters of trait associations in a few key regions [14]. Such convergence is encouraging, as it implies complex stress tolerance is not controlled by completely random, trait-specific loci, but rather by some common “hub” genes or genomic regions. For breeders, markers in these 1B, 2A, 2B (and possibly 6B, 7A) regions could serve as valuable predictors of multi-trait performance under stress. It is worth noting that the allelic effects we detected were often consistent across related traits—for example, the favorable allele at the 1B locus tended to increase chlorophyll and protein (beneficial under stress) and also increase kernel weight, whereas the favorable allele at the 2A locus tended to reduce plant height under stress (dwarfing can be advantageous for stress as it often correlates with earlier maturation or better water use). Such consistent directions suggest these loci might underlie adaptive trait syndromes for stress tolerance.

Candidate genes underlying stress resilience

To probe the biological mechanisms, we examined candidate genes located in the key GWAS regions and also compared their expression in a stress-tolerant versus a stress-susceptible wheat genotype under combined drought–heat conditions. Several plausible candidates emerged that align with known stress response pathways. On chromosome 1B, a Trehalose-6-Phosphate Phosphatase (T6PP) gene (TraesCS1B02G350600) was identified near the peak SNPs. T6PP is involved in the trehalose biosynthetic pathway, converting trehalose-6-phosphate to trehalose. Trehalose and its intermediates act as signaling molecules and osmoprotectants in plants, known to improve stress tolerance by stabilizing proteins and membranes [72]. The presence of a T6PP candidate at a major tolerance locus suggests that enhanced trehalose cycling might be a feature of stress-primed plants. In fact, our expression analysis showed that this T6PP gene was rapidly and strongly induced in the tolerant genotype under combined drought and heat stress; its transcript increased significantly by 12 h into stress and peaked by 36 h in the tolerant line, whereas the sensitive line showed a much smaller increase. This timing indicates that the tolerant plants quickly activate trehalose metabolism, likely to mobilize energy and protect cells as dehydration and heat set in. The early induction of T6PP in tolerant plants could facilitate the accumulation of trehalose or related sugars, contributing to better osmotic adjustment and carbon supply during stress [73]. On chromosome 2A, a notable candidate was a Vacuolar Iron Transporter (VIT) gene(TraesCS2A02G434100/200). Iron is a double-edged sword in stress: it’s essential for photosynthesis and respiration, but free iron can catalyze harmful ROS via Fenton reactions. Under drought and heat, cellular damage often releases free iron, exacerbating oxidative stress. Vacuolar iron transporters help sequester excess iron into vacuoles, reducing cytosolic iron levels and thus limiting ROS generation [74]. The identification of a VIT at the 2A locus (which affected plant height, chlorophyll, and yield) aligns well with a mechanism of improved stress tolerance: genotypes with efficient iron sequestration could protect chloroplasts and other organelles from iron-induced oxidative damage, thereby retaining more chlorophyll and sustaining growth under stress. Our RT-qPCR data support this – the tolerant genotype showed a significant induction of the VIT gene by 24 h into stress (much stronger than in the sensitive genotype). By quickly upregulating iron storage, the tolerant plants may prevent some of the oxidative damage that would otherwise impair photosynthesis and stunt growth during combined drought/heat. This swift activation of iron homeostasis genes in tolerant plants is a hallmark of effective stress response, as also noted by [75], who found that drought-tolerant wheat activated metal transporter genes to mitigate oxidative stress. Several cell wall-related genes were implicated. On 2A, a xyloglucan endotransglucosylase/hydrolase (XTH) gene (TraesCS2A02G434300) was present. XTHs modify xyloglucans in the cell wall, loosening or restructuring the wall, which can be crucial for growth regulation under stress. Meanwhile, on 2B, a pectinesterase gene was a candidate at the GR/GSH locus. Pectinesterases alter the methylation status of pectin in cell walls, which not only affects wall porosity and rigidity but can also release methanol and pectic acid, molecules that can act in stress signaling and ROS modulation. The involvement of these enzymes suggests that cell wall dynamics are part of the plant’s response to combined stress. In a drying soil and high heat, cell walls may need to strengthen in some tissues (to prevent wilting) and loosen in others (to allow continued root growth or floral expansion), and these enzymes mediate such changes. We found that the pectinesterase candidate was upregulated early (by 12 h) in the tolerant genotype, indicating a prompt cell wall adjustment response. Early pectin modification could help maintain cell turgor and perhaps prime the tissue for recovery when stress is alleviated. This rapid response was largely absent in the susceptible genotype, which could contribute to its poorer performance (cells in the sensitive plants may have been less able to adjust their wall mechanics, suffering greater injury). A CCT-domain transcription factor was noted on 2B as well. CCT (CONSTANS, CO-like, and TOC1) domain proteins often relate to photoperiodic flowering and circadian regulation, but some also have roles in stress responses and development timing. A variation in a CCT factor might influence flowering time or developmental rate under stress, which can affect yield (e.g., earlier flowering to escape drought, or sustained growth if conditions allow). Another interesting candidate found (likely on a different chromosome 6B region linked to ascorbate content) was an F-box protein gene. F-box proteins are components of the ubiquitin–proteasome system, targeting specific proteins for degradation, including many signaling proteins [76]. Stress-induced F-box proteins can modulate hormone signals (like ABA or ethylene) and the turnover of misfolded proteins. Identification of an F-box suggests protein quality control or signaling regulation is a factor in stress memory. Additionally, a serine/threonine protein phosphatase in that region implies changes in phosphorylation signaling cascades during stress. These could be involved in ABA signaling or MAPK pathways that are central to drought and heat responses [77]. While the exact function is unclear, it might be an RNA or small-molecule methyltransferase. Either way, such enzymes can be part of epigenetic regulation or metabolic fine-tuning in stress conditions. On 6B, we found a glycosyltransferase candidate enzymes that attach sugar moieties to other molecules (like secondary metabolites, hormones, or even antioxidant flavonoids). Upregulation of a glycosyltransferase in tolerant plants (as we observed by 24–36 h) could indicate the plants are detoxifying reactive compounds or storing signaling molecules in an inactive glycosylated form until needed. For example, glycosylation of ABA or phenolics can modulate their activity during stress. Taken together, the candidate gene analysis paints a picture of stress resilience as a concerted, multi-layered response. Importantly, when we examined the expression of six candidate genes (spanning categories like carbohydrate metabolism, cell wall, metal homeostasis, signaling, epigenetics, and secondary metabolism), we found that the tolerant wheat genotype mounted a faster and stronger induction of all these genes compared to the susceptible genotype. By 12 h into combined drought–heat stress, the tolerant plants had already significantly upregulated genes for trehalose phosphate metabolism (T6PP), cell wall modification (pectinesterase), and an arginine/serine-rich splicing factor (involved in RNA processing), whereas the sensitive plants showed minimal changes by that time [78]. By 24 h, the tolerant genotype showed heightened expression of the vacuolar iron transporter, the methyltransferase, and the glycosyltransferase, none of which were as strongly induced in the susceptible line. This timing corresponds to the phase when stress effects become more pronounced (prolonged water deficit and heat), so the tolerant plants by then are activating metal sequestration, possibly chromatin/histone modification, and scavenging/detoxification of reactive intermediates (via glycosylation), all critical for coping with sustained stress. The susceptible genotype lagged, with some of these defenses only modestly rising later. At 36 h, all examined genes in the tolerant plants reached their peak or remained highly expressed, indicating a full deployment of defense programs, whereas the susceptible plants never attained comparable expression levels [79]. The coordinated induction of genes in the tolerant wheat spans many functional categories from osmolyte production (trehalose via T6PP), ROS/metal detox (VIT, glycosyltransferase), structural fortification (pectinesterase, possibly actin bundling), to regulatory processes (splicing factor, methyltransferase). Such a polygenic response likely provides an integrated protection: osmolytes and trehalose protect proteins and membranes, antioxidants and iron sequestration reduce oxidative damage, cell wall and cytoskeleton adjustments preserve cellular integrity, and epigenetic/transcriptional changes ensure that stress-response genes stay active. This mirrors observations in other studies that found tolerant cultivars exhibit a more rapid and robust transcriptional reprogramming under stress than sensitive cultivars [80,81]. For example, drought-tolerant wheat varieties have been reported to show earlier upregulation of protective genes (like dehydrins, heat-shock proteins, and antioxidant enzymes) compared to drought-sensitive varieties, correlating with less cellular damage [71,82]. Our data specifically highlight early upregulation of a splicing factor and a methyltransferase in the tolerant line, hinting that post-transcriptional and epigenetic regulatory mechanisms are at play in transgenerational stress tolerance. This aligns with the growing evidence that stress memory involves epigenetic marks and altered RNA processing that prime certain genes for quicker activation [83,84]. In conclusion, the results from this multigenerational study demonstrate that wheat’s exposure to combined drought and heat stress can induce heritable enhancements in its defense arsenal, antioxidants, osmolytes, and stress proteins while partially preserving growth and yield functions in later generations. These effects are underpinned by complex genetic loci, some of which control multiple facets of the stress response.

Transgenerational priming appears to rewire the plant’s regulatory networks so that, when the stress strikes again, tolerant progeny launch a suite of protective measures (from metabolic adjustments to gene expression shifts) more rapidly and intensely than unprimed plants. The outcome is a plant better “prepared” for stress: it suffers damage (chlorophyll loss, yield reduction) but to a lesser extent, and it invests heavily in survival tools (ROS scavengers, osmoprotectants). From an applied perspective, understanding these mechanisms and associated loci could guide breeding programs aimed at climate-resilient crops. For instance, markers on 1B or 2A linked to high antioxidant capacity and maintained grain weight under stress could be used to select for wheat lines with inherent stress memory. Additionally, the candidate genes identified, such as those in trehalose metabolism, iron homeostasis, and epigenetic regulation, provide potential targets for genetic engineering or gene editing. Introducing alleles that confer stronger or quicker induction of these protective genes might endow crops with improved tolerance to the increasingly frequent combination of drought and heat episodes. Finally, our findings reinforce a broader biological concept: plants can “remember” past stresses, and this memory can be encoded in their genes and even passed to subsequent generations. While transgenerational stress memory in crops is a relatively new frontier, evidence is accumulating that mild stresses or priming in one generation can lead to epigenetic modifications (DNA methylation, histone changes, small RNAs) that alter how progeny respond to stress. Our multi-cycle experiment under combined drought and heat extends this paradigm, showing that even extremely harsh stress combinations can induce long-lasting adaptive responses. This offers a hopeful perspective that, through either natural acclimation or human-mediated selection, crops might be equipped to cope with complex stress environments in a future of climate uncertainty. The trade-off, however, is that such memory may come at a cost to maximal productivity under non-stress conditions, as resources are reallocated. Future research should thus aim to uncouple the positive protective effects from growth penalties, ensuring that stress memory enhances resilience without significantly diminishing yield potential in favorable conditions.