Genetic dissection of heat-responsive physiological traits to improve adaptation and increase yield potential in soft winter wheat

Climate change, including higher temperatures (HT) has a detrimental impact on wheat productivity and modeling studies predict more frequent heat waves in the future. Wheat growth can be impaired by high daytime and nighttime temperature at any developmental stage, especially during the grain filling stage. Leaf chlorophyll content, leaf greenness, cell membrane thermostability, and canopy temperature have been proposed as candidate traits to improve crop adaptation and yield potential of wheat under HT. Nonetheless, a significant gap exists in knowledge of genetic backgrounds associated with these physiological traits. Identifying genetic loci associated with these traits can facilitate physiological breeding for increased yield potential under high temperature stress condition in wheat. We conducted genome-wide association study (GWAS) on a 236 elite soft wheat association mapping panel using 27,466 high quality single nucleotide polymorphism markers. The panel was phenotyped for three years in two locations where heat shock was common. GWAS identified 500 significant marker-trait associations (MTAs) (p ≤ 9.99 × 10− 4). Ten MTAs with pleiotropic effects detected on chromosomes 1D, 2B, 3A, 3B, 6A, 7B, and 7D are potentially important targets for selection. Five MTAs associated with physiological traits had pleiotropic effects on grain yield and yield-related traits. Seventy-five MTAs were consistently expressed over several environments indicating stability and more than half of these stable MTAs were found in genes encoding different types of proteins associated with heat stress. We identified 500 significant MTAs in soft winter wheat under HT stress. We found several stable loci across environments and pleiotropic markers controlling physiological and agronomic traits. After further validation, these MTAs can be used in marker-assisted selection and breeding to develop varieties with high stability for grain yield under high temperature.


Background
Worldwide, wheat is grown on more than 218 million hectares of land and provides approximately 20% of dietary calorie needs [1]. Although, there has been a substantial increase in yield since the Green Revolution, the pace of increase in yield production is not predicted to match demand resulting from increase in human population and changing weather pattern [2,3]. High temperature (HT) stress is one of the major consequences of climate change and poses a serious threat to wheat production [4]. Global temperature has increased by 0.5°C in the twentieth century [5] and this warming trend is expected to continue up to 1.5-4.5°C by the end of twenty-first century, resulting in elevated daytime maximum (HDT) and nighttime minimum temperatures. (HNT) [6]. Post-anthesis heat stress is very common in wheat growing areas and can cause large reductions in grain yield [7]. Although, some researchers suggest that HDT and HNT cause damage of a similar magnitude to winter wheat [8], others report a stronger negative impact on yield of HNT compared to HDT [9].
Genetic improvement of yield under HT via direct selection is hindered by the quantitative nature of grain yield and large genotype by environment interaction. Optimizing carbohydrate partitioning using traits such as spike fertility (SF), internode partitioning, and spike organ partitioning is essential to overcome sink limitations and increase harvest index (HI) [10,11]. Although higher grain yield potential of wheat has been largely a result of increased HI under HT, increased biomass and crop adaptation are equally important. Physiological trait (PT) selection insures development of stress resilient genotypes with functioning metabolic activities including photosynthesis and respiration under HT [12]. Moreover, previous studies have reported strong correlations of these traits with grain yield. Therefore, PTs can serve as indirect selection tools to select superior genotypes from large numbers of breeding lines for stress environments and to compensate for a large genotype by environment interaction. This is important to wheat breeders since it can save substantial amounts of labor, time, and money and permits rapid screening of a large number of genotypes in relatively short time [13,14]. Selection of desirable physiological trait (PTs) associated with heat adaptation and combinations is essential for future improvement of crops and provides opportunities for optimizing genetic yield gain. A model proposed to improve yield in wheat under heat stress includes partitioning of assimilates, radiation use efficiency (RUE), and light interception (LI) [15]. Some of the candidate PTs associated with these components are well documented as being heat adaptive traits, including higher leaf chlorophyll content measured as SPAD (soil-plant analyses development) value, intact leaf greenness measured as normalized difference vegetation index (NDVI), membrane thermostability (MT), and canopy temperature (CT). While these PTs are good candidates for improving heat tolerance and yield potential in wheat, limited knowledge on their genetic basis prevent full exploitation [11,16]. Use of these indirect selection tools is limited in breeding programs due to their complex evaluation procedure [17]. Therefore, identifying novel genetic loci (QTLs) associated with PTs under heat stress and using them as selection tools can results in a cumulative genetic effect on yield, which is the basis of maker assisted physiological trait breeding [18].
Association mapping is a powerful approach that utilizes genetic diversity and historical recombination events to provide a high resolution map of trait-linked loci [19]. Although genome-wide association studies (GWAS) have been used to identify quantitative trait loci (QTL) in wheat for various simply inherited traits like disease resistance [20], currently, limited information are available for GWAS of complex PTs, particularly under HT. Recently, the International Wheat Genome Sequencing Consortium (IWGSC) published a full chromosome-anchored reference genome which allows more precise curation of marker trait associations (MTAs) identified by GWAS. In this study, GWAS was performed on 236 advanced soft red winter wheat accessions using 27,466 SNPs generated by GBS. The panel was phenotyped at two heat stress locations over three years. The objectives of this study were: i) to identify novel MTAs linked to NDVI, CT, SPAD and MT under HT and ii) to identify candidate genes for these MTAs and investigate their underlying function.

Phenotypic analyses
There was significant genotypic variation (P < 0.001) for all measured traits (Additional file 3), as expected given the diverse genetic backgrounds of the SWAMP. Environments (growing years and locations) and their interaction were all significant (P < 0.05) determinants of phenotypic traits except for MT (Additional file 3). Trait means and a summary in response to each environment are provided (Table 1). All traits had moderate heritabilities, ranging from 60% for MT to 35% for CT (Table 1).
Pearson correlation coefficients (r) among PTs and their relationship with GY, GN, SF, HI, SHI and TGW were calculated using the combined dataset. Genetic data, LD decay, and population structure Population structure analysis of the SWAMP was performed in our previous study using 27,466 high quality GBS-derived SNP markers (minor allele frequency; MAF > 0.05 and missing data < 20%) [10]. Briefly, these SNP markers were distributed throughout the A (9958, 36%), B (9,968,~36%) and D (6,954,~25%) genomes. A total of 686 SNPs where found on unplaced scaffolds and thus were classified as unmapped SNPs. Chromosome 2B had the highest number of SNPs (1960) and chromosome 4D had the lowest (571). The population structure analysis grouped the 236 SWAMP lines into three genetic demes containing 49, 144, and 43 lines respectively (Additional file 6). PC analysis revealed substantial admixture among lines in the SWAMP, with the first and second PC explaining only 4.7 and 3.1% of the total genotypic variance, respectively (Additional file 6).
LD was computed using the "LDcorSV" package in R to determine the approximate marker density required for GWAS [10]. The LD decay below the line of critical value (r 2 = 0.2) was estimated at 1182, 1920 and 2916 bp for ranges of 30,000, 40,000, and 50,000 bp, respectively, across the whole genome (Additional file 7). The magnitude of change in LD decay between a sample range of 30,000-50,000 bp was 1734 bp. Population structure was investigated to avoid false positive associations in GWAS (Additional file 6) [19].

Marker-trait association
The GWAS identified novel MTAs for all measured PTs and explained a large portion of phenotypic variances from 5 to 23% (Additional file 8 Co-localized MTAs controlling multiple PTs were detected in the study. Ten pleiotropic SNP markers on chromosomes 1D, 2B, 3A, 3B, 5A, 6A, 7B, and 7D were detected across different environments ( Table 3). Seven of them were associated with NDVIa and NDVIg indicating common MTAs for NDVI expressed during anthesis and grain filling period. SNP S7D_635578722 had a positive allelic effect for SPAD, NDVIa and NDVIg. SNP S2B_466014434 was associated with SPAD and MT  and had a negative allelic effect on both traits. S3A_ 609909640 was associated with SPAD and CT and had a negative allelic effect on CT and a positive allelic effect on SPAD (Table 3). Interestingly, we detected five significant MTAs on chromosome 3A (S3A_12554694 and S3A_12554700), 5A (S5A_590056740), 6B (S6B_ 149148874), and 7D (S7D_18808932) for PTs, which were associated with GY and other yield related traits in our previous study [10]. Seventy-five out of 500 MTAs were expressed in multiple environments and were considered as stable MTAs (Table 4). We identified 23 stable markers for SPAD on chromosomes 1B, 1D, 2B, 2D, 3B, 4A, 5B, 5D, 6B, 6D, 7A, and 7D with PVEs ranging from 13 to 19% (Table 4, Additional file 8). For MT, we identified 20 stable MTAs on chromosomes 1B, 1D, 2A, 2D, 3A, 3D, 5B, 5D, 6A, 6B, and 7B with PVEs ranging from 5 to 14%. Similarly, 17 stable MTAs were detected for CT with PVE ranging from 8 to 13%. For NDVIa and NDVIg, we identified 3 and 12 stable markers respectively that were unique to corresponding growth stages (Table 4).

Gene annotation
Functional annotation of all stable MTAs was carried out using the IWGSC v1.0 reference genome sequence assembly. Forty-one out of 75 stable MTAs associated with PTs were anchored within functional genes. These MTAs had a wide range of functional annotations and are potential candidate genes for QTLs of interest (Table 5). Candidate genes were further investigated using past literatures to understand their possible functions. We discovered that these candidate genes encode different classes of proteins including F-box family proteins, RNA-binding proteins, disease resistance protein and protein kinase that have suggestive roles in response to biotic and abiotic stresses ( Table 5). In addition, gene annotation was also carried out for all the significant MTAs (Fig. 2). Interestingly, we identified several MTAs (linked to different traits) in different chromosomes that had common genes with exact same annotation ( Fig. 2 and Additional file 8).

Discussion
Growth and development of wheat is very sensitive to HT during anthesis and grain filling [7]. The impact of HDT and HNT on wheat growth and grain yield is well documented in many studies [7,9,21]. In this study, the SWAMP was evaluated in two heat prone environments of southeast USA with a goal to identify significant MTAs for use in breeding to improve adaptability and optimize yield potential. Short episodes of HDT (> 30°C) and HNT (> 21°C) were common during from anthesis to grain filling period in both environments (Additional file 1), thus, the MTAs identified in this study can provide useful information to understand the genetic bases of PTs under HDT and HNT.
HT reduces leaf area index and increases senescence rate which subsequently impairs photosynthesis and reduces grain yield [15,22]. Measuring chlorophyll content using SPAD, as a proxy for the entire photosynthetic complex, indicates photosynthetic potential. Higher expression of SPAD values during reproductive stages in wheat have been associated with heat tolerance resulting in higher grain yield potential [23][24][25]. In this study, genotypes showed significant genotypic variation in SPAD values with moderate broad-sense heritability (0.49) ( Table 1). Previous studies reported similar heritability for SPAD values [20]. Pearson's correlation showed strong positive relationship of SPAD with GY, SF, GN, HI, SHI, TGW, and MT. The result was supported by  (Table 4). Stable MTAs were located in chromosomes 1B, 1D, 2B, 2D, 3B, 4A, 5B, 5D, 6B, 6D, 7A, and 7D (Fig. 1a). In this study, SPAD shared common MTAs with MT (S2B_466014434), CT (S3A_ 609909640), and NDVIa and NDVIg (S7D_635578722) ( Table 3). These are potential new targets for multi-trait improvement and marker assisted breeding. Thirteen of these stable MTAs had functional annotation suggesting their involvement in abiotic stress including heat stress (Table 5). Three markers within 9 bp range on chromosome 6D (S6D_16178496, S6D_16178499, S6D_ 16178505) were within gene, TraesCS6D01G038900, with functional annotation of 2-oxoglutarate (2OG) and Fe (II)-dependent oxygenase superfamily protein. This gene has been reported to increase oxidative stress and reduce flower and pod number in canola when affected by heat stress [26]. One MTA on chromosome 4A (S4A_625244392) within gene TraesCS4A01G346600 had functional annotation as F-box protein. A gene (TaFBA1) encoding homologous F-box protein was reported to regulate gene expression and improve enzymatic antioxidant levels in response to heat stress in tobacco [27]. These proteins area involved in regulating many other biological processes, including biotic and abiotic stresses, floral development, embryogenesis, hormonal responses, and senescence [28]. Another MTA (S2D_574118879) within gene TraesCS2D01G469000 (GDSL-like Lipase/Acylhydrolase superfamily protein) has a predicted functional role in response to thermal stress in sorghum [29]. We also found MTAs associated with drought stress. Moreover, we found several other MTAs whose annotations suggest different roles including response to drought stress (TBC1 domain family member, SNF1-related protein kinase regulatory subunit beta-2) plant senescence (E3 ubiquitin-protein ligase ORTHRUS 2), salt stress (Glutamyl-tRNA (Gln)   amidotransferase subunit A), disease resistance (Disease resistance protein (NBS-LRR class) family), metabolism (Adenosine kinase-like protein) and antioxidant defense (Telomere repeat-binding factor like-protein) [30][31][32][33][34][35][36]. Our study identified novel MTAs associated with SPAD measurements in US soft wheat under HT conditions which contributes to a better understanding of the genetic basis of SPAD traits in wheat. Upon further validation, these MTAs can be used in future marker assisted breeding programs to overcome sink limitations, improve HI and ultimately increase yield potential. Another important trait that provides rapid measurement of crops to characterize the canopy for LAI and leaf greenness is NDVI. High leaf chlorophyll content at anthesis and the ability to retain greenness (delayed senescence) during grain filling stages is associated with higher heat tolerance [37]. Studies have confirmed that NDVI could be used to predict grain yield in wheat [14]. We found significant genotypic variation among SWAMP genotypes for NDVI at anthesis and grain filling period. NDVIa and NDVIg showed moderate broad sense heritability of 0.56 and 0.40 respectively (Table 1) which aligned with the results from previous studies [38]. NDVI showed significant positive correlation with GY, GN, HI, TGW, and MT at anthesis and grain filling stages. Some studies have reported a strong correlation of NDVI with wheat grain yield at any growth stages [14]. Others have reported a varied relationship of NDVI with grain yield, depending on growth conditions [39]. Here, we found several unique as well as common MTAs associated with NDVIa and NDVIg. We identified total 102 MTAs for NDVIa with PVEs ranging from 8 to 23% (Table 2, Fig. 1b). Six out of 95 MTAs were expressed in several environments indicating stability of these markers under different environments (Table 4). For NDVIg, 99 significant MTAs were detected with PVE ranging from 5 to 17% (Table 2, Fig. 1b) with 14 stable MTAs ( Table 2, Table 4, Fig. 1a). One MTA for NDVIa on chromosome 2B (S2B_717098540) is annotated as SAURlike auxin-responsive protein family (TraesCS2B01G522200) and has been reported to be upregulated under heat stress in Arabidopsis [40]. For NDVIg, MTA on chromosome 1B (S1B_677572998) within gene TraesCS1B01G467900 is annotated as Methyltransferase which has a predicted role of genetic or epigenetic regulation of heat responses in plants [41]. Another MTA (S5B_583295527) within gene TraesCS5B01G407600 had functional annotation for Myb family transcription factor-like protein. Myb transcription factor were reported to be significantly induced by heat treatment in rice and wheat and thus play important roles in response to high temperature stress [42,43]. Moreover, we found several other MTAs whose annotations suggest different roles including response to abiotic stress including drought stress (Multidrug resistance protein ABC transporter family protein, Myb/SANT-like DNAbinding domain protein, RING finger protein 13) [44][45][46][47]. Eight common MTAs were detected for NDVI at anthesis and grain filing period on chromosomes 1D, 3B, 5A, 6A, 7B, and 7D which can be potentially important targets for marker assisted selection. Three MTAs (S3A_699988530, S7B_701649275, and S7D_  [27,28,48]. Our study identified novel MTAs associated with NDVI which provides insight into the genetic basis of this trait in wheat at different growth stages in US soft wheat under HDT and HNT conditions. The selective permeability of plasma membrane is highly sensitive to heat stress affecting growth and development of a plant [49]. Therefore, MT is another important PT for understanding heat tolerance in wheat where lower expression of solute leakage leaf tissue indicates stability (tolerance) of cell membrane to elevated temperature [24]. We measured solute leakage from heat stressed plant tissue to estimate damage to cell membrane [50]. Overall, we found significant genotypic variation in SWAMP for MT with moderate broad sense heritability (0.60) ( Table 1) in agreement with previous studies [51]. MT was positively correlated with GY, SF, GN, HI, SHI, TGW, NDVIa, NDVIg and SPAD indicating that MT can be used as an additive component trait to improve yield potential in wheat under HT (Additional file 4). We identified 95 MTAs for MT with PVEs ranging from 5 to 15% (Table 2, Fig. 1b) indicating its quantitative nature. Twenty out of 95 MTAs were expressed in multiple environments indicating stability of these markers under different environments (Table 4). Twenty stable MTAs were located in chromosomes 1B, 1D, 2A, 2D, 3A, 3B, 3D, 5B, 5D, 6A, 6B, and 7B (Fig. 1a). Seven of these stable MTAs had functional annotation suggesting their involvement in abiotic stress including heat stress (Table 5). Two MTAs (S3D_590224603, and S3D_590224620) within 17 bp of each other were detected within gene TraesCS3D01G501200 with a functional annotation of protein kinase (Table 5). Protein kinases have been found to play a role in plant defense and adaptation responses and were reported to be upregulated by heat stress in durum wheat [52]. Another MTA (S1D_ 262475151) within gene TraesCS1D01G190700 (Heavy metal transport/detoxification superfamily protein) was reported to have an important role in growth and development of canola under heat stress conditions [53]. An MTA on chromosome 6B (S6B_42215716) within gene TraesCS6B01G063500 is annotated as peroxidase. Heat stress triggers the production and accumulation of harmful reactive oxygen species like hydrogen peroxide and their detoxification by antioxidant systems is important for protecting plants against heat stress [54,55]. A significant increase in the activity of peroxidase (antioxidant) under short term heat stress has been reported in heat tolerant genotypes indicating efficient antioxidative defense system in wheat [56]. Two MTAs within 23 bp on chromosome 2B (S5B_509105168 and S5B_509105191) were in gene TraesCS5B01G325000 with annotation as Fbox family protein.
Another parameter that has been frequently used to estimate heat tolerance in wheat is CT [57]. Lower canopy temperature indicates water status and transpiration rate in controlling temperature to avoid dehydration under stress [37,58]. In this study, there was significant genotypic variation in CT with moderate broad-sense heritability (0.35) ( Table 1). CT showed a significant negative correlation with GY, GN, SPAD, and MT. Genotypes with cooler canopies are presumed to have better root systems and retain chlorophyll content and membrane stability resulting in higher yield under high Fig. 1 Summary of GWAS. a genome-wide distribution of stable markers trait associations and b) range of percentage of variation explained for physiological traits in SWAMP. SPAD, soil-plant analyses development; MT, cell membrane thermostability; CT, canopy temperature (°C); NDVIa, normalized difference vegetation index at GS65; NDVIg, normalized difference vegetation index at grain filling temperature. We identified 110 MTAs for CT with PVEs ranging from 8 to 13% (Table 2, Fig. 1b). Seventeen out of 110 MTAs were expressed in multiple environments indicating stability of these markers under different environments (Table 4). Stable MTAs were located in chromosomes 2B, 5B, 5D, 6B, and 7B (Fig. 1a). Ten of the stable MTAs were located in five genes with annotated functions. An MTA on chromosome 5B (S5B_ 610295429) is within gene TraesCS5B01G436700 which is annotated as a lipid transfer protein. Lipid transfer proteins are low molecular weight proteins that are involved in many biological roles, such as anther development, different signaling pathways and heat stress both at the seedling and the grain-filling stages [59]. One MTA (S5B_621237427) on chromosome 5B was found in gene (TraesCS5B01G448700). This gene is annotated as mitochondrial transcription termination factor-like protein and is reported to be involved in heat tolerance in Arabidopsis [60]. Another MTA detected on chromosomes 5B (S5B_601343966), has gene annotation for Zinc finger CCCH domain-containing protein 16 (TraesCS5B01G425500). This protein was found to be over expressed upon high temperature stress in bread wheat [61]. We found several other MTAs whose annotations suggest different roles including abiotic stress response (BTB/POZ and MATH domain-containing protein 2), senescence (Maintenance of telomere capping protein 2, Protein Phosphatase 2C) and salinity stress (Peptidase M50 family protein) [62][63][64].
In summary, we detected 500 MTAs located in different chromosomes out of which 81 MTAs were linked to the same trait in multiple environments (suggesting stability) and ten MTAs linked to multiple traits (suggesting pleiotropy) ( Table 3-4). Notably, MTAs associated with multiple PTs within different genomic regions had the same functional annotation (Fig. 2; and Additional file 8). For instance, 13 MTAs for SPAD, MT, CT, NDVIa and NDVIg were annotated as F-box family proteins (Fig. 2). Similarly, the genes annotated as zinc finger proteins harbored MTAs for SPAD, CT and NDVIa. This result indicated the likely gene families that are important for physiological traits to improve yield potential under heat stress.
Some MTAs associated with PTs had pleiotropic effects with GY and other yield related traits (Additional file 9). These MTAs have been described in our previous studies. In a previous study [10], we found two MTAs (S3A_ 12554694 and S3A_12554700) within 6 bp were associated with TGW and GY. In this study, we found the same MTA associated with MT indicating that plasma membrane thermostability may contribute to increased TGW and GY under heat stress condition. Another MTA on chromosome 7D (S7D_18808932) also had a pleiotropic effect on MT and TGW. MTAs with pleiotropic effect on CT and TGW were detected on chromosome 6B (S6B_ 149148874) under multiple HT environments. One MTA (S5A_590056740), associated with HI in our previous study [10] was also linked to NDVIa and NDVIg. The co- Fig. 2 Potential candidate gene functions harboring SNPs affecting physiological traits under heat stress. The traits and count of marker-trait associations (for two or more traits) located within genes that have the same gene annotation is shown inside a bar. SPAD, soil-plant analyses development; MT, cell membrane thermostability; CT, canopy temperature (°C); NDVIa, normalized difference vegetation index at GS65; NDVIg, normalized difference vegetation index at grain filling localization of observed MTAs for HI with NDVIa and NDVIg suggests that HI is highly dependent upon greenness at anthesis as well as grain filling period under heat stress. This might indicate that higher photosynthetic reserves at anthesis can later be translocated to the developing grain.
In recent years, rapid progress in ground-based and unmanned aerial high-throughput field phenotyping have resulted in a variety of non-invasive imaging techniques which can lead to significant improvements in precision and speed of phenotyping for PTs in large plant populations with high resolution and high precision. These platforms can have large impact on validating as well as finding new genetic loci that are relevant to heat tolerance.

Conclusion
A large number of MTAs were identified in soft red winter wheat under the environments with HDT and HNT conditions. The MTAs were detected for four PTs: SPAD, MT, CT, NDVI for which there is reasonable evidence of being heat adaptive. We found several stable loci across environments and pleiotropic markers controlling multiple traits among PTs and also associated with yield-related traits under hot environments. We identified candidate genes affecting several important biological processes in plants including response to heat stress. Identifying regulatory loci associated with traits like PTs and yield-related traits can assist in developing ideotypes that can maximize the amount of assimilates and conversion of enhanced carbon capture and biomass growth for improving yield potential. Further validation of these MTAs in different controlled environmental conditions is required before they can be used extensively in marker assisted selection and breeding for heat tolerance in wheat.

Plant materials and experimental design
Field evaluations were conducted on 236 advanced genotypes of a soft red winter wheat association mapping panel (SWAMP) that are well adapted to the warm and humid south and southeastern regions of the USA. These lines were developed by public and private soft wheat breeding programs in the south and southeastern USA and the list is available in the NCBI database with accession number PRJNA578088 (https://www.ncbi.nlm.nih.gov//bioproject/ PRJNA578088). The seeds were collected from different soft wheat breeding programs (University of Arkansas, University of Georgia and Louisiana State University) from in the south and southeastern USA. The SWAMP was evaluated for PTs in five trials at the two heat stressed locations in Florida: Citra and Quincy. Citra had more frequent episodes of HT (> 30°C) during the grain filling period than Quincy (Additional file 1). In Citra, the SWAMP was phenotyped for three growing seasons Traits assessed in each year and used in GWAS can be found in (Additional file 2). All yield trials were planted in six row plots (3 m length × 1.5 m width) at the seeding rate of 100 kgh − 1 . The SWAMP was planted in randomized augmented block design [65] in all trials with 236 unreplicated entries and three repeated check varieties (SS8641, AGS2000 and Jamestown).

Trait measurement
Four physiological traits (PTs) including CT, SPAD, NDVI, and MT were measured at different time points. CT was measured three times during grain filling period (Zadoks stages 67, 72, and 77) using a handheld infrared thermometer (Fluke 572-2 IR thermometer, Fluke Corporation, Everett WA) on sunny days when the temperature reached the daily high between 1300 and 1500 h [66]. CT data were taken from the same side of each plot at 50 distance from the edge and approximately 50 cm above the canopy at an angle of 30 o to the horizontal. The average of three time point readings was used for the association analysis. Chlorophyll content was measured on flag leaves from eight random main stems for each plot at anthesis plus seven days (Zadoks stage 72) using a handheld self-calibrated SPAD chlorophyll meter (Minolta SPAD-502 Spectrum Technologies Inc., Plainfield, IL, US). The SPAD-502 instrument provides a convenient means of assessing relative leaf chlorophyll concentration. The chlorophyll content was measured on intact flag leaves one third of the way from the base of the abaxial surface. The average of eight readings was used for further statistical analysis.
SPAD Chlorophyll meter data were taken on the same day or the closest possible day coinciding with CT and MT. NDVI was measured at anthesis (NDVIa; Zadoks stage 65) and grain filling stage (NDVIg; Zadoks stage 77). A GreenSeeker handheld crop sensor (Trimble Navigation Limited 935 Stewart Drive Sunnyvale, California 94,085) was used for collecting NDVI readings. The GreenSeeker handheld crop sensor was hold 50 cm above the canopy facing the center of the bed. A 30-40 NDVI readings were recorded/plot and the mean value of those readings represented the NDVI value of the respective plot. To determine MT, flag leaves were collected from ten random main stems at anthesis plus seven days (Zadoks stage 72) from each plot. One cm diameter leaf disks from each leaf were extracted from midway between the base and the tip of the leaf blade using a leaf disc puncher and placed in glass vials containing 20 ml deionized water. The vials were placed in shaker for 24 h at room temperature to ensure diffusion of electrolytes. After 24 h, electrolyte leakage was measured using a conductometer (Thermo Scientific Orion Star A212) followed by autoclaving the vials (0.10 MPa pressure, 121°C for 15 min) to release all the electrolytes from plant tissue. The vials were placed in shaker for 24 h and electrolyte leakage was measured again. MT was expressed in percentage units as the reciprocal of relative leakage [50].
where T 1 is the conductivity reading after heat treatment, and T 2 is the conductivity reading after autoclaving. Grain yield (GY) and yield-related traits including grain number (GN), harvest index (HI), thousand grain weight (TGW), spike fertility (SF), and spike harvest index (SHI) were also calculated to determine correlation among traits. The details of the calculation of these traits is described in previous study [10].

Phenotypic data analysis
Combined analysis of variance (ANOVA) was conducted assuming a mixed linear model. The 'lme4' package [67] and the R software program (v3.5.1, R Development Core Team) were used to calculate best linear unbiased estimates (BLUEs) assuming a fixed genotypic effect (all other effects were random): where the phenotypic response (Y ijk ) is a function of the overall mean (μ), ith genotype (G i ), jth environment, genotype-environment interaction (GE ij ), kth block (B k ) nested within the jth environment (E j ), and the residual error (ε ijk ). BLUEs for combined as well as individual locations were also calculated and therefore will be discussed hereafter as BLUEC (BLUE values estimated from Citra), BLUEQ (BLUE values estimated from Quincy) and BLUEA (BLUE values estimated from all environments). Broad sense heritability was calculated assuming random genotypic effect (all other effects were random) and was obtained by: where H 2 , broad-sense heritability estimate; σ 2 G , genetic variance; σ 2 G × E , genotype-by-environmental variance; σ 2 e, residual variance; n, number of environments; and r, number of replications.
Pearson's correlations were calculated from BLUEs in R using the "corrplot" package (v3.5.1, R Development Core Team) and used to determine the direction and magnitude of measured trait associations. Associations between traits were also explored in principle component (PC) biplot analysis using the package "factoextra" in R [68].

Genotyping
The detail of genotyping process, SNP discovery and filtering criteria has been described in previous study [10]. In brief, DNA was isolated from fresh, young leaves using a modified Cetyl trimethylammonium bromide (CTAB) protocol [69]. The GBS libraries were prepared using MspI and PstI-HF restriction enzymes and pooled together in 96-plex and sequenced in an Ion Torrent Proton sequencer (Thermo Fisher Scientific, Waltham, MA, USA) at the USDA Central Small Grain Genotyping Lab, Kansas State University, Manhattan, KS, USA. Prior to analysis, 80 poly-A bases were appended to the 3′ end of all sequencing reads so that TASSEL 5.0 would attempt to use reads shorter than 64 bases rather than discarding short reads. SNP calling was performed in TASSEL v5.0 GBS v2.0 discovery pipeline [70] and aligned to the Chinese Spring wheat (RefSeq v1) genome sequence [71] using the default settings of BWA (version 0.6.1). The markers were filtered based on the criteria of minor allele frequency (MAF > 5%) and missing values (< 20%).

Linkage disequilibrium, population structure and GWAS analysis
Linkage disequilibrium (LD) and population structure analysis of the SWAMP has been described in detail previously [10]. Briefly, "LDcorSV" package [72] in R (v3.5.1, R Development Core Team) was used to estimate LOESS (Locally weighted scatterplot smoothing) regressions of mean r 2 (coefficient of LD) between pairs of SNPs sampled at the range of 30,000, 40,000, and 50, 000 bp. The intersection between critical value (r 2 = 0.2) and LOESS line was considered as the distance beyond which LD starts to decay. Population structure was observed using Bayesian information criterion (BIC) score provided by discriminant analysis of principal components (DAPC, "adegenet" package, R Development Core Team 2013) [73] to determine the optimum number of demes supported by the results. The principal component analysis was performed using "prcomp" ("stats" package) to investigate the genetic differentiation among and within demes.
GWAS was performed in three BLUE datasets (BLUEC, BLUEQ, BLUPA) for each trait to identify significant MTAs in SWAMP using Fixed and random model Circulating Probability Unification (FarmCPU) model [74,75] executed in the Genome Association Prediction Integrated Tool (GAPIT) package in R software package [76]. The first three principal components were used as covariates by observing model fit in Q-Q (quantile-quantile) plots, and kinship was determined from FarmCPU [74]. A uniform value of -log10(p) = 4.00 (p = 9.99 × 10 − 4 ) was used as the cut-off to define significant MTAs based on Q-Q plots [77,78]. Candidate genes associated with significant MTAs and their annotation were identified using Chinese Spring wheat reference genome (IWGSC RefSeq v1.0) [71]. The putative genes were further investigated in past literature to determine their association with phenotypic traits under heat stress.