Heat Stress-Tolerant Quantitative Trait Loci Identified Using Backcrossed Recombinant Inbred Lines Derived from Intra-Specifically Diverse Aegilops tauschii Accessions

In the face of climate change, bringing more useful alleles and genes from wild relatives of wheat is crucial to develop climate-resilient varieties. We used two populations of backcrossed recombinant inbred lines (BIL1 and BIL2), developed by crossing and backcrossing two intra-specifically diverse Aegilops tauschii accessions from lineage 1 and lineage 2, respectively, with the common wheat cultivar ‘Norin 61′. This study aimed to identify quantitative trait loci (QTLs) associated with heat stress (HS) tolerance. The two BILs were evaluated under heat stress environments in Sudan for phenology, plant height (PH), grain yield (GY), biomass (BIO), harvest index (HI), and thousand-kernel weight (TKW). Grain yield was significantly correlated with BIO and TKW under HS; therefore, the stress tolerance index (STI) was calculated for these traits as well as for GY. A total of 16 heat-tolerant lines were identified based on GY and STI-GY. The QTL analysis performed using inclusive composite interval mapping identified a total of 40 QTLs in BIL1 and 153 QTLs in BIL2 across all environments. We detected 39 QTLs associated with GY-STI, BIO-STI, and TKW-STI in both populations (14 in BIL1 and 25 in BIL2). The QTLs associated with STI were detected on chromosomes 1A, 3A, 5A, 2B, 4B, and all the D-subgenomes. We found that QTLs were detected only under HS for GY on chromosome 5A, TKW on 3B and 5B, PH on 3B and 4B, and grain filling duration on 2B. The higher number of QTLs identified in BIL2 for heat stress tolerance suggests the importance of assessing the effects of intraspecific variation of Ae. tauschii in wheat breeding as it could modulate the heat stress responses/adaptation. Our study provides useful genetic resources for uncovering heat-tolerant QTLs for wheat improvement for heat stress environments.


Introduction
With the rising global demand for food and frequent occurrences of abiotic stresses, such as heat stress, it is crucial to develop wheat varieties with high yield potential and heat tolerance [1,2].Therefore, breeding for heat stress tolerance is essential to enhance wheat grain yield and adaptation.Typically, crop breeders utilize genetic variation to improve crops against environmental stress.However, due to intensive selection for grain yield and other desirable traits, most wheat varieties exhibit narrow genetic diversity, which limits their potential for further improvement of heat stress tolerance [3].Exploring and utilizing the adaptive potential of wild progenitors in a systematic approach is a powerful strategy Plants 2024, 13, 347 2 of 18 and a practical solution to identify promising material and broaden the genetic base of wheat [4].
Wild relatives of wheat, including the Aegilops species, are valuable resources for developing new genetic materials.Numerous studies have reported the tolerance of the species to abiotic stresses [5][6][7].Among the Aegilops species, Aegilops tauschii is the most promising due to the similarity of its D genome to that of bread wheat.No special cytological technique is required to induce meiotic recombination [6,8,9].
Globally, attempts have been made to broaden the genetic diversity of wheat utilizing Ae. tauschii.Several QTLs and marker-trait associations (MTAs) have been identified for abiotic stress-adaptive traits [10][11][12][13].However, a limited number and diversity of Ae. tauschii has been explored and utilized in the development of synthetic wheat.To ensure systematic exploration and extensive utilization of the tremendous diversity of Ae. tauschii in the breeding programs, multiple synthetic derivatives (MSDs) have been developed, utilizing 43 Ae.tauschii accessions [8,14].The MSD population was systematically evaluated to identify useful stress-adaptive traits, MTAs, and QTLs.The MSD population displayed significant genetic diversity when exposed to heat, drought, and heat-drought combinations under field conditions in Sudan.Several MTAs were identified to be associated with heat, drought, heat-drought, and grain characteristics and quality [15][16][17][18][19][20].
The MSDs, as a mixture population of several BC 1 F 7 lines derived from crosses of 43 Ae.tauschii accessions, are expected to have high linkage disequilibrium; therefore, the previously identified MTAs would not be useful for direct use in breeding unless validated.A better approach for efficient and precise mapping and validation of the QTLs is to reduce the genomic contribution of Ae. tauschii in the progeny by utilizing one or more backcrosses and developing backcross-derived inbred lines (BILs).Therefore, we developed BILs targeting traits associated with heat stress tolerance to be evaluated in field experiments under heat and normal conditions and select lines fixed for desired donor alleles of QTLs.
We selected two MSD lines previously identified by Elbashir et al. [21] based on their high heat tolerance efficiency.Two backcross recombinant inbred lines (BILs) were developed utilizing the primary synthetic derivatives of the two MSD lines crossed and backcrossed to 'Norin 61' (N61).Then, crosses were advanced to BC 1 F 5 using the singleseed descent method.
The two BILs were used in this study to identify QTLs associated with heat stress tolerance under field conditions.We evaluated the heat stress response of the two BIL populations (107 lines from BIL1 and 164 lines from BIL2) across four environments in Sudan with temperature gradients ranging from relatively cool temperature in the north (Dongola), to a continuous heat stress condition in the central region (Wad Medani).Utilizing genotyping by random amplicon sequencing direct (GRAS-Di) markers, we identified novel QTLs for heat stress tolerance and verified the MTAs previously reported using the MSD population.Our study offers novel genetic resources and QTLs for breeding wheat with enhanced tolerance to heat stress conditions.

Climate Conditions during the Growing Seasons
Temperature data were recorded in season 2020/2021 at Wad Medani (WM1) and in season 2021/2022 at Dongola (DN), Waha (WA), and Wad Medani (WM2).The mean maximum and minimum temperatures for each environment during the growing seasons were 30.0 and 11.6 • C at DN, 32.9 and 16.4 • C at WA, 36.3 and 17.9 • C at WM1, and 35.7 and 15.9 • C at WM2, respectively (Figure 1a-d).In general, DN was the coolest among the four environments during the growing season (Figure 1a), whereas WM1 was the hottest.During the grain-filling period, the maximum temperature ranged from 21.5 to 36.8 • C at DN, 25.0 to 43.0 • C at WA, 27.0 to 44.0 • C at WM1, and 27.0 to 42.0 • C at WM2.At DN, maximum temperatures ≥40.0 • C during the grain-filling period were not reported; however, at WA, WM1, and WM2, 7, 19, and 11 days were reported with maximum temperatures ≥40.0 • C, respectively.Likewise, the number of days with mean temperatures ≥20.0 • C during the grain filling period was 33, 57, 75, and 67 days at DN, WA, WM1, and WM2, respectively.
Plants 2024, 13, x FOR PEER REVIEW 3 of 19 maximum temperatures ≥ 40.0 °C during the grain-filling period were not reported; however, at WA, WM1, and WM2, 7, 19, and 11 days were reported with maximum temperatures ≥ 40.0 °C, respectively.Likewise, the number of days with mean temperatures ≥ 20.0 °C during the grain filling period was 33, 57, 75, and 67 days at DN, WA, WM1, and WM2, respectively.

Impact of Heat Stress on the BIL Populations
The analysis of variance (ANOVA) revealed significant effects of genotype, environment, and their interaction on most traits in both BILs, with the exception of genotypic effect on plant height (PH) and biomass (BIO) in BIL1, environmental effect on thousand kernel weight (TKW) in BIL1, and the genotype and environment (G × E) interaction on PH in BIL2 (Table S1).
High-to-moderate broad-sense heritability (h 2 ) estimates were found for GY, DH, and DM in both BILs (Table S1).Moderate-to-low h 2 estimates were found for GFD, BIO, and HI in both BILs.The h 2 estimated for PH and TKW varied between the two BILs (Table S1).
We calculated the two stress tolerance indices (STIs) for GY.The STI1-GY was calculated taking DN as the non-stress environment, whereas the STI2-GY was calculated considering WM2 as a non-stress environment.Then, we performed regression analysis between the GY at DN and the STI1-GY at WA, and the STI1-GY at WM1 (Figure 2a,b).We also performed a regression analysis between the GY at WM2 and the STI2-GY at WM1

Impact of Heat Stress on the BIL Populations
The analysis of variance (ANOVA) revealed significant effects of genotype, environment, and their interaction on most traits in both BILs, with the exception of genotypic effect on plant height (PH) and biomass (BIO) in BIL1, environmental effect on thousand kernel weight (TKW) in BIL1, and the genotype and environment (G × E) interaction on PH in BIL2 (Table S1).
High-to-moderate broad-sense heritability (h 2 ) estimates were found for GY, DH, and DM in both BILs (Table S1).Moderate-to-low h 2 estimates were found for GFD, BIO, and HI in both BILs.The h 2 estimated for PH and TKW varied between the two BILs (Table S1).
We calculated the two stress tolerance indices (STIs) for GY.The STI1-GY was calculated taking DN as the non-stress environment, whereas the STI2-GY was calculated considering WM2 as a non-stress environment.Then, we performed regression analysis between the GY at DN and the STI1-GY at WA, and the STI1-GY at WM1 (Figure 2a,b).We also performed a regression analysis between the GY at WM2 and the STI2-GY at WM1 (Figure 2c).In BIL1, six genotypes exhibited higher STI values than the recurrent parent, N61, in both WA and WM1 (Figure 2a,b).In BIL2, five genotypes exhibited higher STI and GY values than N61 at WA.At WM1, three genotypes exhibited higher STI and GY values than N61 (Figure 2a,b).For STI2, sixteen and seven genotypes in BIL1 and BIL2, respectively, exhibited higher STI2 values relative to N61 (Figure 2c).Three genotypes in both BILs showed higher GY and STI2 values than N61 (Figure 2c).
(Figure 2c).In BIL1, six genotypes exhibited higher STI values than the recurrent parent, N61, in both WA and WM1 (Figure 2a,b).In BIL2, five genotypes exhibited higher STI and GY values than N61 at WA.At WM1, three genotypes exhibited higher STI and GY values than N61 (Figure 2a,b).For STI2, sixteen and seven genotypes in BIL1 and BIL2, respectively, exhibited higher STI2 values relative to N61 (Figure 2c).Three genotypes in both BILs showed higher GY and STI2 values than N61 (Figure 2c).Comparison of grain yield at Dongola (DN) for the two populations (BIL1 and BIL2) using the stress tolerance index (STI).The STI-GY was calculated twice.First, STI1-GY was calculated regarding the GY values at DN as the non-stress environment and WA or WM1 as the hot environment.Second, STI2-GY was calculated considering the GY of WM2 as the non-stress environment and WM1 as the stress environment.The GY at DN was regressed against STI1-GY at (a) WA and (b) WM1.The GY at WM2 was regressed against STI2-GY at (c) WM1.The red circle and diamond represent the recurrent parent N61 in BIL1 and BIL2, respectively.The horizontal dashed red line is plotted based on the GY mean plus the Least Significant Difference (LSD 0.05) of both BIL1 and BIL2.Similarly, the horizontal dashed blue line corresponds to the GY mean plus the LSD of BIL1, and the dashed black line reflects the GY mean plus the LSD of BIL2.The vertical blue line is calculated based on the STI mean plus the LSD of BIL1, while the vertical black line is calculated based on the STI mean plus the LSD of BIL2.The horizontal blue arrows represent the heat tolerance ranges (HTR) for both BILs.

Relationship among Traits
In both BILs, GY significantly and consistently correlated (p < 0.05) with BIO and HI in all environments (Table S2).However, at DN, GY significantly correlated with seed number per spike (SN) and TKW in BIL1 and BIL2, respectively.At WM1, GY significantly correlated with SN in BIL1 and BIL2, whereas at WM2, GY correlated only with TKW in BIL2.The STI1-GY at WA and WM1 significantly correlated with BIO and HI in both BILs.Likewise, the STI2-GY at WM1 significantly correlated with SN, BIO, and HI in both BILs (Table S2).Comparison of grain yield at Dongola (DN) for the two populations (BIL1 and BIL2) using the stress tolerance index (STI).The STI-GY was calculated twice.First, STI1-GY was calculated regarding the GY values at DN as the non-stress environment and WA or WM1 as the hot environment.Second, STI2-GY was calculated considering the GY of WM2 as the non-stress environment and WM1 as the stress environment.The GY at DN was regressed against STI1-GY at (a) WA and (b) WM1.The GY at WM2 was regressed against STI2-GY at (c) WM1.The red circle and diamond represent the recurrent parent N61 in BIL1 and BIL2, respectively.The horizontal dashed red line is plotted based on the GY mean plus the Least Significant Difference (LSD 0.05 ) of both BIL1 and BIL2.Similarly, the horizontal dashed blue line corresponds to the GY mean plus the LSD of BIL1, and the dashed black line reflects the GY mean plus the LSD of BIL2.The vertical blue line is calculated based on the STI mean plus the LSD of BIL1, while the vertical black line is calculated based on the STI mean plus the LSD of BIL2.The horizontal blue arrows represent the heat tolerance ranges (HTR) for both BILs.

Relationship among Traits
In both BILs, GY significantly and consistently correlated (p < 0.05) with BIO and HI in all environments (Table S2).However, at DN, GY significantly correlated with seed number per spike (SN) and TKW in BIL1 and BIL2, respectively.At WM1, GY significantly correlated with SN in BIL1 and BIL2, whereas at WM2, GY correlated only with TKW in BIL2.The STI1-GY at WA and WM1 significantly correlated with BIO and HI in both BILs.Likewise, the STI2-GY at WM1 significantly correlated with SN, BIO, and HI in both BILs (Table S2).

Linkage Maps for the BILs
Both BIL1 and BIL2 were genotyped using genotyping by random amplicon sequencing direct (GRAS-Di).The details of the linkage map construction of BIL1 were mentioned in Ahmed et al. [22].Briefly, the high-density linkage map was developed using 2882 markers.The markers were distributed unevenly across multiple chromosomes and subgenomes.The D-subgenome recorded the highest marker density followed by B-and A-subgenomes.Chromosome 3D displayed the highest number of markers, whereas chromosome 6B showed the lowest number of markers.
In BIL2, 19,765 markers were used for GRAS-Di genotyping, with 6504 exhibiting polymorphisms between the synthetic wheat donor parent, Syn44, and the recurrent parent, N61.Of these polymorphic markers, 3404 (52.3%) were of high quality, with an average of 162 markers per chromosome.A high-density linkage map was constructed utilizing the 3404 markers spread across 21 linkage groups, covering a total genetic distance of 5673.33 cM (Figure 3 and Table S3).The average distance per chromosome was 270.16 cM.The markers were not uniformly distributed among the chromosomes and subgenomes.Most markers were mapped to the B-(1207, 35.46%) and D-(1201, 35.28%) subgenomes, which had total genetic lengths of 1978.99 and 1943.64 cM, respectively.A total of 996 (29.26%) markers were mapped to the A-subgenome with a total genetic length of 1750.7 cM (Table S3).The B-and D-subgenomes had the highest marker density, with one marker per 1.7 cM, whereas the A-subgenome had one marker per 2.0 cM (Table S3).Chromosome 7D had the highest number of markers (264) with a genetic distance of 435.36 cM, whereas chromosome 5A had the lowest number of markers (69) with a genetic distance of 207.37 cM.Chromosomes 2A, 1B, and 2B had marker gaps greater than 30 cM (Table S4).

Linkage Maps for the BILs
Both BIL1 and BIL2 were genotyped using genotyping by random amplicon sequencing direct (GRAS-Di).The details of the linkage map construction of BIL1 were mentioned in Ahmed et al. [22].Briefly, the high-density linkage map was developed using 2882 markers.The markers were distributed unevenly across multiple chromosomes and subgenomes.The D-subgenome recorded the highest marker density followed by B-and Asubgenomes.Chromosome 3D displayed the highest number of markers, whereas chromosome 6B showed the lowest number of markers.
In BIL2, 19,765 markers were used for GRAS-Di genotyping, with 6504 exhibiting polymorphisms between the synthetic wheat donor parent, Syn44, and the recurrent parent, N61.Of these polymorphic markers, 3404 (52.3%) were of high quality, with an average of 162 markers per chromosome.A high-density linkage map was constructed utilizing the 3404 markers spread across 21 linkage groups, covering a total genetic distance of 5673.33 cM (Figure 3 and Table S3).The average distance per chromosome was 270.16 cM.The markers were not uniformly distributed among the chromosomes and subgenomes.Most markers were mapped to the B-(1207, 35.46%) and D-(1201, 35.28%) subgenomes, which had total genetic lengths of 1978.99 and 1943.64 cM, respectively.A total of 996 (29.26%) markers were mapped to the A-subgenome with a total genetic length of 1750.7 cM (Table S3).The B-and D-subgenomes had the highest marker density, with one marker per 1.7 cM, whereas the A-subgenome had one marker per 2.0 cM (Table S3).Chromosome 7D had the highest number of markers (264) with a genetic distance of 435.36 cM, whereas chromosome 5A had the lowest number of markers (69) with a genetic distance of 207.37 cM.Chromosomes 2A, 1B, and 2B had marker gaps greater than 30 cM (Table S4).

Identified QTLs in All Environments
QTL analysis was performed in both BILs using DH, DM, GFD, PH, GY, BIO, TKW, HI, and SN.In addition, stress tolerance indices for GY, BIO, and TKW and their Best Linear Unbiased Predictor (BLUP) were used.Using inclusive composite interval mapping of QTLs with additive and dominance effect analysis (ICIM-ADD), we identified 40 QTLs in BIL1 for the studied traits in the four environments.The identified QTLs were mapped on all chromosomes except 3A, 4A, 5A, 3B, 4B, and 6B.The LOD scores of the identified QTLs ranged from 2.50 to 5.22, and the phenotypic variation explained ranged from 5.14 to 15.43%.The numbers of QTLs detected at WM1, DN, WA, and WM2 were seventeen, ten, nine, and four, respectively.The highest number of QTLs was identified in the D-subgenome (17, 42.5%) for DH, GFD, GY, HI, PH, STI1-GY, STI1-TKW, STI2-TKW, Figure 3. Genetic maps constructed using 3404 GRS-Di markers in BIL2.For BIL1, refer to Ahmed et al. [22].The positions of the QTLs are colored.

QTLs Associated with Heat Stress Response in Both BILs
In BIL1, 14 QTLs associated with stress tolerance indices for GY, BIO, and TKW were detected on chromosomes 1A, 2B, 1D, 5D, and 7D (Table 1, Figure 4).Three QTLs for STI1-GY were only detected in WA on chromosomes 1B, 1D, and 6D, explaining a phenotypic variation of 7.50, 15.43, and 8.67, respectively.A single QTL associated with STI1-BIO was detected in WM1 on chromosome 1A, explaining 5.14% of the phenotypic variation.For STI1-TKW and STI2-TKW, 10 QTLs were detected in both heat stress environments (WA and WM1) on chromosomes 1A, 2B, 1D, 5D, and 7D, with phenotypic variation ranging from 8.29 to 14.54%.Of these ten QTLs, two on chromosome 1A and one on chromosome 2B were stable and overlapped across the two environments.One of the two stable QTLs of STI1-TKW on chromosome 1A was detected at 173 cM between the flanking markers AMP0035547 and AMP0004300, and the phenotypic variations explained were 9.89 and 14.54% at WA and WM1, respectively.The other QTL was detected at 118-119 cM between the flanking markers AMP0036610 and AMP0034796 and the phenotypic variations explained were 11.37 and 14.37% at WA and WM1, respectively.Meanwhile, the QTL of STI1-TKW on chromosome 2B was detected at 62 and 63 cM between the flanking markers AMP0009891 and AMP0006464 and explained 8.66 and 10.41% of the phenotypic variation at WA and WM1, respectively (Table 1, Figure 4).
In BIL2, 25 QTLs associated with stress tolerance indices for GY, BIO, and TKW were identified on chromosomes 3A, 5A, 2B, 4B, and all D-subgenome chromosomes, except for chromosome 1D.Six QTLs associated with STI1 or STI2 for GY were found on chromosomes 3A, 5A, 4B, and 6D.Among the six QTLs, a QTL was stable and consistently identified at WA and WM1 at 88 cM on chromosome 5A between the flanking markers AMP0011577 and AMP0030240, explaining 10.49 and 13.70% of the phenotypic variation, respectively (Tables 1 and 2).For STI-BIO and STI2-BIO, 12 QTLs were detected on chromosomes 3A, 5A, 2B, 3D, 4D, 5D, and 7D.Among them, one QTL on chromosome 5A was stable.The stable QTL was detected at 3 cM between flanking markers AMP0003832 and AMP0029058 and explained 3.10 and 6.20% of phenotypic variations.Seven QTLs were detected at WM1 associated with STI1 and ST2 for TKW on chromosomes 3A, 2D, 4D, and 6D (Table 1).
Chromosomes 1A, 2A, 7A, 2B, 5B, and 7B, as well as all D-subgenome chromosomes except 7D, were common regions that harbored different QTLs in both BILs.BIL1 exhibited a higher number of QTLs on chromosomes 1A, 2B, and 1D, whereas BIL2 had a greater number of QTLs on all D-subgenome chromosomes, except 1D and 7D, as well as 7A and 7B (Figure 4).
Interestingly, some of the identified QTLs in the two BILs co-located in the same region of the chromosome with some MTAs identified in previous studies conducted using the related population of the MSDs (Figure 4).In addition, most of the co-located QTLs in both BILs were identified in the D-subgenome with the chromosome 5D contributing the most (Figure 4).In BIL2, 25 QTLs associated with stress tolerance indices for GY, BIO, and TKW were identified on chromosomes 3A, 5A, 2B, 4B, and all D-subgenome chromosomes, except for chromosome 1D.Six QTLs associated with STI1 or STI2 for GY were found on chromosomes 3A, 5A, 4B, and 6D.Among the six QTLs, a QTL was stable and consistently identified at WA and WM1 at 88 cM on chromosome 5A between the flanking markers AMP0011577 and AMP0030240, explaining 10.49 and 13.70% of the phenotypic variation, respectively (Tables 1 and 2).For STI-BIO and STI2-BIO, 12 QTLs were detected on chromosomes 3A, 5A, 2B, 3D, 4D, 5D, and 7D.Among them, one QTL on chromosome 5A was

Discussion
In the face of climate change, evaluating wheat lines with limited genetic diversity or with similar genetic backgrounds to the elite breeding lines would not be useful for delivering superior climate-smart varieties.Broadening genetic diversity and bringing more useful alleles and genes from wild relatives of wheat has proven to be successful [6,7,10,27,28].Therefore, a population with a wide diverse genetic background has been developed utilizing 43 Ae.tauschii accessions and named multiple synthetic derivatives (MSDs) [8,14].The MSD population was evaluated over several seasons in multiple environments ranging from normal to continuous heat stress conditions.The results showed that a number of superior lines could be identified under various conditions and for different climateresilient traits [15][16][17][18]21].Moreover, GWAS analysis identified several MTAs under heat, drought, and heat-drought conditions [15][16][17]19,20].Following the identification of several heat-tolerant lines from the MSD population, backcross inbred lines (BILs) were developed in the background of those heat-tolerant lines with exotic alleles from two Ae.tauschii accessions collected from two different locations.
In this study, we evaluated two BILs comprising 271 lines in four environments ranging from relatively cooler environments at Dongola in northern Sudan to continuous heat stress environments at Wad Medani in central Sudan.The analysis of variance showed that the lines of the two BILs exhibited significant variation in grain yield and some agronomic traits under the relatively cool and continuous heat stress conditions, confirming those previously reported for the MSD population and other bread wheat genotypes carrying exotic alleles [10,15,21].The great variation in key traits such as grain yield coupled with a high broad-sense heritability demonstrated the possibility of exploring them in breeding for heat stress tolerance.
In this study, heat stress reduced grain yield compared to the relatively cooler environment at DN, especially at WA and WM1.Compared to DN, the mean grain yield in BIL1 reduced by about 58% at both WA and WM1.Similarly, in BIL2, reductions in GY of 52 and 56% were found at WA and WM1, respectively, compared to that at DN.However, the mean grain yields of the two BILs at WM2 were always comparable with those at DN despite the fact that significant reductions in the grain yield of some individual lines were observed.Figure 1 shows that the two BILs had longer grain-filling periods at the two WM environments than DN, although the heading time of the two BILs was earlier at the WM environments.This might be due to the significant G × E interactions observed for different traits, especially DH, DM, GFD, GY, and BIO (Table S1).For instance, about 72% of the BIL1 lines that showed a reduction in grain yield at WM2 compared to DN were either early (with DH less than 60 days) or late (DH more than 70 days).This might be related to the vernalization gene doses in these populations, which deserve further investigation.
Compared to the recurrent parent, N61, several lines in both BILs showed better performance under both stressed and non-stressed environments.In both BILs at the heat stress environment of WM1, twelve lines exhibited higher STI-GY, out of which seven lines had higher grain yield than N61.Similarly, at WA, sixteen lines showed higher STI-GY than N61, of which eight lines also had higher grain yield.The high grain yield and stress tolerance index were closely associated with the accumulation of high biomass as indicated by the strong correlation of grain yield and STI-GY with the biomass at both heat-stressed environments.High biomass accumulation was found to be important under the same heat stress environments [29], as well as when exotic-derived lines from Ae. tauschii were compared under heat stress with elite wheat lines [10,30].The high-yielding heat-tolerant lines identified in both BILs are promising and could be integrated into wheat heat stress tolerance breeding programs.

The High-Resolution Linkage Maps
The GRAS-Di platform generates many chromosome-spanning genetic markers, allowing the construction of high-resolution linkage maps [31].We constructed two genetic maps using 2882 and 3404 GRAS-Di markers in BIL1 and BIL2, respectively, that were polymorphic between the synthetic parents and N61.The map lengths in both BILs were relatively long, possibly because of the high density of GRAS-Di markers, which generally increases the total length of the linkage map [31].Nevertheless, the length of our maps was comparable to those reported in wheat [32,33].The D-subgenome was longer in BIL1 than that of the A-and B-subgenomes [22].However, in BIL2, the D-subgenome map was similar in length to that of the B-subgenome.The similar map length of the B-and D-subgenomes in BIL2 might indicate higher polymorphism between N61 and the synthetic donor of BIL2 (Syn44) than that of BIL1 (Syn32).Additionally, the higher number of high-quality markers used for map construction in BIL2 than in BIL1 suggests that the D-subgenome of Syn44 exhibited higher interaction with the A-and B-subgenomes of N61.The D-subgenome has been reported previously as the shortest with the lowest number of markers [34,35] due to its low level of polymorphism [36].Unlike these reports, in this study, a higher polymorphism was recorded in the D-subgenome, which might be attributed to its origin from the wild D-genome introduced from Ae. tauschii.Therefore, the two BILs used here could provide a unique opportunity to study the D-subgenome and could be a valuable breeding material for climate change-resilient wheat improvement.

QTLs Identified in All Environments
QTL analysis performed using the two BILs identified 193 QTLs on all chromosomes except 3A, 4A, 5A, 3B, 4B, and 6B under stress and non-stress environments.Although QTLs were detected for all traits studied, the highest number of QTLs was associated with TKW and PH (26 QTLs each), followed by the phenological traits.A total of nine QTLs for grain yield were identified in both BILs.A total of 156 QTLs (80.8%) were identified under the heat stress environments of WA, WM1, and WM2.The high number of QTLs identified under heat stress could be a good indicator of the usefulness of such populations in mining QTLs associated with stress-adaptive traits, as has been previously reported [15,37,38].The D-subgenome contributed 36.9% of the total QTLs identified, including five out of nine GY QTLs and thirteen out of forty QTLs associated with the stress tolerance index.Our results confirmed the importance of the D-subgenome derived from Ae. tauschii for mining stress-adaptive traits [18].

Identification of Stable Major QTLs for Yield-and Heat Stress Tolerance-Related Traits
The present study identified a total of 39 QTLs related to heat stress tolerance in both BILs, of which 25 were from BIL2.The higher number of QTLs identified in BIL2 for heat stress tolerance suggests that BIL2 has greater heat adaptability than BIL1.When the results of this study are linked with that reported by Mahjoob et al. [44] on leaf hair density in two intra-specifically diverse Ae. tauschii accessions, it is clear that different traits have evolved independently in different lineages of the species.Apparently, KU-2124, collected from TauL2 in Iran, contributed more heat stress adaptation traits than KU-2039, collected from TauL1 in Afghanistan.Thus, to further improve the resilience of wheat to climate change, the intraspecific diversity and lineage differences of Ae. tauschii should be explored and exploited.
Out of the 39 heat stress-related QTLs, 12 were found to be major QTLs due to their high contribution to phenotypic variation (10.1-15.5%).The 12 major QTLs were detected on chromosomes 1A, 3A, 5A, 2B, 1D and 5D.Out of the twelve QTLs, eight were for STI-TKW on chromosomes 1A, 3A, 2B, 1D, and 5D.Three of these QTLs were stable across two to three environments and in the BLUP.Additionally, four QTLs were found for STI-GY on chromosomes 3A, 5A, 1D, and 6D with one stable across two environments (Tables 2 and S5).Following the validation of these major stable QTLs, they could be employed to bolster the heat tolerance of wheat with the aid of marker assistance selection.

Common and Specific Regions of Detected QTLs in BIL1 and BIL2
The QTLs identified under heat stress conditions on chromosomes 7A and 7B for DH, DM, GFD, GY, and TKW were specific to BIL2.We found a polytropic effect for DM, DH, and PH under heat stress on chromosome 7A, and for DM and TKW on 7B.A QTL for SSI-KW was co-localized with QTLs of TKW on chromosome 7B [40]; however, a QTL associated with heat stress was not previously reported on chromosome 7A.Thus, these regions in 7A and 7B may be necessary for improving DH, DM, PH, and TKW for heat stress adaption.
The uniqueness of BIL populations used in this study was demonstrated by QTLs identified on chromosomes 3D, 4D, 5D, and 6D.We also found a pleiotropic effect between several important traits in the four D-chromosomes.To our knowledge, no heat stressrelated QTLs have been reported in these regions of the chromosomes.Previously, by utilizing the MSD population, several MTAs have been identified in these regions [15][16][17], which were colocalized with the QTLs identified in this study.This implies that QTLs in these regions are novel and could be specific to the exotic alleles introduced from Ae. tauschii in the MSD and BIL populations.
On chromosome 6D, we identified two QTLs for STI-GY in BIL1 and BIL2 at 53 and 92 cM, respectively.The close position of the two QTLs suggests that they might be the same.No QTL associated with STI has been reported in this region; thus, it is mostly a novel QTL.Following validation of this QTL within the identified hot-spot regions, it could be utilized to enhance the heat tolerance of wheat via marker-assisted selection (MAS).
The pleiotropic and new QTLs identified on chromosomes 7A, 7B, 3D, 4D, 5D, and 6D suggest that these regions are valuable for heat stress tolerance-related traits.Therefore, they could be targeted and prioritized for mining heat stress tolerance alleles and genes and marker-assisted selection.Utilizing the MSD population, we demonstrate the efficacy of the systematic approach we employed in establishing the BIL population.
In this study, we conducted QTL analysis for agronomic traits under heat stress and non-heat stress conditions using two BIL populations.We identified 39 novel QTLs associated with GY, BIO, and TKW stress tolerance indices.More importantly, we identified nine novel STI-GY QTLs on chromosomes 3A, 5A,1B, 4B, 1D, and 6D.Notably, we found a critical region in chromosome 6D where we identified an STI-GY QTL in both BILs, located only 39 cM apart.This region might be of high significance and deserves attention for potential targeted marker development.We also discovered new QTLs associated with important agronomic traits, including grain yield and thousand-kernel weight under heat stress conditions.These QTLs could be employed for marker-assisted selection and gene discovery after validation and can be given priority in breeding for heat stress tolerance.

Plant Materials
In this study, we used two populations of backcrossed recombinant inbred lines (BILs, BC 1 F 5 ).Through intensive evaluation of multiple synthetic derivatives (MSDs), developed utilizing 43 accessions of Ae. tauschii and the durum wheat cultivar 'Langdon', converted to kg ha −1 for further analysis.TKW (g) and SN were determined from random samples of 10 spikes.HI was measured as the ratio of BIO to GY (GY/BIO × 100).
To identify heat stress-tolerant genotypes, the stress tolerance index (STI) for GY (STI-GY), BIO (STI-BIO), and TKW (STI-TKW) was calculated according to Fernandez [47] as where Y N is GY under non-stress cconditions Y S is GY under heat stress conditions, and Y N is the mean GY under non-stress conditions.For BIO and TKW, GY was replaced with BIO and TKW, respectively.We calculated STI-GY twice.First, STI1-GY was calculated regarding the GY values at DN as the non-stress environment and WA or WM1 as the hot environment.Second, STI2-GY was calculated considering the GY of WM2 as the non-stress environment and WM1 as the stress environment.

Statistical Analysis of Phenotypic Data
Analysis of variance (ANOVA) for all traits was performed in Plant Breeding Tools v. 1.4.2[48].The genotype (G), environment (E), and (G × E) interaction were considered random effects.For each trait, the Best Linear Unbiased Predictor (BLUP) was computed and subsequently employed for additional analysis.We used Tukey's honestly significant difference (HSD) test for between-environment comparisons and the t-test for between-BIL comparisons.Pearson's correlation coefficient among traits in each environment was calculated using IBM SPSS Statistics for Windows v. 28.0.1.1(15)[49].Broad-sense heritability was estimated in Plant Breeding Tools.Total genomic DNA was extracted from 2-week-old leaves of BILs using the CTAB method [50].The DNA samples (20 µL; 50-100 ng µL −1 ) were sent to Eurofins Genomics Company, Tokyo, Japan (https://eurofinsgenomics.jp/jp (accessed on 20 October 2023)), for a whole-genome scan with genotyping by random amplicon sequencing direct (GRAS-Di) markers.

Maps Construction
The linkage map for both BILs was constructed using GRAS-Di markers.The construction of GRAS-Di libraries followed the protocol outlined in Hosoya et al. [51].The libraries were sequenced using the Illumina HiSeq series.GRAS-Di software (TOYOTA, Aichi, Japan, https://eurofinsgenomics.jp/jp/service/ngs/gras-di/ (accessed on 21 January 2024)), which is commercially accessible, was used for genotyping.The software assesses the quality of markers by applying empirical criteria for genotyping reproducibility [31].This assessment is based on the number of reads and the consistency of genotyping results across samples (presence and absence of the reads).To determine reproducibility, the software uses trial data from the GRAS-Di database for crop species (Patent ID P2018-42548A).The markers are then ranked using A, B, C, D, and E labels, with A representing the highest reproducibility (≥99.99%),B reproducibility between 99.98% and 99.99%, C reproducibility between 99.90% and 99.98%, and D reproducibility between 99.80% and 99.90%.The A-, B-, C-, and D-ranked markers are of high quality to be used for genotyping.
Because the E-ranked markers contain tested samples with missing values, they are not recommended to be used for genotyping [31].
Previously, we genotyped the BIL1 population using GRAS-Di markers and used this material for QTL analysis of seed dormancy [22].In this study, we genotyped BIL2 using the same method as BIL1.Therefore, both BILs' genetic maps were utilized in the current study.The details of BIL1 map construction have been mentioned in Ahmed et al. [22].
In the same way, the 164 lines of BIL2 were genotyped with 19,765 GRAS-Di markers.We removed markers amplified in all samples from all parents, markers of low quality (E), and markers with at least one mismatch.The remaining 6504 markers were used to construct the linkage maps.In the first step, we implemented the BIN tool algorithm in the IciMapping software version 4.2 [51].The remaining markers were binned according to their segregation patterns.After binning, we grouped the markers using a logarithm of odds (LOD) threshold of 3.0 [52].Linkage groups were assigned according to the genomic position of the SNP markers determined during SNP calling.Recombination frequencies between markers were converted to centiMorgans (cM) using the Kosambi mapping function [53].We inspected the initial linkage map for duplicate lines, segregation distortion, switched alleles, and single and double crossovers (genotyping errors) using the R/qtl and R/ASMap packages available in the R Statistical Computing Environment v 2.15.1 [54,55].Finally, after removing low-quality markers and correcting for genotyping errors, the genotypic data of BIL1 and BIL2 with 2882 and 3404 high-quality markers, respectively, were used to construct the final QTL map in IciMapping software version 4.2.

QTL Analysis
The studied traits for BIL1 and BIL2 were used for QTL mapping in QTL IciMapping.We used the mean of the traits at each environment as well as their BLUPs.An inclusive composite interval mapping of QTL with additive and dominance effect (ICIM-ADD) analysis was performed.The significant LOD threshold (2.5) for declaring a QTL (α = 0.05) was determined from 10,000 permutations.We reported a LOD score of significant QTLs ranging from 2.5 to 5.2 in BIL1 and 2.5 to 18.6 in BIL2.

Conclusions
The systematic exploration, utilization, and evaluation of synthetic wheat developed using 43 Ae.tauschii accessions enabled the identification of several heat-tolerant lines and MTAs associated with important heat adaptation traits.In this study, we utilized Ae. tauschii donors of two heat-tolerant MSD lines to develop two BILs.The evaluation of the two BILs under field conditions with continuous heat stress revealed 16 heat-tolerant lines as well as 39 novel QTLs associated with heat stress.Specifically, the use of the two BILs permitted us to pinpoint desirable QTLs absent in elite wheat cultivars.The major and stable QTLs identified in this study, associated with the stress tolerance indices of GY and TKW, hold considerable promise and warrant further investigation to fully elucidate their potential applications.Efforts are underway to develop near-isogenic lines (NILs) using the selected heat-tolerant BIL lines to validate the identified QTLs and develop markers that will assist in wheat genetic improvement under heat stress.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/plants13030347/s1, Figure S1: Box plot illustrating the environmental effects on the nine evaluated traits in four environments: Dongola (DN), Waha (WA), Wad Medani season one (WM1), and Wad Medani season two (WM2).The boxes indicate the medians, interquartile range, and whiskers range.Tukey's honestly significant difference test was used to compare environments.Different letters indicated significant differences.Table S1.Overall mean, range, and broad sense heritability (h 2 ) for eight evaluated traits in the two-backcross recombinant inbred lines populations (BIL1 and BIL2) grown at Dongola (DN), Waha (WA), Wad Medani first season (WM1) and Wad Medani second season (WM2).Table S2.Correlation coefficients between the traits measured in both recombinant inbred lines populations BIL1 (above triangle) and BIL2 (down triangle) at Dongol (DN), Waha (WA), Wad Medani 1st season (WM1), and Wad Medani 2nd season (WM2).Table S3.Distribution of GRAS-Di markers used for map construction across A, B, and D-subgenomes in the backcrossed recombinant inbred lines (BIL2).Table S4.Marker alignment in the sub-genomes A, B, and D with their genetic and physical positions on each chromosome in BIL2 population.Table S5.QTLs detected by Inclusive composite interval mapping with additive and dominance effect (ICIM-ADD) analysis in both BIL1 and BIL2 under heat (H, WA, and WM1) and optimum (C, DN, and WM2) environments and the BLUP.

Figure 1 .
Figure 1.The daily maximum and minimum temperatures recorded during the growing season at (a) Dongola (DN), (b) Waha (WA), (c) Wad Medani first season (WM1), and (d) Wad Medani second season (WM2).The heading ranges (HRs, days) and grain filling duration (GFD, days) of both BIL1 and BIL2 are indicated by the dotted green and black vertical lines, respectively.

Figure 1 .
Figure 1.The daily maximum and minimum temperatures recorded during the growing season at (a) Dongola (DN), (b) Waha (WA), (c) Wad Medani first season (WM1), and (d) Wad Medani second season (WM2).The heading ranges (HRs, days) and grain filling duration (GFD, days) of both BIL1 and BIL2 are indicated by the dotted green and black vertical lines, respectively.

Figure 2 .
Figure 2. Comparison of grain yield at Dongola (DN) for the two populations (BIL1 and BIL2) using the stress tolerance index (STI).The STI-GY was calculated twice.First, STI1-GY was calculated regarding the GY values at DN as the non-stress environment and WA or WM1 as the hot environment.Second, STI2-GY was calculated considering the GY of WM2 as the non-stress environment and WM1 as the stress environment.The GY at DN was regressed against STI1-GY at (a) WA and (b) WM1.The GY at WM2 was regressed against STI2-GY at (c) WM1.The red circle and diamond represent the recurrent parent N61 in BIL1 and BIL2, respectively.The horizontal dashed red line is plotted based on the GY mean plus the Least Significant Difference (LSD 0.05) of both BIL1 and BIL2.Similarly, the horizontal dashed blue line corresponds to the GY mean plus the LSD of BIL1, and the dashed black line reflects the GY mean plus the LSD of BIL2.The vertical blue line is calculated based on the STI mean plus the LSD of BIL1, while the vertical black line is calculated based on the STI mean plus the LSD of BIL2.The horizontal blue arrows represent the heat tolerance ranges (HTR) for both BILs.

Figure 2 .
Figure 2. Comparison of grain yield at Dongola (DN) for the two populations (BIL1 and BIL2) using the stress tolerance index (STI).The STI-GY was calculated twice.First, STI1-GY was calculated regarding the GY values at DN as the non-stress environment and WA or WM1 as the hot environment.Second, STI2-GY was calculated considering the GY of WM2 as the non-stress environment and WM1 as the stress environment.The GY at DN was regressed against STI1-GY at (a) WA and (b) WM1.The GY at WM2 was regressed against STI2-GY at (c) WM1.The red circle and diamond represent the recurrent parent N61 in BIL1 and BIL2, respectively.The horizontal dashed red line is plotted based on the GY mean plus the Least Significant Difference (LSD 0.05 ) of both BIL1 and BIL2.Similarly, the horizontal dashed blue line corresponds to the GY mean plus the LSD of BIL1, and the dashed black line reflects the GY mean plus the LSD of BIL2.The vertical blue line is calculated based on the STI mean plus the LSD of BIL1, while the vertical black line is calculated based on the STI mean plus the LSD of BIL2.The horizontal blue arrows represent the heat tolerance ranges (HTR) for both BILs.

Figure 3 .
Figure 3. Genetic maps constructed using 3404 GRS-Di markers in BIL2.For BIL1, refer to Ahmed et al. [22].The positions of the QTLs are colored.

Figure 4 .
Figure 4.The detected QTLs in common regions in both backcrossed recombinant inbred line populations, BIL1 and BIL2, with co-localized previously reported QTLs.The QTLs detected in BIL1 are in black while those in BIL2 are in red.Normal fonts represent the optimum conditions, while bold fonts indicate the heat condition for both BIL1 and BIL2.The blue color represents the QTLs previously reported in our studies using the multiple synthetic derivative (MSD) line population.The green color represents QTLs reported by other researchers.

Figure 4 .
Figure 4.The detected QTLs in common regions in both backcrossed recombinant inbred line populations, BIL1 and BIL2, with co-localized previously reported QTLs.The QTLs detected in BIL1 are in black while those in BIL2 are in red.Normal fonts represent the optimum conditions, while bold fonts indicate the heat condition for both BIL1 and BIL2.The blue color represents the QTLs previously reported in our studies using the multiple synthetic derivative (MSD) line population.The green color represents QTLs reported by other researchers.

Table 1 .
The QTLs associated with stress tolerance indices (STIs) of grain yield (GY), biomass (BIO), and thousand-kernel weight (TKW) in two BILs grown in four environments.

Table 2 .
Key stable/major QTLs overlapped at least in two heat stress environments or the two BILs.