Heat Stress Tolerant-QTL Identified Using Backcrossed Recombinant Inbred Lines Derived From Intra-specifically Diverse Aegilops tauschii Accessions

1. 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 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, mostly a limited number and diversity of Ae. tauschii have 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 (MSD) 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 combination under field conditions in Sudan. Several MTAs were identified associated with heat, drought, heat-drought, and grain characteristics and quality [15,16,17,18,19,20].

The MSD as mixture population of several BC₁F₇ 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 the 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₁F₅ using the single-seed 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.

2. Results

2.1. Climate conditions during the growing seasons

Temperature data were recorded during the 2020/2021 growing season at Wad Medani (WM1) and during the 2021/2022 growing season 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 grain filling period, the maximum temperature ranged from 21.5-36.8°C at DN, 25.0-43.0°C at WA, 27.0-44.0°C at WM1 27.0-42.0°C at WM2. At DN, maximum temperatures ≥ 40.0°C during grain filling period were not reported, however, at WA, seven days were reported with maximum temperature ≥ 40.0°C. At WM1 and WM2, 19 and 11 days were reported with maximum temperatures ≥ 40.0°C, respectively. Likewise, number of days with mean temperatures ≥ 20.0°C during the grain filling period were 33, 57, 75, and 67 days at DN, WA, WM1, and WM2, respectively.

2.2. 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 by environment (G × E) interaction on PH in BIL2 (Table 1).

In BIL1, the mean grain yield (GY) was 5,553 kg ha^-1 in Dongola (DN), 2340 kg ha^-1 in Waha (WA), 2426 kg ha^-1 in Wad Medani first season (WM1), and 5,211 kg ha^-1 in Wad Medani second season (WM2). The mean GY in BIL2 was 5,408 kg ha^-1 in DN, 2,531 kg ha^-1 in WA, 2,451 kg ha^-1 in WM1, and 6,136 kg ha^-1 in WM2. The GY, BIO, days to heading (DH), days to maturity (DM), and PH were significantly reduced at WM1, and WA compared to DN and WM2 in both BILs (Supplementary Figure 1). The longest grain filling duration was recorded at WM1 followed by WM2, WA, and DN. The highest TKW was always recorded at WA in both BILs followed by DN, whereas the WM environments recorded the lowest kernel weight (Supplementary Figure 1)

High to moderate broad sense heritability (h²) estimates were found for GY, DH, and DM in both BILs (Table 1). Moderate to low h² estimates were found for GFD, BIO, and HI in both BILs. The h² estimated for PH and TKW varied between the two BILs (Table 1).

We calculated the stress tolerance index (STI) for GY, to identify heat-tolerant genotypes, and then performed regression analysis between GY at DN and STI1-GY at WA, and STI1-GY at WM1 (Figures 2a and b). We also performed regression analysis between the GY at WM2 and the STI2-GY at WM1 (Figure 2c). In BIL1, six genotypes (5.6%) exhibited higher STI values than the recurrent parent, N61, in both WA and WM1 (Figure 2a and b). Among the six genotypes, three and four lines showed higher GY than N61 at WA and WM1, respectively. In BIL2, 10 (6%) genotypes exhibited higher STI values than N61 at WA, out of which five lines showed higher GY than N61. At WM1, six (3.7%) genotypes exhibited higher STI values than N61, out of which three lines showed higher GY than N61 (Figures 2a and b). For STI2, 16 (15%) and seven (4.3%) genotypes in BIL1 and BIL2, respectively, exhibited higher STI2 values relative to N61 (Figure 2c). Only two genotypes in BIL1 and one in BIL2 showed higher GY than N61 (Figure 2c).

2.3. Relationship among traits

In both BILs, GY significantly and consistently correlated (P< 0.05) with BIO and HI in all environments (Supplementary Table 1). At DN, GY significantly correlated with seed number per spike (SN) (0.247, P<0.001) in BIL1 and with TKW (0.169, P <0.05) in BIL2. At WM1, GY significantly correlated with SN in BIL1 and BIL2 (0.437 and 0.353, respectively, P<0.001), whereas at WM2, GY correlated with TKW in BIL2 (0.244, P<0.001). The stress tolerance index for GY calculated taking DN as the non-stress environment (STI1-GY) at WA and WM1, significantly correlated with SN (0.318 and 0.395, respectively, P<0.001), BIO (0.573 and 0.702, respectively, P<0.001) and HI (0.424 and 0.405, respectively, P<0.001) in BIL1. In BIL2 at WA and WM1, the STI1-GY significantly correlated with BIO (0.687 and 0.688, respectively, P<0.001) and HI (0.472 and 0.348, respectively, P<0.001). Likewise, the stress tolerance index for GY calculated taking WM2 as non-stress environment (STI2-GY) at WM1 significantly correlated with SN (0.385, P<0.001), TKW (0.278, P<0.001), BIO (0.616, P<0.001) and HI (0.477, P<0.001) in BIL1, and with SN (0.278, P<0.001), BIO (0.691, P<0.001), and HI (0.365, P<0.001) in BIL2 (Supplementary Table 1).

2.4. Linkage maps for the BILs

Both BIL1 and BIL2 were genotyped by genotyping by random amplicon sequencing direct (GRAS-Di). The details of linkage map construction of BIL1 were mentioned in Ahmed et al. [22]. Briefly, the high-density linkage map was developed using 2,882 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 6,504 exhibiting polymorphisms between the synthetic wheat donor parent, Syn44, and the recurrent parent, N61. Of these polymorphic markers, 3,404 (52.3%) were of high quality, with an average of 162 markers per chromosome. A high-density linkage map was constructed utilizing the 3,404 markers spread across 21 linkage groups, covering a total genetic distance of 5,673.33 cM (Figure 3 and Supplementary Table 2). 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 (1,207, 35.46%) and D (1,201, 35.28%) subgenomes, which had total genetic lengths of 1,978.99 and 1,943.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 (Supplementary Table 2). 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 (Supplementary Table 2). 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 (Supplementary Table 3).

2.5. 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 were used. Using inclusive composite interval mapping of QTL 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 number of QTLs detected at WM1, DN, WA, and WM2 were 17, 10, 9, and 4, 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, and TKW. A total of 13 QTLs (32.5%) were identified in the A-subgenome for BIO, DH, GY, PH, STI1-BIO, STI1-TKW, STI2-TKW, and TKW, whereas 10 QTLs (25%) were identified in the B-subgenome for DM, GFD, HI, PH, SN, STI1-TKW, and TKW. The QTLs associated with GY were detected on chromosomes 6A, 1D, 3D, and 5D, whereas QTLs of both STI-GY were detected on chromosomes 1B, 1D, and 6D (Supplementary Table 4).

In BIL2, 153 QTLs were identified for the studied traits in the four environments. The identified QTLs were mapped on all chromosomes except 1B and 1D. The LOD score of the identified QTLs ranged from 2.50 to 18.63, and the phenotypic variation explained ranged from 0.32 to 18.91%. The highest number of QTLs was identified at WM1 (72), whereas 28, 27, and 26 QTLs were detected at WA, DN, and WM2, respectively. Fifty-six (36.60%) QTLs were identified in the A-subgenome associated with all studied traits, while 54 (35.29%) QTLs were identified in the D-subgenome for BIO, DH, DM, GFD, GY, HI, PH, STI2-BIO, STI2-GY, STI1-TKW, STI2-TKW, and TKW. The lowest number of QTLs (43, 28%) was identified in the B-subgenome associated with BIO, DH, DM, GFD, HI, PH, STI1-BIO, STI2-BIO, STI2-GY, and TKW. The QTLs of GY were detected on chromosomes 5A, 7A, and 5D, while those associated with both STI-GYs were detected on chromosomes 3A, 5A, 4B, and 6D (Supplementary Table 4).

2.5. 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 2, 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 10 QTLs, two on chromosome 1A and one on chromosome 2B were stable 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 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 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 2, Figure 4).

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

Chr1	Trait	Pop2	Pos3 (cM)	Left Marker	Right Marker	LOD4	PVE (%)5	Add6	Co-localized with
1A	STI1-TKW	BIL1	119	AMP0036610	AMP0034796	4.37	14.37	0.12	Guan et al., 2011
1A	STI1-TKW	BIL1	173	AMP0035547	AMP0004300	2.98	9.89	-0.10	Guan et al., 2011
1A	STI1-TKW	BIL1	118	AMP0036610	AMP0034796	3.74	11.37	0.11	Guan et al., 2011
1A	STI1-TKW	BIL1	173	AMP0035547	AMP0004300	4.27	14.02	-0.11	Guan et al., 2011
1A	STI2-TKW	BIL1	173	AMP0035547	AMP0004300	4.82	14.54	-0.10	Guan et al., 2011
1A	STI1-BIO	BIL1	148	AMP0034796	AMP0020845	2.64	5.14	-0.12
1B	STI1-GY	BIL1	144	AMP0017578	AMP0023142	2.60	7.49	0.07
1D	STI1-GY	BIL1	117	AMP0027815	AMP0029085	5.23	15.44	-0.10
1D	STI2-TKW	BIL1	236	AMP0005955	AMP0027742	3.52	10.14	-0.09
2B	STI1-TKW	BIL1	63	AMP0009891	AMP0006464	3.21	10.41	0.10
2B	STI1-TKW	BIL1	62	AMP0009891	AMP0006464	2.95	8.66	0.08	Paliwal et al., 2012
2B	STI1-BIO	BIL2	150	AMP0012513	AMP0026808	2.80	3.46	-0.17
2D	STI1-TKW	BIL2	21	AMP0020907	AMP0024533	4.11	3.15	0.19	Guan et al., 2018
3A	STI1-GY	BIL2	49	AMP0010424	AMP0003988	5.33	14.20	0.10
3A	STI1-TKW	BIL2	129	AMP0014988	AMP0016989	18.63	15.53	-0.19
3A	STI1-TKW	BIL2	137	AMP0029972	AMP0030211	10.85	8.05	0.14
3A	STI1-TKW	BIL2	178	AMP0007900	AMP0004728	4.15	2.82	-0.07
3A	STI2-TKW	BIL2	39	AMP0030786	AMP0010424	2.83	6.29	0.33
3A	STI1-BIO	BIL2	49	AMP0010424	AMP0003988	3.03	0.57	0.07
3D	STI2-BIO	BIL2	300	AMP0001446	AMP0012860	4.58	2.97	-0.31
4B	STI2-GY	BIL2	140	AMP0018665	AMP0020290	3.27	8.54	0.07
4B	STI1-BIO	BIL2	176	AMP0025189	AMP0026555	3.92	1.12	-0.28
4B	STI2-BIO	BIL2	178	AMP0026555	AMP0003848	5.10	2.89	-0.36
4D	STI1-TKW	BIL2	93	AMP0031292	AMP0028457	2.56	1.70	-0.08
4D	STI2-BIO	BIL2	105	AMP0009857	AMP0007548	4.21	2.01	-0.41
5A	STI1-GY	BIL2	88	AMP0011577	AMP0030240	3.81	10.50	0.08	Hassuni et al., 2019
5A	STI1-BIO	BIL2	3	AMP0003832	AMP0029058	2.75	6.20	0.20
5A	STI2-GY	BIL2	88	AMP0011577	AMP0030240	7.84	13.70	0.09	Hassuni et al., 2019
5A	STI2-GY	BIL2	168	AMP0008559	AMP0030185	3.26	5.65	0.06	Hassuni et al., 2019
5A	STI1-BIO	BIL2	3	AMP0003832	AMP0029058	2.94	3.09	0.19
5A	STI2-BIO	BIL2	87	AMP0025208	AMP0011577	4.14	1.12	0.08
5A	STI2-BIO	BIL2	227	AMP0001406	AMP0015434	5.74	2.43	-0.39
5D	STI1-TKW	BIL1	7	AMP0022256	AMP0000398	4.39	14.08	-0.12	Wang et al., 2021
5D	STI2-BIO	BIL2	204	AMP0010296	AMP0028613	2.95	3.10	-0.29
6D	STI1-GY	BIL1	53	AMP0036794	AMP0032738	3.03	8.67	0.07
6D	STI2-GY	BIL2	92	AMP0016445	AMP0014713	3.90	6.84	0.07
6D	STI2-TKW	BIL2	305	AMP0003394	AMP0027092	2.76	5.28	-0.10	Guan et al., 2018
7D	STI2-TKW	BIL1	89	AMP0019618	AMP0017004	2.96	8.30	-0.08	Paliwal et al., 2012
7D	STI2-BIO	BIL2	343	AMP0018976	AMP0002072	2.50	2.56	-0.36

¹ Chromosome, ² Populations, ³ Position, ⁴ Logarithm of odds, ⁵ Phenotypic variance explained, ⁶ Additive effect.

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

Chr1	Trait	Pop2	Pos3 (cM)	Left Marker	Right Marker	LOD4	PVE (%)5	Add6	Co-localized with
1A	STI1-TKW	BIL1	119	AMP0036610	AMP0034796	4.37	14.37	0.12	Guan et al., 2011
1A	STI1-TKW	BIL1	173	AMP0035547	AMP0004300	2.98	9.89	-0.10	Guan et al., 2011
1A	STI1-TKW	BIL1	118	AMP0036610	AMP0034796	3.74	11.37	0.11	Guan et al., 2011
1A	STI1-TKW	BIL1	173	AMP0035547	AMP0004300	4.27	14.02	-0.11	Guan et al., 2011
1A	STI2-TKW	BIL1	173	AMP0035547	AMP0004300	4.82	14.54	-0.10	Guan et al., 2011
1A	STI1-BIO	BIL1	148	AMP0034796	AMP0020845	2.64	5.14	-0.12
1B	STI1-GY	BIL1	144	AMP0017578	AMP0023142	2.60	7.49	0.07
1D	STI1-GY	BIL1	117	AMP0027815	AMP0029085	5.23	15.44	-0.10
1D	STI2-TKW	BIL1	236	AMP0005955	AMP0027742	3.52	10.14	-0.09
2B	STI1-TKW	BIL1	63	AMP0009891	AMP0006464	3.21	10.41	0.10
2B	STI1-TKW	BIL1	62	AMP0009891	AMP0006464	2.95	8.66	0.08	Paliwal et al., 2012
2B	STI1-BIO	BIL2	150	AMP0012513	AMP0026808	2.80	3.46	-0.17
2D	STI1-TKW	BIL2	21	AMP0020907	AMP0024533	4.11	3.15	0.19	Guan et al., 2018
3A	STI1-GY	BIL2	49	AMP0010424	AMP0003988	5.33	14.20	0.10
3A	STI1-TKW	BIL2	129	AMP0014988	AMP0016989	18.63	15.53	-0.19
3A	STI1-TKW	BIL2	137	AMP0029972	AMP0030211	10.85	8.05	0.14
3A	STI1-TKW	BIL2	178	AMP0007900	AMP0004728	4.15	2.82	-0.07
3A	STI2-TKW	BIL2	39	AMP0030786	AMP0010424	2.83	6.29	0.33
3A	STI1-BIO	BIL2	49	AMP0010424	AMP0003988	3.03	0.57	0.07
3D	STI2-BIO	BIL2	300	AMP0001446	AMP0012860	4.58	2.97	-0.31
4B	STI2-GY	BIL2	140	AMP0018665	AMP0020290	3.27	8.54	0.07
4B	STI1-BIO	BIL2	176	AMP0025189	AMP0026555	3.92	1.12	-0.28
4B	STI2-BIO	BIL2	178	AMP0026555	AMP0003848	5.10	2.89	-0.36
4D	STI1-TKW	BIL2	93	AMP0031292	AMP0028457	2.56	1.70	-0.08
4D	STI2-BIO	BIL2	105	AMP0009857	AMP0007548	4.21	2.01	-0.41
5A	STI1-GY	BIL2	88	AMP0011577	AMP0030240	3.81	10.50	0.08	Hassuni et al., 2019
5A	STI1-BIO	BIL2	3	AMP0003832	AMP0029058	2.75	6.20	0.20
5A	STI2-GY	BIL2	88	AMP0011577	AMP0030240	7.84	13.70	0.09	Hassuni et al., 2019
5A	STI2-GY	BIL2	168	AMP0008559	AMP0030185	3.26	5.65	0.06	Hassuni et al., 2019
5A	STI1-BIO	BIL2	3	AMP0003832	AMP0029058	2.94	3.09	0.19
5A	STI2-BIO	BIL2	87	AMP0025208	AMP0011577	4.14	1.12	0.08
5A	STI2-BIO	BIL2	227	AMP0001406	AMP0015434	5.74	2.43	-0.39
5D	STI1-TKW	BIL1	7	AMP0022256	AMP0000398	4.39	14.08	-0.12	Wang et al., 2021
5D	STI2-BIO	BIL2	204	AMP0010296	AMP0028613	2.95	3.10	-0.29
6D	STI1-GY	BIL1	53	AMP0036794	AMP0032738	3.03	8.67	0.07
6D	STI2-GY	BIL2	92	AMP0016445	AMP0014713	3.90	6.84	0.07
6D	STI2-TKW	BIL2	305	AMP0003394	AMP0027092	2.76	5.28	-0.10	Guan et al., 2018
7D	STI2-TKW	BIL1	89	AMP0019618	AMP0017004	2.96	8.30	-0.08	Paliwal et al., 2012
7D	STI2-BIO	BIL2	343	AMP0018976	AMP0002072	2.50	2.56	-0.36

¹ Chromosome, ² Populations, ³ Position, ⁴ Logarithm of odds, ⁵ Phenotypic variance explained, ⁶ Additive effect.

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. 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 2).

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 MSD (Figure 5). 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).

3. Discussion

In the face of climate change, evaluating wheat lines with limited genetic diversity or with similar genetic backgrounds from 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,21,23]. Therefore, a population with a wide diverse genetic background has been developed utilizing 43 Ae. tauschii accessions and named multiple synthetic derivatives (MSD) [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 climate-resilient traits [15,16,18,24,25]. Moreover, GWAS analysis identified several MTAs under heat, drought, and heat-drought conditions [15,16,19,20,25]. 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 MSD population and other bread wheat genotypes carrying exotic alleles [10,15,24]. 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 that 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 WM environments. This might be due to the significant G × E interactions observed for different traits, especially DH, DM, GFD, GY, and BIO (Table 1). 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 the stressed and non-stressed environments. In both BILs at the heat stress environment of WM1, 12 lines exhibited higher STI-GY, out of which seven lines had higher grain yield than N61. Similarly, at WA, 16 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 [26], as well as when exotic-derived lines from Ae. tauschii were compared under heat stress with elite wheat lines [10,27]. The high-yielding heat-tolerant lines identified in both BILs are promising and could be integrated into wheat heat stress tolerance breeding programs.

3.1. The high-resolution linkage maps

The GRAS-Di platform generates many chromosome-spanning genetic markers, allowing the construction of high-resolution linkage maps [28]. We constructed two genetic maps using 2,882 and 3,404 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 [28]. Nevertheless, the length of our maps was comparable to those reported in wheat [29,30]. The D-subgenome was longer in BIL1 than that of 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 [31,32] due to its low level of polymorphism [33]. 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.

3.2. QTL 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,34,35]. The D-subgenome contributed by 36.9% of the total QTLs identified including five out of nine GY QTLs and 13 out of 40 QTLs associated with stress tolerance index. Our results confirmed the importance of D-subgenome derived from Ae. tauschii for mining stress adaptive traits [18,25].

As shown in Figure 4, a number of the QTLs identified in this study were co-located with some of the previously reported QTLs under heat stress conditions on chromosomes 1A, 2A, 2B, 5B, 7B, 1D, and 2D [36,37,38,39,40,41,42,43]. For instance, a number of STI-TKW QTLs were co-localized with previously identified QTLs for HSI-TGW (heat susceptibility index of thousand-grain weight) and TKW under heat and combined heat-drought on chromosomes 1A, 2B, and 2D.

3.2.1. Identification of stable major QTL 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, the 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.

We used STI for GY, BIO, and TKW, and identified new QTLs for heat stress tolerance, including STI-GY on chromosomes 1B, 4B, 1D, and 6D; STI-TKW on chromosomes 3A, 2B, and 4D; and STI-BIO on chromosomes 1A, 3A, 3D, 4B, 4D, 5A, and 7D in both BILs. Recently, Farhad et al. [45] reviewed QTLs identified under heat stress in wheat and reported several QTLs related to heat stress tolerance in wheat identified so far. Several QTLs for heat and stress susceptibility indices (HSI, SSI) for kernel weight, grain filling duration, and grain yield have been identified. For instance, QTLs for SSI-GY on chromosomes 2A, 5A, and 5B [45,46], and HSI-TKW on chromosomes 1A, 2B, 2D, 5B, 5D, 6A, 6B, 6D, and 7D [37,39,47] have been reported. These newly identified heat-stress QTLs have the potential to improve GY under heat stress after further validation.

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 12 QTLs, eight were for STI-TKW on chromosomes 1A, 3A, 2B, 1D, and 5D, while four QTLs were found for STI-GY on chromosomes 3A, 5A, and 1D.

3.2.2. 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 [43]; however, 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 BILs 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. Up to our knowledge, no heat stress related QTLs have been reported in these regions of the chromosomes. Previously, utilizing MSD population several MTAs have been identified in these regions [15,16,18], 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 the 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.

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 BILs 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.

4. Materials and Methods

4.1. Plant Materials

In this study, we used two populations of backcrossed recombinant inbred lines (BILs, BC₁F₅). Through intensive evaluation of multiple synthetic derivatives (MSD), developed utilizing 43 accessions of Ae. tauschii and the durum wheat cultivar ‘Langdon’, then crossed and backcrossed to the hexaploid wheat cultivar N61 [8,14], under field heat stress conditions, two heat tolerant MSD lines were identified [24]. We utilized the two Ae. tauschii parents (KU-2039 and KU-2124) of the two synthetic derivatives (Syn32 and Syn44) for the development of the two BILs. KU-2039 was collected from Afghanistan, whereas KU-2124 was collected from Iran. The two accessions were classified into two distantly related lineages, TauL1 for KU-2039 and Taul2 for KU-2124.

First, the two Ae. tauschii accessions (DD genome) were crossed with durum wheat cultivar ‘Langdon‘ (AABB genome)’ to develop two primaries synthetic hexaploid wheat (SHW). Then the two SHW were crossed and backcrossed to N61 to generate BC₁F₁. The single seed descent (SSD) method was used to develop the two BILs (BC₁F₅). Initially, the first BIL, generated from KU-2039 (BIL1) consisted of 166 lines, whereas the second BIL, generated from KU-2124 (BIL2), consisted of 236 lines. All lines with very late headings or those that failed to flower due to the vernalization requirement were excluded after evaluation under the heat stress conditions at Wad Medani, Sudan. The final numbers used for the multi-environment evaluation in BIL1 and BIL2 were 107 and 164 lines, respectively.

4.2. Experimental Sites and Design

We conducted field experiments in two seasons (2020/21 and 2021/22) in Sudan. In the first season, the experiment was conducted at the Gezira Research Station Farm (GRSF), Agricultural Research Corporation (ARC), Wad Medani (14°24′N, 29°33′E, 407 m a.s.l. abbreviated as WM1). In the second season, in addition to GRSF (abbreviated as WM2), the experiment was conducted at Dongola Research Station Farm, ARC, Dongola, Northern Sudan (19°08′N, 30°27′E, 239 m a.s.l.; abbreviated as DN) and Waha Farm, south of Khartoum (32° 33’ 35.6328’’ E, 380 m a.s.l.; abbreviated as WA). For more details on the characteristics of the experimental sites at DN and GRSF, refer to Elbashir et al. [23] and Mohamed et al. [19]. The Waha site is characterized as a semi-arid area with mostly haplusterts soil type, clay loam in texture, with pH range of 7.5-8.0, low N (0.05-0.20%), and low available P (30-60 kg ha^-1) (http://susis.sd/south-region-0-30/, accessed on November 12, 2023). All experiments were arranged in an augmented randomized complete block design with five blocks and six checks (five Sudanese cultivars and the recurrent parent, N61.

At all sites, seeds of each genotype were manually sown in a plot consisting of two rows, 1.0 m long and 0.2 m apart. At WM1 and WM2, sowing was done during the 4^th week of November, whereas at DN and WA, sowing was done during the 2^nd week of December. All crop management including seed treatments, irrigation, fertilizer application, weed control, etc. followed the recommendations of the ARC for wheat production in Sudan as previously described [16,19,24]. During the winter season, there is usually no rainfall at any of the experimental sites. Therefore, wheat is fully irrigated with irrigation intervals of 10-14 days, depending on weather conditions, mainly temperature, to avoid water stress.

4.3. Phenotyping of BIL populations

4.3.1. Traits evaluation

At all sites, phenological, yield and some yield component traits were recorded. Days to heading (DH), was recorded as the number of days from first irrigation until 50% of the spikes merged from the leaf sheath. Days to maturity (DM) was the number of days from first irrigation until 90% of the plants had lost the green color of the glumes. Grain filling duration (GFD) was calculated as the interval in days between DH and DM. Plant height (PH, cm) was taken at maturity by measuring the distance from the ground to the top of the spike, excluding awns.

Grain yield (GY) and related traits, viz biomass (BIO), thousand kernel weight (TKW), harvest index (HI), and seed number/spike (SN) were recorded. GY and BIO were determined as grain weight and above-ground dry weight per plot, respectively, and then 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 et al. [49] as

$$STI-GY=\frac{\left({\mathrm{Y}}_N\right)\left({\mathrm{Y}}_S\right)}{{\left(Y_{\overline{N}}\right)}^2}$$

where $Y_N$ is GY under non-stress condition, $Y_S$ is GY under heat stress conditions, and $Y_{\overline{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.

4.3.2. Statistical analysis of phenotypic data

Analysis of variance (ANOVA) for all traits was performed in Plant Breeding Tools v. 1.4.2 [50]. The genotype (G), environment (E), and (G x E) interaction were considered as random effects. We used Tukey’s honestly significant difference (HSD) test for between-environment comparisons and the T-test for between-BILs comparisons. Pearson’s correlation coefficient among traits in each environment was calculated using IBM SPSS Statistics for Windows v. 28.0.1.1(15) [51]. Broad sense heritability was estimated in Plant Breeding Tools.

4.4. Genotyping of the BILs, map construction, and QTL analysis

4.4.1. DNA extraction

Total genomic DNA was extracted from 2-week-old leaves of BILs using the CTAB method [52]. The DNA samples (20 µL; 50-100 ng µL^-1) were sent to Eurofins Genomics Company, Japan (https://eurofinsgenomics.jp/jp/home/) for whole genome scan with genotyping by random amplicon sequencing direct (GRAS-Di) markers.

4.4.2. Maps construction

The linkage map for both BILs was constructed using GRAS-Di markers. The details of BIL1 map construction have been mentioned in Ahmed et al. [22]. The same procedure and steps were also followed with BIL2. In brief, 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 6,504 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 [53]. 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 [54]. 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 [55]. 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 [56]. Finally, after removing low-quality markers and correcting for genotyping errors, the genotypic data of BIL1 and BIL2 with 2,882 and 3,404 high-quality markers, respectively, were used to construct the final QTL map in IciMapping software version 4.2.

4.4.3. QTL analysis

The studied traits for BIL1 and BIL2 were used for QTL mapping in QTL IciMapping. 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 QTL ranging from 2.5 to 5.2 in BIL1 and 2.5 to 18.6 in BIL2.

5. 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. 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.