Identification of Single-Nucleotide Polymorphisms (SNPs) Associated with Heat Tolerance at the Reproductive Stage in Synthetic Hexaploid Wheats Using GWAS

The projected rise in global ambient temperature by 3–5 °C by the end of this century, along with unpredicted heat waves during critical crop growth stages, can drastically reduce grain yield and will pose a great food security challenge. It is therefore important to identify wheat genetic resources able to withstand high temperatures, discover genes underpinning resilience to higher temperatures, and deploy such genetic resources in wheat breeding to develop heat-tolerant cultivars. In this study, 180 accessions of synthetic hexaploid wheats (SHWs) were evaluated under normal and late wheat growing seasons (to expose them to higher temperatures) at three locations (Islamabad, Bahawalpur, and Tando Jam), and data were collected on 11 morphological and yield-related traits. The diversity panel was genotyped with a 50 K SNP array to conduct genome-wide association studies (GWASs) for heat tolerance in SHW. A known heat-tolerance locus, TaHST1, was profiled to identify different haplotypes of this locus in SHWs and their association with grain yield and related traits in SHWs. There was a 36% decrease in grain yield (GY), a 23% decrease in thousand-grain weight (TKW), and an 18% decrease in grains per spike (GpS) across three locations in the population due to the heat stress conditions. GWASs identified 143 quantitative trait nucleotides (QTNs) distributed over all 21 chromosomes in the SHWs. Out of these, 52 QTNs were associated with morphological and yield-related traits under heat stress, while 15 of them were pleiotropically associated with multiple traits. The heat shock protein (HSP) framework of the wheat genome was then aligned with the QTNs identified in this study. Seventeen QTNs were in proximity to HSPs on chr2B, chr3D, chr5A, chr5B, chr6D, and chr7D. It is likely that QTNs on the D genome and those in proximity to HSPs may carry novel alleles for heat-tolerance genes. The analysis of TaHST1 indicated that 15 haplotypes were present in the SHWs for this locus, while hap1 showed the highest frequency of 25% (33 SHWs). These haplotypes were significantly associated with yield-related traits in the SHWs. New alleles associated with yield-related traits in SHWs could be an excellent reservoir for breeding deployment.


Introduction
Wheat is a major staple food for almost 40% of the global population, cultivated on almost 220 million hectares of the land, with 750 million tons of production and a trade associated with HST at the reproductive stage in synthetic hexaploid wheat accessions, and (iii) characterize the TaHST1 locus and identify haplotype variation associated with heat stress at the reproductive stage.

Phenotypic Data Analysis
Phenotypic data showed significant variation among all traits in both conditions (control and HS). The analysis of variance (ANOVA) for all of the morphological traits is presented in Table S2. All of the traits showed significant variations for genotypes, genotype × location, genotype × year, genotype × treatment, and other possible combinations, with the exception of genotype × year interaction for some traits (e.g., PH and TKW).
The ambient environmental variables-such as daily maximum temperature (T max ), daily minimum temperature (T min ), daily average temperature (T avg ) in • C, and relative humidity (%) during the field trials-are described in Table S3. For example, among all locations, Sindh received the hottest days (>30 • C) compared to Bahawalpur and Islamabad. The descriptive data presented in Table 1 show the variations in SHWs under control and HS conditions. In control conditions, the mean DH was 119 days, with a range from 102 days (AUS33402) to 157 days (AUS34259), while in HS conditions the mean DH was 89 days, with a range from 77 days (AUS33403, AUS33415, AUS33421, AUS34453, and AUS34458) to 115 days (AUS34257). In control conditions, the mean value for GY was 9.4 t/ha, with a range from 7 t/ha to 12 t/ha, while in HS conditions the mean value for GY was 6 t/ha, with a range from 4 t/ha (AUS30672, AUS33386, AUS33397, AUS33404, AUS34235, and AUS34257) to 8 t/ha (AUS30627, AUS30642, AUS30648, AUS33414, and AUS34231). For TKW in control conditions, the mean value was 47 g, with a range from 19 g (AUS30660) to 68 g (AUS33417), while in HS conditions the mean value was 36 g, with a range from 18 g (AUS34257) to 49 g (AUS30296). Phenotypic data are presented in Supplementary Table S2. Pearson's correlation coefficients are given in Figure 1 for all traits in control and HS conditions, along with their distribution graphs. Positive correlation was observed between yield-related traits such as DH and TKW (r = 0.566 ***), GY and DH (r = 0.718 ***), SW and SL (r = 0.518 ***), PH and TKW (r = 0.496 ***), SW and GY (r = 0.590 ***), GY and TKW (r = 0.510 ***), GY and SL (r = 0.784 ***), and SL (r = 0.752) and DM and SL (r = 0.785 ***). TKW (r = 0.510 ***), GY and SL (r = 0.784 ***), and SL (r = 0.752) and DM and SL (r = 0.785 ***).

Marker Statistics
In total, 66,836 markers were genotyped using the 50 K SNP chip. Thirty-eight markers were removed after quality control due to missing data. A further 27,485 markers were removed due to having MAF < 0.05. Subsequently, 39,313 markers were used for the GWASs. These SNPs were distributed on all chromosomes, with the maximum number of SNPs on chr6A (n = 2484) and the minimum number of SNPs on chr4D (n = 794). The B-genome contained the greatest number of SNPs (n = 15,190), followed by the A-genome (n = 14,269) and the D-genome (n = 9926).

QTNs Associated with Yield Traits under Heat Stress
In total, 143 QTNs were identified, out of which 48 QTNs were identified in control conditions, 52 QTNs were identified in HS, and 43 QTNs were identified for HSI. The greatest number of QTNs was identified for GpS (n = 24), while the fewest QTNs were identified for SpSP (n = 6). On chr2B and 5B, the highest number of QTNs (n = 12) was identified, while the lowest number of QTNs (n = 1) was found on chr5A and chr7A. In total, two QTNs were consistently identified in control and HS conditions. A QTN on chr7B at 437 Mb was consistently associated with DM in control and HS conditions. Another QTN on chr6A at 390 Mb was consistently associated with DM both in control and HS conditions.
In total, 15 QTNs were pleiotropically associated with multiple traits. A QTN on chr1A at~560 Mb was associated with multiple yield-related traits, including GWpS, TpP, and SpSP, in both control and HS conditions. A QTN on chr2A at~44 Mbs was associated with GpS and TKW in both HS and control conditions. On chr2B, a QTN found at~6 Mb was associated with GWpS and DM in both HSI and control conditions. Another QTN on the same chromosome at~731 Mb was associated with TpP and PH under HS conditions. A QTN on chr2D at~598 Mb was associated with GpS and TKW in both HS and control conditions; another QTN on the same chromosome at~606 Mb was found to be associated with SL and SW in control and HSI conditions. A QTN was found on chr3A at~727 Mb, associated with SL and TKW in both HSI and control conditions. On chr3B, three QTNs were found to be associated with multiple yield traits: one at~625 Mb, associated with TKW and DM in control conditions; another at~637 Mb was associated with GpS and TKW in both HS and HSI conditions; and a QTN at~822 Mb was associated with SL and GpS in both HS and control conditions. A QTN was found on chr3D at~602 MB, linked with DH and GY in HS conditions. On chr5B, a QTN was found at~588 Mb, associated with GY and GpS in both HS and HSI conditions. On chr6A, two QTNs were found: one at~270 Mb, associated with DH and DM in both HS and HSI conditions; and the other at~380 Mb, associated with DH, TpP, and DM in both HS and HSI conditions. A QTN on chr6D at 458 Mb was linked with DH, PH, TKW, and GY in both HS and HSI conditions. On chr7B, a QTN at~437 Mb was found to be associated with DH and DM in both HS and control conditions. All of the QTNs associated with multiple traits in HS conditions are presented in Figure 2 and Table 2.    Likewise, the QTNs associated with multiple traits were designated as pleiotropic QTNs (Table 3). In total, 14 pleiotropic QTNs were identified, out of which 3 QTNs on chr2D, chr3D, and chr6D were contributed by the D-sub genome and may carry novel alleles from Ae. tauschii. A QTN on chr2A at ~44 Mb was associated with GpS and TKW under both HS and control conditions. Another QTN on chr3D was associated with DH and GY and could be an important locus because plants tend to flower early to escape heat stress, and simultaneous control of such traits is an important finding. Likewise, the QTNs associated with multiple traits were designated as pleiotropic QTNs (Table 3). In total, 14 pleiotropic QTNs were identified, out of which 3 QTNs on chr2D, chr3D, and chr6D were contributed by the D-sub genome and may carry novel alleles from Ae. tauschii. A QTN on chr2A at~44 Mb was associated with GpS and TKW under both HS and control conditions. Another QTN on chr3D was associated with DH and GY and could be an important locus because plants tend to flower early to escape heat stress, and simultaneous control of such traits is an important finding.

Co-Localization of QTNs and Heat Shock Proteins
In total, eight QTNs were found in proximity to HSPs under HS conditions (Table 4). Two QTNs on chr2B were found in proximity to HSPs, but the QTN (AX-94760904) associated with TpP was the closest to TaHSP70.24, at~731 Mb. A QTN (AX-95237322) on chr3D at~45 Mb associated with PH was found in proximity to HSP TaHSP40.95. A QTN (AX-109911103) on chr5A at~611 Mb and a QTN (AX-94886530) on chr5B at~186 Mb, associated with GWpS, were found in proximity to HSPs TaHSP40.150 and TaHSP40.155, respectively. Based on the SNP effect, favorable and unfavorable alleles were identified, and their frequencies were determined. For DH, seven favorable and unfavorable alleles were identified in HS conditions. The scatterplot in Figure 3 shows the effects of the favorable and unfavorable alleles on DH; with the increase in the number of favorable alleles, the number of days to heading decreased significantly, while an increase in the number of unfavorable alleles reduced the DH. A similar pattern was observed in GY: with the increase in the number of favorable alleles, GY increased, while with the increase in the number of unfavorable alleles, GY decreased. The coefficients of determination (R 2 ) indicated that effect of favorable alleles ranged from R 2 = 0.77 (DH) to R 2 = 0.66 (GY) (Figure 3), while for unfavorable alleles the range of effect was from R 2 = 0.75 (DH) to R 2 = 0.67 (GY) (Figure 3).

Characterization of the TaHST1 Locus in the Diversity Panel
Characterization of the TaHST1 locus in synthetic wheat accessions was performed using PCR-based primers. Based on the results of these primers, 15 haplotypes were identified. Hap-1 showed the highest frequency, at 25%, and all markers showed positive amplification. Hap-2 was present in 16% of accessions, and only two markers showed positive amplification, while Hap-3 and Hap-4 were present in 15% and 10% of accessions, with one and two deleted sites, respectively. The effect of top Hap1 to Hap4 was significant on HSI for TpP ( Figure 4). All of the haplotypes, along with their frequency and phenotypic data under HS, are presented in Table 5.

PH
AX   identified. Hap-1 showed the highest frequency, at 25%, and all markers showed positive amplification. Hap-2 was present in 16% of accessions, and only two markers showed positive amplification, while Hap-3 and Hap-4 were present in 15% and 10% of accessions, with one and two deleted sites, respectively. The effect of top Hap1 to Hap4 was significant on HSI for TpP ( Figure 4). All of the haplotypes, along with their frequency and phenotypic data under HS, are presented in Table 5.

Phenotypic Variation for Agronomic Traits in SHWs under Heat Stress
In the past, a small number of SHWs and their derivatives were characterized to explain why they are superior to the conventional bread wheat genotypes, and they were evaluated for high-temperature tolerance [22][23][24][25]. Here, we describe the performance of a larger collection of 180 SHWs, which were sown at the normal time and late to characterize their variation in many agronomic traits and their relationships under heat stress conditions. ANOVA illustrated significant variation of accessions under normal and heat stress conditions, and the outcomes were consistent with those previously documented for bread wheat when identifying the genotype by environmental interactions and considerable variance in agronomic traits [26,27]. A greater variation was shown by SHWs for all of those traits that could be exploited in breeding for heat tolerance [25]. Overall, the highest broad-sense heritability was observed for DH, indicating that this trait has a high response to selection, and under HS conditions the large variations observed in SHWs could be very useful in heat tolerance breeding programs for wheat germplasm to fight high-temperature stress. Early heading and prematurity are preferable for high yield gains under heat-stressed conditions, due to the stress avoidance mechanism.
As expected, the GY of SHWs was reduced significantly under HS conditions. Conversely, late heading and late maturity resulted in low GY, demonstrating that the enhanced production is a result of their adaptation and capability to avoid late heat stress. To avoid the late high-temperature stress by the end of the season, early heading and maturity enable the SHWs to fill their grains casually. PH showed high-to-moderate heritability estimates and could also be a useful trait in HS tolerance breeding. To bring improvement in wheat, GY has always been the most important criterion for direct breeding, and the further improvement of GY through direct selection may not be easy, due to the low heritability estimates suggested by the findings of the present study. The low heritability estimates obtained for GY suggest that direct selection of the GY for further betterment may be difficult in wheat. Therefore, indirect selection for increased GY at elevated temperatures using yield-related traits is an appropriate criterion for screening for tolerance to abiotic stresses, including heat and drought.

Effects of Heat Stress in Synthetic Hexaploid Wheat
It was observed that under HS conditions, all of the traits showed a decreasing trend, except for SpSP. This is likely due to the fact that spikes emerge at the onset of the reproductive stage when temperature conditions are optimal. Some SHW accessions (e.g., AUS30284, AUS30286, AUS30297, AUS30637, AUS30645, and AUS30660) showed higher resilience to HS conditions for yield-related traits.
GY is highly affected by the environment, as described by Sharma et al. [28]. During heat stress conditions, a noticeable decrease in yield and other developmental traits was observed, and many other reports of the heat stress responses of wheat also drew the same conclusion [28]. High-temperature stress has shown significant decreases in PH [29], maturity time [30], flowering time [31], number of grains per spike and weight of grains per spike [32], TKW [33,34], and GY [35]. The reduction in GY was observed more during the year 2015 as compared to 2014 (39% and 36%, respectively), because of higher temperatures. During this experiment, we observed that the 2 • C rise in temperature in the year 2015 reduced the GY by 5% compared to 2014. A previous study observed a 3 to 4% decrease in wheat production under controlled conditions for every 1 • C increase above 15 • C, and the number of kernels declined by 12.5% for every 1 • C increase from 25/20 • C to 35/20 • C [33]. Similarly, it has been indicated by global simulations that the production of wheat will decrease by 6% on average for every 1 • C rise in temperature, equating to a yield decline of almost 42 million tons [36].

GWASs for the Identification of Loci Associated with Heat Stress Tolerance
Identification of QTNs associated with HST is crucial for breeders. Therefore, a number of gene mapping studies were conducted at the flowering and reproductive stages [13,37]. Recently, Tanin et al. [21] compiled all of the abiotic-stress-related QTLs in wheat as meta-QTLs and provided a detailed framework for comparative studies. In total, 11 meta-QTLs were aligned with those found in our work, including 3QTLs for GY on chr1B, chr1D, and chr2B. Two QTLs were identified by Sangwan et al. [38] on chr2A and 2D for DH. In our study, no QTNs were identified on chr2A or2D. However, eight QTNs were identified on chr1A, chr3D, chr5B, chr6A, and chr7B in HS. Similarly, two QTLs were identified on chr2D for DM [38]. In contrast, four QTNs were identified in our study on chr3D, chr6A, and chr7B. For PH, four QTNs were identified on chr2B, chr3A, chr3D, and chr6D. Meanwhile, in previous studies, two QTLs were identified on chr2A and chr4A for PH [38]. For SL, three QTLs were identified by Maulana et al. [15] on chr3B and chr7D.
In our studies, a QTN was identified on chr3B for SL in HS. Three QTNs were identified for SpSP on chr1A, chr2B, and chr2D. No QTL had been identified previously for SpSP on these chromosomes. In contrast, eight QTLs were identified on chr1D, chr3A, chr4B, chr4D, chr5A, chr6A, and chr6B [39]. Three QTLs were identified on chr2B, chr3A, and chr5B for GpS [40]. In our study, three QTNs were identified on same chromosomes for GpS in HS. For GY, two QTLs were identified on chr2B and chr2D [41]. In our studies, five QTNs were identified on chr1B, chr1D, chr3D, chr6D, and chr7D for GY in HS. For TKW, in our study, only one QTN was identified on chr6D in HS, while five QTLs were identified by Li et al. [40]-on chr1D, chr4D, chr5A, chr5D, and chr6B. As discussed by Ogbonnaya et al. [42], favorable alleles in a single environment show little improvement in multiple-stress environments. We also observed a similar phenomenon in the present study, where environment-specific MTAs and MTAs for the average environments did not significantly overlap. Differences in loci controlling stability and environment-specificity suggest that there may be separate evolutionary trajectories for them.

Co-Localization of QTNs and HeatShock Proteins
All of the HS-specific QTNs in proximity to HSPs could be ideal candidates for genes underpinning heat tolerance, and due to the use of SHWs they are likely to carry new alleles of HSPs. Several HSPs (e.g., TaHSP100, TaHSP70, and TaHSP90) are induced underheat stress conditions [26]. In this context, the QTN for SW on chr4B is likely to be a good candidate for further studies, because its closet HSP is TaHSP70.52, which was reported to be induced during heat stress. There is a need for detailed analysis of the HSPs reported in this study for further validation.

Characterization of the TaHST1 Locus in the Diversity Panel
The TaSHT1 locus plays a significant role in HST at both the seedling and reproductive stages. Five markers were used to detect the presence of theTaSHT1 locus in SHW [17]. Recently, Khan et al. (2022) [43] identified 24 haplotypes of TaHST1 in Pakistani wheat cultivars. Fifteen haplotypes were determined based on the results of these markers. The results indicated that the presence of this QTL has significant effects on wheat under heat stress conditions. Two markers for this locus (Xhau-2 and -3) have shown significant associations with phenotypic traits under heat stress conditions. It is likely that haplotype distribution in SHW is different from that in cultivated wheats, and that new alleles are available due to the use of durum wheat in developing SHWs. The combination of 127 bp alleles at Xhau-4 and the presence of alleles at Xhau-5 have a positive impact on yieldrelated traits and could be used as a positive marker for selecting heat-tolerant accessions.
Conclusively, the data generated here on a collection of SHWs represent a promising source for selecting heat-tolerant candidates for use in breeding, and the loci identified here can enhance our knowledge of the distribution of heat-tolerant alleles in wheat.

Plant Material
A collection of 180 SHWs was used for this experiment (Table S1). These SHWs were derived from the combinations of 44 durum wheat cultivars and 100 Ae. Tauschii accessions as a subset of the panels previously characterized for grain morphology [44], biotic stress resistances [45,46], and grain quality [47].

Field Evaluation
To allow exposure of the genotypes to different heat treatments for the duration of flowering and ripening, six field tests were carried outat three different sites for 2 years (2014 and 2015), with two sowing time treatments (normal planting, termed as control; and late sowing, termed as heat stress (HS)). The three locations were the Regional Agricultural Research Institute at Bahawalpur in Punjab Province (BWP), the National Agricultural Research Centre (NARC) at Islamabad (ISB), and the Nuclear Institute of Agriculture at Tando Jam in Sindh (SND). The dates of normal sowing and late sowing were 15 November and 20 December, respectively. The sowing of genotypes was performed in two-row plots that were 2 m long and spaced 30 cm apart, with 15cm between rows, and in a randomized complete block design with three replications. The weight of seeds sown per plot was adjusted to obtain a constant seeding rate of 30 viable seeds per row by considering measurements of TKW and germination level. A small-plot grain drill was used for sowing purposes. Before plantation, triple superphosphate (4.3 g m −2 of P 2 O 5 ) was applied by furrow placement, and urea was applied earlier than the second irrigation (8.6 g m −2 of N). During the crop season, irrigation was carried out four times at the following growth stages: (1) from germination to the appearance of seedlings, (2) from emergence to the double-ridge stage, (3) from double-ridge to anthesis, and (4) from anthesis to maturity [48]. The following traits were studied: plant height (PH, cm), measured as the distance from the soil level to the top of the spike, excluding awns; days to heading (DH); days to maturity (DM), considered as the number of days from sowinguntil the plant reached 50% maturity; grain number per square meter (GN), counted from 20 spikes randomly from every plot; tillers per plant (TpP), which were randomly selected from each genotype; spike length (SL) from the origin of the spike to the tip of the spike, without awns; spikelets per spike (SpSP), which were collected from each accession by picking three spikes and then counting the spikelets; weight of spikes per spike (SW), which was determined by weighing 3 spikes from each line; grain weight per spike (GWpS), which was calculated by measuring the weight of three randomly selected spikes from each row; thousand-kernel weight (TKW, g), which was measured after harvest by weighing double samples of 500 kernels from each plot; and grain yield (GY, kg m −2 ), which was calculated as the weight of grains produced per unit area.
The heat susceptibility index (HSI) was calculated using the following formula: for example, for grain yield, HSIGY = (1 − Y/Yp)/D, where Y represents the yield of a genotype during the late planting, Yp shows the mean yield of genotypes at the time of normal planting, and stress intensity D = 1 − X/Xp, where X indicates the mean Y of all genotypes and Xp is the mean Yp of all genotypes.

DNA Extraction and Genotyping
Extraction of DNA was carried out using an amended CTAB method [49]. Synthetic wheat accessions were genotyped by using a 50 K SNP chip containing 66,798 SNPs. After filtering with a minor allele frequency of <5% and missing data of <20%, a total of 39,383 SNPs were then used for GWASs.

Analysis of Phenotypic Data
The phenotypic data were collected and arranged in Microsoft Excel. To test the statistical significance among different sources of variation, analysis of variance (ANOVA) was applied for all 11 traits. Phenotypic effects in the ANOVA model were divided into overall mean, genotypic effect, effect of random error, effect on replicate (i.e., block) in the environment (including both year and location), effect of environment, genotype with effect of treatment, genotype with environmental effect, and effect of treatment. Let y lijk be the observed value of the trait of interest, where i is the accession and kis the replicate under the jth environmental condition along with the years and locations in this study, while l is the treatment. The linear model applied in ANOVA is as follows: where l = 1, 2, . . . , L (L = 2 for normal and heat stress treatments), i = 1, 2, . . . , n (n = 203), j = 1, 2, . . . , e (e = 6 with 3 locations and 2 years), k = 1, 2, . . . , r (r = 2) is the total mean value of the entire population, R k/j is the kth effect on the replicate by the jth environment, G i is the effect of the genotype on the ith accession, E j is the effect of the environment on the jth environment, GE ij is the effect of the interaction between the ith accession and the jth environment, GD il is the effect of the interaction between the ith accession and the lth treatment, ε lijk is the random error effect (which is supposed to be distributed normally, with a mean value of zero), and the variance is σ 2 ε . The ANOVA illustrated above was utilized in SAS software (SAS Institute, Cary, NC, USA, 2007) together with the GLM methods.
The genotypic values from the best linear unbiased prediction (BLUP) under normal and heat stress treatments for all of the accessions were considered as the subsequent phenotypes for comparison. The BLUP values were deliberated as follows: the trait experimental value was defined as y ijk , where i is the accession, k is the replicate, and j is the environment, with the locations and years used in this study. The BLUP values for the mixed model were determined as follows: where i = 1, 2, . . . , n (n = 203), j = 1, 2, . . . , e (e = 6 by three locations and two years), k = 1, 2, . . . , and r (r = 2), E j ,GE ij, G i , µ, and R k/j , are as described above. Otherwise, normal distribution was followed, and all of the effects were observed as random effects GE are the variances representing the replicate, genotype, environment, and genotypeenvironment interaction, respectively. The values of BLUP were measured with the help of the MIXED procedure in the SAS software (SAS Institute, Cary, NC, USA, 2007). Pearson's correlation was computed using the R program.

Genotypic Data Analysis and GWASs
GWASs for 11 phenotypic traits were conducted using a recently developed model selection algorithm-the Fixed and Random Model Circulating Probability Unification (FarmCPU; Liu et al., 2016) [50]. FarmCPU takes into account the confounding problem between covariates and test markers using a fixed-effects model (FEM) and a randomeffects model (REM). The first five principal components calculated using TASSEL was used as covariates. The default p-value threshold that FarmCPU used was a Bonferroni-corrected threshold of 0.01. This Bonferroni-corrected threshold was overly restrictive when the LD among genotypic markers was large, so the threshold was calculated using the formula "p = 0.05/number of markers", with 1000 permutations. In this function, the phenotypes were permuted to break the relationship with the genotypes. A vector of minimum pvalues of each experiment was outputted, and the 95% quantile value of the vector was recommended for the p threshold in the FarmCPU model. The quantile-quantile (Q-Q) plot was used for assessing the fitness of the model to the population structure.
All of the SNPs associated with phenotypes under heat stress conditions were compared with the heat shock protein (HSP) framework of wheat [51], and those SNPs found within ã 5 Mb region of HSPs were reported to be candidate regions associated with TaHSPs. Favorable and unfavorable alleles were identified for all marker-trait associations (MTAs). For each MTA, the allele with an increasing effect was deemed to be a favorable allele, while the alternate allele was deemed to be an unfavorable allele (i.e., allele with decreasing effect). For some traits, such as DH and PH, the allele with decreasing effect was deemed to be a favorable allele and the alternate allele was deemed to be an unfavorable allele (i.e., allele with increasing effect). In each accession, the frequency of favorable and unfavorable alleles was counted based on the number of trait-associated SNPs. The coefficient of determination (R 2 ) was determined between the number of favorable/unfavorable alleles and the relevant phenotype as a measure of the combined allelic effects of all trait-associated SNPs on the phenotype.

Supplementary Materials:
The following supporting information can be downloaded at https:// www.mdpi.com/article/10.3390/plants12081610/s1: Supplementary Table S1: Names and pedigree of the 180 genotypes used in this study; Supplementary Table S2: Phenotypic data on synthetic hexaploid wheats in control and heat stress treatments; Table S3a: Meteorological information of Bahawalpur (BWP) for normal planting and late planting during the two cropping seasons (2013/14 and 2014/15); Table S3b: Meteorological information of Sindh (SND) for normal planting and late planting during the two cropping seasons (2013/14 and 2014/15); Table S3c: Meteorological information of Islamabad (ISD) for normal planting and late planting during the two cropping seasons (2013/14 and 2014/15); Supplementary Table S4: Phenotypic data on synthetic hexaploid wheats in control and heat stress treatments; Supplementary Table S5: Genotypic data as a HapMap file of all synthetic hexaploid wheats used in this study; Figure