Genome-Wide Association Analyses Identify QTL Hotspots for Yield and Component Traits in Durum Wheat Grown under Yield Potential, Drought, and Heat Stress Environments

Understanding the genetic bases of economically important traits is fundamentally important in enhancing genetic gains in durum wheat. In this study, a durum panel of 208 lines (comprised of elite materials and exotics from the International Maize and Wheat Improvement Center gene bank) were subjected to genome wide association study (GWAS) using 6,211 DArTseq single nucleotide polymorphisms (SNPs). The panel was phenotyped under yield potential (YP), drought stress (DT), and heat stress (HT) conditions for 2 years. Mean yield of the panel was reduced by 72% (to 1.64 t/ha) under HT and by 60% (to 2.33 t/ha) under DT, compared to YP (5.79 t/ha). Whereas, the mean yield of the panel under HT was 30% less than under DT. GWAS identified the largest number of significant marker-trait associations on chromosomes 2A and 2B with p-values 10−06 to 10−03 and the markers from the whole study explained 7–25% variation in the traits. Common markers were identified for stress tolerance indices: stress susceptibility index, stress tolerance, and stress tolerance index estimated for the traits under DT (82 cM on 2B) and HT (68 and 83 cM on 3B; 25 cM on 7A). GWAS of irrigated (YP and HT combined), stressed (DT and HT combined), combined analysis of three environments (YP + DT + HT), and its comparison with trait per se and stress indices identified QTL hotspots on chromosomes 2A (54–70 cM) and 2B (75–82 cM). This study enhances our knowledge about the molecular markers associated with grain yield and its components under different stress conditions. It identifies several marker-trait associations for further exploration and validation for marker-assisted breeding.


INTRODUCTION
Durum wheat (2n = 28, AABB, Triticum turgidum L. ssp. durum) is the most commonly cultivated form of allotetraploid wheat, and is grown on 8% of the world's wheat area (FAOStat, 2016). It originated in the Mediterranean region, and is used to make pasta and semolina products (Ren et al., 2013). Approximately 75% of durum wheat is still grown in the Mediterranean basin in irrigated and rainfed environments, which contributes to 50% of the worldwide production (Li et al., 2013;Kabbaj et al., 2017). Climate change will significantly impact the Mediterranean region; temperatures are predicted to increase by 3-5 • C, and annual precipitation is likely to decrease by 4-27% during the cropping season (Flato et al., 2013). The frequency and duration of dry spells and heat waves are also likely to increase in dryland areas (Parry et al., 2005;Bates et al., 2008).
Terminal drought and heat stresses negatively affect wheat grain weight and yield (Araus et al., 2003;Slafer et al., 2005). Dissecting the genetic bases of durum wheat responses to drought and heat stress is a prerequisite for breeding future genotypes (Graziani et al., 2014), and can be accomplished through the complementary approaches of association mapping and QTL mapping (Zhu et al., 2008). Association mapping, which is based on linkage disequilibrium (LD), is a powerful approach to genetic mapping that provides high resolution of detected loci, due to the presence of high genetic diversity and historic recombination of alleles in the assembled association mapping population (Sukumaran and Yu, 2014). A new multi-parent approach to map quantitative trait locus (QTL)-Multi-parent advanced generation inter-cross (MAGIC)-is an alternative to traditional linkage mapping but has less genetic diversity than the diverse association mapping panel (Milner et al., 2016). Several association mapping studies have been conducted to dissect the genetic basis of grain yield in durum wheat (Mengistu et al., 2016;Kidane et al., 2017). Maccaferri et al. (2011) evaluated a collection of 189 elite durum wheat lines in 15 environments with varying water availability during the cropping cycle. They identified 56 markers that explained 3.5-4.2% of the variation in grain yield, but the number of marker-trait associations (MTA) under drought stress were less in number compared with irrigated conditions. Several studies have focused on genetic diversity and molecular characterization of durum wheat landraces using different markers systems, but not many have used the DArTseq marker system in durum wheat (Yildirim et al., 2011;Kabbaj et al., 2017;Monostori et al., 2017).
Genetic improvement of durum wheat for drought and heat stresses can be achieved by direct or indirect selection for yield in target environments, or environments like target environments. Direct selection involves selecting for yield, whereas indirect selection-i.e., physiological breeding-selects for yield components or other associated traits such as canopy temperature, normalized difference vegetation index, grain number, and thousand grain weight (TGW; Araus et al., 2008;Reynolds et al., 2009b;Tuberosa, 2012;Reynolds and Langridge, 2016). The difficulty is in knowing which trait combinations should be selected to produce stable high yielding genotypes under varying environmental stress conditions (Habash et al., Abbreviations: BLUP, best linear unbiased prediction; DTA, days to anthesis; DTH, days to heading; DTM, days to maturity; GNO, grain number m −2 ; GWAS, genome-wide association study; IWIN, International Wheat Improvement Network; IWIS, International Wheat Information System; LD, linkage disequilibrium; MTA, marker-trait association; NDVI, normalized difference vegetative index, NDVIllg, NDVI at grain filling stage; NDVIvg, NDVI at vegetative stage; PH, plant height; SSI, stress susceptibility index; TOL, stress tolerance; SNP, single nucleotide polymorphism; STI, stress tolerance index; TGW, thousand grain weight; YLD, grain yield. 2009). It is therefore essential to understand the genetics and gene action of these traits. This study phenotyped a durum panel under yield potential (YP), drought stress (DT), and heat stress (HT) conditions. We studied the genetic diversity of the panel and identified molecular markers associated with grain yield and its components under YP, DT, and HT. In addition, common genomic regions (QTL hotspots) associated with grain yield (YLD), TGW, and grain number m −2 (GNO) under YP and DT, and markers for stress indices under DT and HT, were identified.

Germplasm
The 208-line durum wheat panel used in this study is a subset of the 15,000 CIMMYT gene bank accessions that were previously evaluated and characterized for use in breeding for heat and drought tolerance. The panel consisted of durum lines from different International Wheat Improvement Network (IWIN) nurseries: 2IDYN, 3IDYN, 15IDYN, 33EDUYT, 34IDSN, and 24EDYT-SA (

Phenotyping
Phenotyping was conducted at the Campo Experimental Norman E Borlaug in Cd. Obregon, Sonora, Mexico, under YP, DT, and HT conditions. Planting date and irrigation schedule were managed to create different stress conditions ( Table 1). The experiment was planted in an alpha lattice design with two replications, 2 m plots, 2 rows per plot in raised bed system 75 cm wide, with a 5 g m −2 seed planting density, and flooded irrigation in critical developmental stages. The genotypes were arranged in different "blocks" for alpha lattice design based on days to flowering. Irrigation frequency was the same under YP and HT to avoid the confounding effect of drought under HT conditions. N application was dependent on moisture availability and varied from 150 to 200 kg ha −1 . Fungicides and pesticides were applied to control local diseases and pests. Over two cropping seasons (2014-15 and 2015-16), the following traits were measured in accordance with established protocols (Pask et al., 2012): YLD, GNO, TGW, days to heading (DTH; under HT), days to anthesis (DTA; under YP and DT), days to maturity (DTM), plant height (PH), and normalized difference vegetation index (NDVI) at vegetative (NDVIvg), and grain filling stages (NDVIllg). Based on our previous research (data not shown), DTH and DTA are highly correlated under YP and DT conditions, but DTH and DTA might occur at the same time under HT conditions, which makes it difficult to measure DTA under HT.

High-Throughput Genotyping
A modified CTAB (cetyltrimethylammonium bromide) method (Saghai-Maroof et al., 1984) was used to extract genomic DNA from fresh leaves collected from the 208 entries. DNA quality and concentration were determined by electrophoresis in 1% agarose gel. High-throughput genotyping was conducted in 96

Data Availability
The genetic and phenotypic data used in the present study are available at http://hdl.handle.net/11529/11053. We used the consensus map from diversity arrays; physical positions of the markers are available in the above link.

Linkage Disequilibrium (LD)
Understanding the LD pattern in germplasm is important for selecting the marker density required for GWAS and for defining identified QTL regions (Siol et al., 2017). We computed LD decay using the open source R package, "sommer" (Covarrubias-Pazaran, 2016).

Population Structure Analysis
Population structure of the durum panel was inferred using the model-based clustering method implemented in STRUCTURE software (Pritchard and Przeworski, 2001;Falush et al., 2007), along with the delta K approach to statistically test the results (Evanno et al., 2005). We used 1,300 random SNPs, at least 5 cM apart, to estimate the population structure. Simulations were run by inferring K from 2 to 10, with 20,000 iterations and 5,000 burn-ins. The results were entered into the structure harvester (http://taylor0.biology.ucla.edu/structureHarvester/) to obtain the panel's delta K statistics (Earl and vonHoldt, 2012). Principal component analyses of the SNP data and neighbor joining tree construction were conducted in MEGA7 software (Kumar et al., 2016) and results from different approaches were compared to deduce population structure.

Statistical Analysis and Stress Indices
Analyses of variance were performed and best linear unbiased predictions (BLUPs) were obtained using META-R software (Vargas et al., 2013). When estimating BLUPs, genotypes, genotype-by-environment interactions, and environments were considered as random factors, while location, block, and replication were fixed factors. In addition, BLUPs for YLD and other traits were calculated for each treatment using flowering time as a covariate. We also calculated BLUPs for irrigated environments (YP and HT combined), stressed environments (DT and HT combined), and combined analysis of all environments (YP, DT, and HT). Pearson correlation coefficients (r) among different traits and locations were calculated using the cor() function in R, and the corrplot() package was used to plot the results. Broad sense repeatability (H 2 ) was estimated using: where σ 2 g is the genotypic variance, σ 2 ge is the genotype by environment interaction variance, σ 2 e is error variance, r is the number of replications, and l is the number of environments.

Genome-Wide Association Study (GWAS)
MTA analyses were conducted using the TASSEL 5.2.38 software, where generalized linear models and mixed linear models were fitted using SNP data, population structure matrix (Q), kinship matrix (K), coefficient of parentage matrix, and principal components. Principal components and Q matrix were fitted as fixed effects, while coefficient of parentage matrix and K matrix were fitted as random effects (in different combinations) to estimate the best model for each trait (Yu et al., 2006;Zhang et al., 2010;. Model fitting and the best model for each trait was based on the quantile-quantile (Q-Q) plots (Sukumaran et al., 2012;Sukumaran and Yu, 2014). We used a GWAS threshold of -log(p) = 3 to declare significant associations, which was determined based on the Q-Q plots and distribution of p-values. The Bonferroni correction has a more stringent threshold but-when tested-it did not result in many significant MTAs, so it was not followed.

Comparison of MTAs and Loci through Blast Searches
We compared significant MTAs to detect association patterns, common loci and unique loci for environments and traits, traits per se, and stress indices. We also conducted comparative analyses of the most significant loci with the bread wheat genome assembly TGACv1 (Triticum aestivum; http://plants. ensembl.org/Triticum_aestivum/Info/Index), and barley genome assembly Hv_IBSC_PGSB_v2 (Hordeum vulgare; http://plants. ensembl.org/Hordeum_vulgare/Info/Index) through BLAST searches to identify syntenic regions and candidate genes, if any.

Genetic Data, Population Structure, and LD Decay
From the raw genetic data of 53,911 markers, we obtained 6,211 SNPs, after removing monomorphic SNPs and checking minor allele frequency (MAF). SNPs with MAF <5% were excluded, as was missing data (SNPs with >20% missing data). The total length of the genetic map was 2,094 cM. The largest chromosomes were 1A (252 cM) and 1B (286 cM), while the smallest was 6B (84 cM). The average distance between SNPs was 0.34 cM. The length of the A genome was 1,087 cM when compared to the B genome of 1,007 cM. Chromosome 2B had the highest number of SNPs (707), while 4B had the lowest (271).
Population structure was inferred through different methods. The STRUCTURE algorithm and structure harvester results showed a delta K plot peak value of five ( Figure S1). The results from the phylogenetic tree and the principal component analysis plot with principal component (PC)1 vs. PC2, colorcoded from the STRUCTURE results, also showed five distinct groups (Figure 1). Group 1 comprised entries with JUPARE C 2001 pedigree, group 2 comprised entries that were crossed to Frontiers in Plant Science | www.frontiersin.org RASCON 37 or ALTAR 84, group 3 comprised entries from Egypt-Africa, group 4 were elite lines derived from several crosses, and group 5 were intermediates between African and elite lines. The group G3 was distant from other groups.
LD analysis used the "sommer" package in R software. The LOESS curve intercepted the line of critical value (r 2 = 2) at 6-8 cM, indicating that-for GWAS-at least one SNP is sufficient within 8 cM region in each chromosome ( Figure S2). We observed only for nine gaps >8 cM, which did not affect the GWAS results.
Agronomic Performance of the Panel YLD was highest under YP (5.79 t/ha; H 2 = 0.80), followed by DT (2.33 t/ha; H 2 = 0.47) and HT (1.64 t/ha; H 2 = 0.30). TGW also varied among the different environments at 44.4 g (H 2 = 0.87), 40.8 g (H 2 = 0.69), and 31.8 g (H 2 = 0.63), under YP, DT, and HT, respectively. The same trend was observed for all traits (Figure 2). Phenological measurements also followed a similar pattern; the crop had a shorter duration under HT (DTA = 55.5 days; DTM = 84.1 days), compared to DT (DTA = 71 days; DTM = 100 days) and YP (DTA = 76 days; DTM = 113 days    We also estimated 18 stress indices (SSI, TOL, and STI for three traits for two stress environments) and their associations with each other. STI and SSI were negatively associated in several cases (r = −1). The highest positive association (r = 1) was between TGW SSI vs. TGW STI under HT conditions ( Figure S3).

Yield Potential
Under YP, 121 MTAs were identified with p < 10 −03 . The trait variation explained by each marker (R 2 ) varied from 0.07 to 0.11 ( Table S2). The highest number of MTAs were identified for PH (30), followed by TGW (25). A comparison of the MTAs identified two loci associated with YLD, TGW, and GNO under YP on chromosomes 2A (54-70 cM) and 2B (78-82 cM). The locus on 2A was not associated with DTA whereas the locus on 2B was associated with DTA, DTM, PH, NDVIvg, and NDVIllg. A locus on chromosome 6A (34 cM) was associated only with YLD, whereas a locus on 5B was associated with PH and YLD, and one on 6B was associated with YLD, DTA, and DTM. DTA and DTM were highly correlated and had several common MTAs, meanwhile NDVIvg and NDVIllg had common MTAs with YLD in chromosome 2B (75-78 cM). MTAs for TGW and GNO that did not affect YLD were identified in chromosomes 5B (82 cM) and 7B (95 cM), respectively (Table S2).
GWAS identified fewer MTAs when DTA was used as a covariate in the BLUPs (compared to without DTA), but the consistent MTAs remained the same. The most consistent region associated with multiple traits was on chromosome 2B (74-82 cM), for YLD, GNO, TGW, DTA, DTM, and PH. The second most consistent region was on chromosome 2A (54-74 cM), for GNO, TGW, DTA, DTM, and PH ( Table 3). The highest number of MTAs was for PH, followed by DTM. MTAs for GNO and TGW on chromosome 2A (61-70 cM) were not associated with YLD, potentially indicating a compensation effect (Table S3).

Drought Stress
Under DT, 159 significant MTAs were identified for eight different traits (Table S4). Of these traits, DTA had the highest number of MTAs (44). Common regions for YLD, TGW, and GNO were not identified under DT, but MTAs for TGW were identified on chromosome 2A (65-70 cM), which was also associated with GNO and NDVIllg. A locus on chromosome 2B (70-82 cM) had significant MTAs for DTM, NDVIllg, DTA, GNO, PH, NDVIvg, and TGW. Meanwhile a locus on chromosome 4A (95-102 cM) was significantly associated with YLD and NDVIvg. MTAs for YLD were identified for six loci, of which five were independent of phenology: chromosome 1A (140-145 cM), 4A (95-96 cM), 5B (30 and 60 cM), 7B (134 cM) ( Table S4). Using DTA as a covariate in the calculation for BLUPs and GWAS under DT identified several consistent SNPs ( Table S5). The most consistent among them was on 7B (36-40 cM), which was associated with YLD, GNO, and NDVIllg ( Table 4). A locus on chromosome 2A (66-70 cM) that was associated with TGW under YP was also associated with TGW under DT. A locus on chromosome 2B (75 cM) was associated with GNO, DTA, PH, and NDVIllg under DT.

Heat Stress
Under HT, 112 MTAs were detected for the seven traits, excluding NDVIllg. TGW was the trait with highest number of MTAs (61) ( Table S6). Eight MTAs were detected for YLD; two of these, on chromosomes 2B (42 cM) and 4A (107-124 cM), were also associated with GNO and YLD. Use of DTA as a covariate showed   For each trait, the genetic position is represented by numbers. YLD, grain yield (t/ha); TGW, thousand grain weight (g), GNO; grain number/m 2 ; DTA, days to anthesis; DTM, days to maturity; PH, plant height (cm); normalized difference vegetative index at vegetative (NDVIvg) and grain filling (NDVIllg) stages.

Comparison of MTAs for YLD, GNO, and TGW under Different Environments
A comparison of the MTAs for YLD, GNO, and TGW identified a locus in chromosome 2B (74-82 cM) as the most common locus for these three traits under YP. Another locus common to GNO and TGW was located on chromosome 2A (61-70 cM). Several loci for YLD, which were independent of TGW and GNO, were For each trait, the genetic position is represented by numbers. YLD, grain yield (t/ha); TGW, thousand grain weight (g), GNO; grain number/m 2 ; DTA, days to anthesis; DTM, days to maturity; PH, plant height (cm); normalized difference vegetative index at vegetative (NDVIvg) and grain filling (NDVIllg) stages.
also identified on chromosomes 3B (17 cM), 7A (40 cM), and 7B (11 cM) ( Figure 4A). Under DT, no common loci were identified for YLD, GNO, and TGW, but independent loci were identified for all three traits ( Figure 4B). Similarly, no common loci were identified under HT, but a locus on chromosome 2B (81-82 cM) was associated with TGW and another at 74-75 cM was associated with GNO ( Figure 4C). For all three environments, the most common locus for various traits was on chromosome 2B (74-82 cM) ( Figure 4D). These results indicate that this region on chromosome 2B (74-82 cM) harbors genes for TGW and GNO. It is a region that is not in high LD, but is a QTL hotspot.

Markers for Stress Tolerance Indices
Three different stress tolerance indices (SSI, TOL, and STI) were calculated for YLD, GNO, and TGW by comparing DT and HT conditions with YP. GWAS detected common genomic regions and 201 MTAs for the stress indices (Tables S8, S9).

Drought Stress Indices
The number of MTAs detected for YLD-SSI, TGW-SSI, and GNO-SSI were 7, 11, and 13, respectively. For TOL, 142 MTAs were detected, while 28 MTAs were detected for STI under DT (Table S8). A comparison of these indicated several common regions ( Table 6) Table 6 and Figure S4).
We detected 93 MTAs for the traits in stressed environments (combined BLUPs of DT and HT conditions). Chromosome 4A (122-125 cM) harbored the highest number of MTAs for YLD. For TGW, the highest number of associations was on chromosome 2A (66-70 cM). As with irrigated conditions, most of MTAs identified for PH were on chromosome 4B (31-32 cM) ( Table S11). No common locus for YLD, TGW, and GNO was identified in this analysis ( Figure 5B). A comparison of MTAs detected by the combined BLUPs of three different environments indicated the most consistent locus for TGW and GNO on chromosome 2A (66-70 cM) ( Figure 5C and Figure S6). The same locus was associated with many other traits (Table S12). Two loci on chromosomes 2B (82 cM) and 5B (48 cM) were associated with TGW, but not with YLD. Similarly, loci on four chromosomes (2A, 5B, 7A, and 7B) were identified for GNO but not YLD (Figure 5C).
A comparison between the MTAs under irrigated, stressed, and combined analysis of three environments indicated a locus on chromosome 2A (66-70 cM) as the most common. Loci common to irrigated and combined analyses were identified on chromosomes 1B (6-8 cM and 71 cM) for YLD, 2B (82 cM) for TGW, and 7A (75 cM) for GNO ( Figure 5D).

QTL Hotspots for YLD, GNO, TGW, and Flowering Time
Further analysis of the GWAS results indicated two loci affecting YLD, TGW, and GNO in different environments (Table S13) ( Figure 6). The highest number of MTAs were detected on chromosome 2A (25%) followed by 2B (13%), with the least on chromosome 1A (1%) (Figure 7A). The variation explained by the markers varied from 7 to 27% ( Figure 7B). DTA had markers explaining the highest variation (>25%), followed by DTM (>24%), and PH (close to 20%) ( Figure 7C). The highest number of MTAs were detected for TGW-TOL combining all stress environments-under YP, DT, and HT environmentsfollowed by TGW and the lowest number of MTAs were detected for YLD tolerance indices estimated as TOL and STI ( Figure 7D).

Candidate Genes and Syntenic Regions
We further analyzed the SNPs in chromosomes 2A (54-70 cM) and 2B (75-82 cM) that were identified as QTL hotspots. There were 28 and 16 unique significant SNPs observed in the 2A and 2B region, respectively. BLASTN analysis of the SNP sequences on Ensemble genome browser for wheat and barley genomes indicated that several SNPs were related to transmembranes or were uncharacterized proteins. Of these, one SNP (100035706) was related to the gene DMAS1-A with a protein characterized as Deoxymugineic acid synthase 1. The physical location of this gene was in chromosome 4AS in bread wheat and in chromosome 4H in barley. In some cases, the hits were on homeologous chromosomes 2A, 2B, and 2D in bread wheat and it were observed in different chromosomes in Barley. We also BLASTed the SNP sequences against the UGRI durum Capelli v1 and Durum Strongfield v1, but the results were on chromosome contigs without further annotation.
conditions. The panel was phenotyped by adjusting planting dates to create YP, DT, and HT conditions in the field. Under HT, mean YLD of the panel reduced by 72%, which was larger than the reduction in YLD under DT (60%), when compared with YP conditions. Similar trends were previously reported in a study on a bi-parental population in similar environments (Pinto et al., 2010). In our study, HT was independent of DT, as irrigation was supplied under HT (similar to yield potential).
We observed moderate to high H 2 estimates for the traits, which is similar to the results achieved by  when they observed spring bread wheat under similar management conditions in Cd. Obregon, Mexico. YLD was negatively associated with DTA and DTM under YP, DT, and HT conditions, indicating a yield advantage for earlier genotypes (Millet et al., 2016). An important question now is whether grainfilling duration affects YLD under these conditions in durum wheat, as earlier studies have shown that grain filling duration and rate is reduced and starch synthesis inhibited in spring wheat when temperatures are above 30 • C . In our study, we did not observe association between grain filling duration and grain yield or TGW when different stress conditions and YP were compared. Another phenology trait, PH, was negatively associated with YLD under YP, but no significant association was observed under DT, and PH was positively associated with YLD under HT. This positive association of PH and YLD under HT may be related to higher biomass of the plants; as biomass and PH are positively associated .
NDVI was highly associated with YLD under HT, indicating that it is an effective selection tool for in-season prediction of wheat YLD (Tattaris et al., 2016). GNO was highly associated (>0.72) with YLD in all environments, whereas TGW was more significantly associated with YLD under HT than under YP and DT. This indicates that there is a tradeoff between TGW, GNO, and other yield components under varying environment conditions and further research is needed .
GWAS is the most popular approach for dissecting the genetic basis of complex traits (Risch and Merikangas, 2007), but this approach is prone to the detection of false positives due to confounding population structure or the effect of phenology genes (Yu et al., 2006;Zhang et al., 2010). Our study used the mixed model framework of Yu et al. (2006), with fixed and random effects to control false positives. Q-Q plots for multiple models were evaluated to select the best models for identifying MTAs . We also observed DTA and DTH were associated with grain yield in this panel. Using phenology as a covariate did not remove the co-localization of QTLs for flowering time and agronomic traits, indicating the strong association between these traits. GWAS was performed with BLUP values estimated with and without DTA/DTH as a covariate. We identified several MTAs for YLD and component traits co-localized with and without phenology genes. Flowering time and PH are associated with adaptation and agronomic performance of traits in several crops, and can be related to drought escape (Shavrukov et al., 2017). This is observed not only in wheat, fine tuning flowering time though earliness per se (Eps) genes (Zikhali et al., 2014Sukumaran et al., 2016), but also in barley, controlled by the EPS2/Eam6 gene with pleotropic effects on agronomic traits (Tondelli et al., 2014). Earlier studies on barley have reported that YLD is defined by the length of different sub-phases of vegetative, flowering, and grain filling stages (Francia et al., 2011).
Co-localizing MTAs were identified for YLD, GNO, and TGW, and individual MTAs were observed for each trait. We also identified MTAs for YLD that were independent of phenology (on chromosomes 6A at 34 cM for YLD under YP and 5B at 30 and 60 cM for YLD under DT). This indicates thatin addition to GNO and TGW-other yield components (e.g., plant density) are also important for increasing YLD (Sukumaran et al., 2015). In addition, the individual traits TGW and GNO can be manipulated, independently of YLD, as unique MTAs were observed for them under different conditions . We identified many genomic regions associated with the traits; some of them-on chromosomes 2A and 2Bwere common to both YP and DT conditions. As far as we know, this is the first comprehensive study focusing on MTAs and their interactions for YLD and its components under YP, DT, and HT conditions in durum wheat. There were fewer MTAs for NDVI, and some co-located with the agronomic traits, though a comparison with an earlier study on nitrogen use efficiency indicated that common QTLs for NDVI may exist on chromosome 3B (Quraishi et al., 2011;Monostori et al., 2017). An earlier study on durum wheat using 10 rainfed and 6 irrigated environments also identified a major YLD QTL on chromosome 2B with significant effects (Maccaferri et al., 2008). Another study used an association mapping approach in durum wheat and identified markers for YLD under different drought stress conditions (Maccaferri et al., 2011). A common genetic map could be used to compare these results with our own and identify common loci.
Common genetic markers for TGW and GNO were detected on chromosomes 2A and 2B under DT and YP conditions but, under HT, the markers associated with TGW and GNO were on chromosome 4B. This indicates that the genetic basis for YLD, TGW, and GNO were determined through a different mechanism in HT, compared to DT and YP conditions. An evolutionary study in wheat comparing wild emmer (T. turgidum subsp. dicoccoides) and durum wheat for TGW and embryo size identified a cluster of loci affecting TGW and grain shape in the long arm of chromosome 2AL, and a novel locus controlling embryo weight in chromosome 2AS (Golan et al., 2015). A BLAST search of the 2A and 2B hotspots indicated that they might be homeologous QTLs, as several SNP hits were observed in chromosomes 2A and 2B. The HvCEN gene with pleotropic effects on flowering time and agronomic traits previously found in chromosome 2H of barley (Francia et al., 2011;Comadran et al., 2012;Tondelli et al., 2014) might be syntenic to the 2A hot spot region (BLAST hits on 2AL:32101-32658; TGACv1 for bread wheat). We detected more MTAs for GNO under DT, though the number of markers for TGW were higher in HT.
Breeding for wide adaption is hindered by high genotype × environment interactions. In this study, we identified loci stable across DT and HT environments using three different stress indices. These stress indices had several common MTAs, which were also associated with YLD and its components. Further studies are required to understand and validate the effect of these markers. The markers on chromosomes 2B (83 cM), 3B (50 cM), 5A (150 cM), 6A (28 cM), and 6B (31 cM) are candidates for further exploration of drought stress tolerance. The loci on chromosomes 2B (74-83 cM) and 3B (68-83 cM) are prominent loci for heat stress tolerance, which may be similar to an earlier study for drought stress (Maccaferri et al., 2008).
More powerful statistical genetics tools would be needed to detect minor alleles associated with the traits. The panel used in this study had a group of lines from West Asia and North Africa, which was genetically different from other entries in the panel, but the low number of lines does not allow for identification of MTAs on minor alleles, if present. Pleiotropy is another aspect of the study where several traits were associated with QTL hotspots in chromosomes 2A and 2B. The locus on chromosome 2A (54-70 cM) affects multiple traits and stress indices indicating QTL hotspots (Figures S7, S8). This QTL region might need to be recombined by mutation or TILLING to see the effect of MTAs on individual traits (Uauy et al., 2009). The variation explained by the MTAs varied from 7 to 27%, indicating that SNP-based GWAS may need to be complemented by haplotypebased GWAS to explain the missing variation (N'Diaye et al., 2017). However, accurate detection of haplotypes based on general LD decay may not give promising results and may need to be substituted by LD block-based haplotypes.
This study contributes large number of MTAs in durum wheat for agronomic traits under YP, DT, and HT conditions. Stable loci across environments were identified that can be further explored for use in marker-assisted selection and gene discovery. Plants' responses to abiotic stresses are complex and it is essential to identify regulatory loci if these traits are to be manipulated. Climate change is predicted to have significant impacts in the Mediterranean region. To mitigate climate change, genomic technologies, and resources with genetic approaches need to be coupled through marker-assisted selection (Habash et al., 2009). Instead of conventional direct selection for YLD, genetic loci for yield components and associated traits should be identified to enable the manipulation of individual traits to obtain the cumulative gene action for enhancing genetic gains in wheat (Reynolds et al., 2009a;Reynolds and Langridge, 2016). Additionally, the recent availability of highquality genome assemblies for tetraploid wheat (Avni et al., 2017), coupled with the analysis of transcriptome profiles under DT conditions (Habash et al., 2014), will facilitate the identification and manipulation of causative loci governing yield stability across a broad range of environmental conditions.

AUTHOR CONTRIBUTIONS
SS and MR: Conceived and designed the study; SS and CS: Genotyped the panel; CS: Aligned the markers on the consensus map and reference genome; SS: Did the genetic and phenotypic analyses; SS: Wrote the manuscript; All authors reviewed the manuscript.

ACKNOWLEDGMENTS
This work was implemented by CIMMYT as part of the TRIGO project, made possible by the generous support of MasAgro, International Wheat Yield partnership (IWYP), and WHEAT CRP. Any opinions, findings, conclusion, or recommendations expressed in this publication are those of the authors and do not necessarily reflect the view of donors. We also acknowledge Dr. Hans Braun, Dr. Karim Ammar, and Dr. Tom Payne for their support to this study and Dr. Roberto Tuberosa for comments improving the manuscript. Editorial support was received from Emma Quilligan. Thanks are due to Araceli Torres Garcia, Nayeli Quiche Morales, Noel Pimienta Valenzuela, Dr. Jacinta Gimeno, and the staff of the Wheat Physiology lab, CENEB, Cd. Obregon, Mexico.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2018. 00081/full#supplementary-material        Table S1 | Description of the entries from the durum wheat panel.        Table S10 | GWAS results for the combined analysis of irrigated conditions (yield potential and heat stress).
Table S11 | GWAS results for the combined analysis of stressed environments (drought and heat stress).
Table S12 | GWAS results of the combined analysis of the traits from yield potential, heat stress, and drought stress.
Table S13 | Comparison between the MTAs identified for yield and yield components with that of stress indices.