Genomic Regions Associated with the Control of Flowering Time in Durum Wheat

Flowering time is a critical stage for crop development as it regulates the ability of plants to adapt to an environment. To understand the genetic control of flowering time, a genome-wide association study (GWAS) was conducted to identify the genomic regions associated with the control of this trait in durum wheat (Triticum durum Desf.). A total of 96 landraces and 288 modern lines were evaluated for days to heading, growing degree days, and accumulated day length at flowering across 13 environments spread across Morocco, Lebanon, Mauritania, and Senegal. These environments were grouped into four pheno-environments based on temperature, day length, and other climatic variables. Genotyping with a 35K Axiom array generated 7652 polymorphic single nucleotide polymorphisms (SNPs) in addition to 3 KASP markers associated with known flowering genes. In total, 32 significant QTLs were identified in both landraces and modern lines. Some QTLs had a strong association with already known regulatory photoperiod genes, Ppd-A and Ppd-B, and vernalization genes Vrn-A1 and VrnA7. However, these loci explained only 5% to 20% of variance for days to heading. Seven QTLs overlapped between the two germplasm groups in which Q.ICD.Eps-03 and Q.ICD.Vrn-15 consistently affected flowering time in all the pheno-environments, while Q.ICD.Eps-09 and Q.ICD.Ppd-10 were significant only in two pheno-environments and the combined analysis across all environments. These results help clarify the genetic mechanism controlling flowering time in durum wheat and show some clear distinctions to what is known for common wheat (Triticum aestivum L.).


Introduction
Flowering induction plays a pivotal role in the plant life cycle, affecting reproductive success and yield depending on the prevailing climatic conditions of the target environment. In durum wheat, heading and flowering times are critical stages in crop development as they play an important role in adaptation, yield potential and grain quality [1]. In addition, climatic stress during anthesis negatively affects pollen fertility [2,3]. Therefore, plant breeders need effective tools to predict flowering time in order to identify promising genotypes adapted to different environmental conditions.
Flowering time in wheat is controlled mainly by three groups of loci, two of which interact with environmental factors, namely photoperiod sensitivity genes (Ppd) and vernalization genes (Vrn) [4], while the third group of loci is defined as 'narrow-sense earliness' or 'earliness per se' (Eps) because these act independently from vernalization and photoperiod [5]. However, it is not certain whether Eps genes act independently of all environmental cues [6,7], since for instance Appendino and Slafer [8] showed that Eps genes could respond to temperature changes. Allelic variation at Ppd genes divides the temperate cereals into photoperiod-sensitive and photoperiod-insensitive classes, whereas differences in the Vrn alleles divide them into winter and spring classes.
Vernalization is the acquisition or acceleration of a plant's ability to flower after exposure to a certain amount of cold temperatures, and it is a strategic adaptation mechanism to postpone heading after the frost-prone winter months [9]. Natural variation in vernalization requirement is mainly associated in common wheat with allelic differences in the Vrn1, Vrn2, and Vrn3 genes. The Vrn genes regulate the transition from vegetative to reproductive phase in response to temperature [4] and thus determine the spring and winter growth habit. Different alleles respond differently to temperatures, but in general prolonged exposures to temperatures at least below 16° C are necessary to achieve the vernalization needs, and colder temperatures tend to accelerate this process, reducing the total number of days needed for flowering. This is one of the main adaptation system that allows winter wheat to survive at lower temperatures than spring wheat [10]. Winter wheat possesses recessive alleles at Vrn-A1, Vrn-B1, Vrn-D1 and Vrn-D5 loci [11,12], while spring wheat has dominant alleles at one or more of them [13]. The dominant allele of Vrn-A1 confers complete insensitivity to vernalization in spring growth habit and is epistatic to dominant alleles of Vrn-B1, Vrn-D1 and Vrn-D5, which confers low sensitivity to vernalization in facultative winter growth habit [11,[14][15][16]. Durum wheat also harbors Vrn-A1 and Vrn-B1 located on the long arms of chromosomes 5A and 5B [17,18]. Recent advances in wheat genomics have allowed for the cloning of Vrn-A1, Vrn-B1 and Vrn-D1 genes [19], and the development of functional SNP markers for their characterization.
Photoperiod response is another important factor influencing the initiation and length of flowering period. Natural variation in response to photoperiod is mainly determined by allelic differences in the Ppd1 gene, a member of the pseudo-response regulator (PRR) gene family [20]. Similarly to vernalization, this a major adaptation mechanism to delay flowering until after the short days of winter have passed to avoid risks of frost damage to the reproductive organs. Photoperiod-sensitive wheat plants initiate flowering only after long days with more than 13 hours of sunlight are perceived, while photoperiod-insensitive types can induce heading irrespective of daylength. Mutation at Ppd-1 locus enables wheat to become photoperiod insensitive and flower irrespective of day length. These mutations have been put under strong selection pressure in the past to enhance yield under certain climatic conditions, via the promotion of early flowering to avoid terminal climatic stresses late in the season. Particularly famous is the case of Norman Borlaug who was able to impose a very stringent selection of photoperiod-insensitive types via the use of shuttle breeding (two breeding cycles per year). This process resulted in widely adapted wheat cultivars that prompted the Green Revolution [21,22].
In durum wheat, photoperiod sensitivity is determined by Ppd-A1 and Ppd-B1 loci, located on chromosomes 2AS and 2BS, respectively [23], while photoperiod insensitivity results from mutations in any of the two Ppd-1 genes. By convention, alleles conferring photoperiod insensitivity are assigned an 'a' suffix (e.g. Ppd-A1a) [24]. In durum wheat, two large deletions within the Ppd-A1 gene designated as alleles 'GS-100' and 'GS-105', were reported to accelerate flowering, which led Wilhelm et al. [25] to conclude that these deletions are the likely causal basis of photoperiod insensitivity in tetraploid wheat.
While the knowledge on gene action to control flowering time in hexaploid wheat is quite accurate and routinely exploited in marker-assisted breeding [26], its application in durum wheat is yet not fully tested and several discrepancies have been found when deploying Vrn and Ppd markers [27]. For that reason, genome wide association study (GWAS) has been adopted to validate the loci associated with flowering time in durum wheat. A QTL associated with Ppd-A1a was shown to significantly reduce heading time in a RIL population derived from the cross 'Kofa' ('GS-100' allele) × 'Svevo' ('GS-105' allele), suggesting that these two alleles have different effects on photoperiod sensitivity in durum [1]. A genome-wide association scan for heading date from 27 field trials in the Mediterranean Basin and in Mexico resulted in the identification of 50 chromosomal regions [28], of which only Ppd-A1 and Ppd-B1 had been previously described.
To better understand the effect of different flowering genes in durum wheat, a GWAS study was conducted using both modern germplasm and landraces, and exposing these to extremely different climatic conditions, with the scope of triggering the expression of different loci and capture their genetic effects. The use of a selected set of contrasting environments was deemed novel to better qualify and define the effect of different loci.

Determination of phenological environments (PhEnv)
Plotting of climatic data ( Figure 2) showed sizable variation among the 13 environments for average temperatures and day length, with Lebanon off-season planted in summer having a clear distinction in daylength and temperatures. Similarly, the experimental sites along the Senegal River (Fanaye and Kaedi) had shorter daylength and higher temperatures. Terbol and Kfardan in Lebanon were the only sites where below zero temperatures were recorded. A clustering analysis of climatic data based on PCA ( Figure S1) revealed four main phenological environments (PhEnv) explaining 95% of the total climate variation divided as follows: PhEnv1 included five environments in Morocco

Analysis of variance for G, G x E, and G x PhEnv
Analysis of variance (Table S1) showed significant (p<0.01) differences for environment, PhEnv, genotype, genotype x environment, and genotype x PhEnv effects for all three tested traits. The GxPhEnv explained 68.7%, 80.7%, and 66.6% of the GxE for DTH, CGDD, and CDL, respectively. DTH ranged from 71 to 107 days with an average of 87 days across environments ( Figure S2). Variation for CGDD and CDL ranged from 1,049° to 2,019° C and 43,871 to 91,760 minutes, respectively ( Figure S2). There was strong genotypic variation among accessions for DTH ( Figure 3), with most of the modern lines flowering between 75-90 days, while landraces reached flowering only after 95-105 days. The high temperature and short photoperiod of PhEnv3 and PhEnv4 resulted in the shortest DTH (65 to 73 days), with the average min temperatures ranging between 14.4 and 21.6 °C, and the average max temperature from 31.0 to 35.5 °C. On the other hand, the short photoperiod and high temperatures prevented few modern lines and most of the landraces to flower at all. To better define the effect on Vrn and Ppd genes at different PhEnv, Table 1 was developed to indicate how these loci will be differently expressed in different PhEnv and between them. PhEnv1 in Morocco experiences only mild cold temperatures so would satisfy only weak Vrn requirements, but daylength increases in the transition to spring so full Ppd requirements can be met. PhEnv2 in Lebanon has the same Ppd conditions as PhEnv1, but the colder winter temperatures favour flowering also in genotypes with stronger vernalization requirements. PhEnv3 along the Senegal River does not experience any cold temperatures and its vicinity to the equator prevents significant changes in daylength during the season. Similarly, PhEnv4 does not experience cold days, and the photoperiod is actually shortening between early summer planting and early fall harvest.

Marker-trait associations
Association analysis identified 41 and 68 markers associated with DTH, 34 and 63 markers for CGDD, and 27 and 66 markers for CDL for landraces and modern lines, respectively ( Figure S3). Regression analysis confirmed 34 significant QTLs, 7 of which were in common between modern lines and landraces, 12 unique of landraces, and 15 unique of modern lines (Table S2 and Table S3). Of these QTLs, six were significant for all the three traits, whereas four were unique for DTH and CDL, nine for both DTH and CGDD as well as for DTH alone in modern lines. For the landraces, eight QTLs were common among all traits and whereas five were unique for DTH and CDL, four for DTH and CGDD, and two were specific for DTH. Based on theses combinations, QTLs associated to DTH and CDL were defined as ppd, to DTH and CGDD as vrn, and to DTH alone and QTL that overlap assigned as eps.

Flowering loci identified among landraces
Among landraces, one QTL was identified only in PhEnv1 conditions, eight QTLs in PhEnv2 and PhEn4 conditions, and two QTLs by PheEnv3, while six QTLs showed significant effects across environments (Table S2)
Several QTLs were disperse among PhEnv and across environments. For instance Q.ICD.Eps-03 (Chr2A) was significant in PhEnv2 and PhEnv3 in modern lines, while in landraces it was significant in PhEnv4 and across environments. Similarly, Q.ICD.Vrn-18 was significant in PhEnv2 in landraces while in modern lines it was prominent in PhEnv3 and PhEnv4. In contrast, Q.ICD.Eps-11, Q.ICD.Ppd-12 and Q.ICD.Vrn-13 were identified in the same PhEnv in both landraces and modern lines, with varying LOD (3.0 to 3.9), r 2 (0.04 to 0.23) and variance (2.0 to 30.3%). Among these, Q.ICD.Eps-11 was significant in PhEnv2 and Q.ICD.Ppd-12 and Q.ICD.Vrn-13 in PhEnv4. However, modern lines also showed significance in PhEnv3 for Q.ICD.Vrn-13. Q.ICD.Vrn-17 was consistently significant across all locations and explained up to 0.04 r 2 and 5.1 to 22% total variance in modern lines while in landraces it was significant only in the combined analysis across environments and accounted for 0.18 r 2 and 3.9% of the phenotypic variation.

Climatic effect on the control of flowering time
From an agricultural perspective, plants are considered better adapted to a specific region when they flower at the appropriate time [29]. Flowering time in durum wheat seem to involve more loci compared to bread wheat, with smaller effects, resulting in higher diversity for flowering time, a driving force for adaptation to the diverse environments of the Mediterranean region [30]. A total of 384 durum entries assessed for flowering time across 13 environments located at different latitudes and temperature regimes (Figure 1 and 3) confirmed a significant effect of all sources of variations for the three flowering traits considered: DTH, CGDD, and CDL. The accessions showed remarkable variation for these traits, with two main subgroups corresponding to landraces and modern lines. These results confirmed that temperature and photoperiod had a significant effect in determining the flowering initiation in addition to a number of genes promoting earliness per se (Eps genes). Four diverse pheno-environments (PhEnv) were determined based on PCA of climatic conditions at the 13 environments. The four PhEnv had contrasting climatic conditions, and their pairwise comparison allowed to distinguish the effects of Ppd and Vrn from those of the Eps loci (Table 1). Interestingly, germplasm tested in PhEnv4 (warm summer of Lebanon) resulted in the shortest average DTH at 62 days, since only germplasm without Ppd and Vrn requirement could flower. Similarly, PhEnv3 along the Senegal River also promoted early flowering time at 65-73 days, preventing germplasm with Ppd and Vrn requirements to flower. PhEnv1 of Morocco allowed for Ppd requirements to be met and mild Vrn, resulting in DTH ranging between 98 and 116 days. PhEnv2 represents Lebanese sites with long cold winter season, which highlighted the role of both Ppd and Vrn requirements, resulting in the longest DTH time from 122 to 137 days. Interestingly, CGDD also varied widely among PhEnv, with the longest season in number of days in Lebanon resulting in the lowest cumulative values (10,490 -11,650), while the other PhEnv overlapped. CDL at flowering increased 5,000 -10,000 minutes from PhEnv3 to PhEnv4, while it was similar in PhEnv1 and PhEnv4. These results suggest that the main climatic effects have a strong role in determining the time to flowering. Past studies have shown that latitude integrates a number of variables affecting wheat development, among which the most important are photoperiod [31], temperature [32] and their interaction [33]. Similar studies in durum wheat [22] suggested photoperiod and temperature associated with different locations were the most important variables in explaining phenological differences. Another research in spring durum wheat indicated that photoperiod and temperature together explained 77% of environmental variation between locations [34]. Our result provide only partial support to these findings, indicating that only up to 39% of phenotypic variation for flowering could be captured by PhEnv determined on the basis of climatic variables.

Known loci involved in the control of flowering time in durum wheat
The present study confirmed a role for four well established loci in controlling flowering time in durum wheat across the four PhEnv. Among them, three loci, Ppd-B1, Vrn1 and Vrn3 showed significant effects in landraces as well as in modern germplasm while Ppd-A1 had strong effect only in modern lines.
Q.ICD.Vrn-13 included the Vrn-A1 locus. It showed a high LOD value and accounted for significant variability (from 10 to 20%) for flowering time in modern lines and landraces in PhEnv3 (Senegal and Mauritania), and only in modern in PhEnv4 (Lebanon off-season), where no vernalization requirements can be met. This is in line with what can be expected since accessions carrying the recessive Vrn allele would not be able to flower under those conditions. Previous studies suggested Vrn-1 to be essential in the control of winter flowering types of common wheat [35][36][37], and durum wheat [22,27,[38][39][40][41][42]. However, results of the present research show only partial agreement with the literature, with a clear role of Vrn1 confirmed only when no cold temperatures occur (see PhEnv3 and 4 in Table 1), but substantially no effect in PhEnv where weak or strong vernalization requirements can be met (see PhEnv1 and 2 in Table 1). Since PhEnv1 and 2 represent true Mediterranean environments where spring durum wheat is primarily cultivated, our results suggest limited importance of Vrn-A1 in determining adaptation of durum wheat in its main area of production.
Earlier studies reported significant role of Vrn3 [4,43,44] suggesting its role in upregulating Vrn-1 above the threshold levels required for flower initiation. Q.ICD.Vrn-18 on Chr 7A herein reported likely corresponds to Vrn3. This QTL had a weak effect in modern germplasm in PhEnv3 and in landraces in PhEnv2. These results suggest a weak role for this locus in the control of flowering time in durum wheat, with no effect in environments where cold temperatures did not occur throughout the season (i.e. PhEnv4), probably due to the masking effect of Vrn-A1.
Locus Q.ICD.Ppd-07 overlap with Ppd-B1 and showed a significant effect only in PhEnv2 (Lebanon), where full photoperiod requirements are presumably met for both modern lines and landraces. Past studies reported the importance of Ppd-B1 in bread and durum wheat [22,28,30,[45][46][47]. Royo et al. [22] reported the strong effect Ppd-B1 in 35 spring durum wheat lines and an interaction between Ppd-B1 and Ppd-A1. Recently, Würschum et al. [30] evaluated European durum genotypes for heading time in five environments and identified Ppd-B1 as the major determinant of heading time in durum wheat, accounting for up to 26.2% of the variance. Our study observed confirmed an effect of Ppd-B1, but it accounted only for 5% of phenotypic variance in landraces as well as modern lines. This difference might be due to the wider genetic diversity assessed in our study as compared to what done previously.
Another major flowering gene reported by many researchers in durum and bread wheat is Ppd-A1. In the present study, Q.ICD.Ppd-21 on Chr 2A spans the Ppd-A1 locus. This was the most significant and stable QTL in our study, with effect in all PhEnv in modern germplasm, while it was monomorphic in landraces (no genetic variation available) as it can be expected. Earlier studies in durum wheat [48] also found strong effect of Ppd-A1 distributed in modern durum wheat lines, suggesting its exploitation started after the green revolution and was then further selected to increase adaptation. Wang et al. [49], Maccaferri et al. [50], and Royo et al. [22,27] also suggested a major role of Ppd-A1 over Ppd-B1 in durum wheat.

Identification of novel loci involved in the control of flowering time in durum wheat
The major loci known to be involved in the control of flowering time explained only partially the remarkable variation observed within the tested panel. Therefore, the association between phenotypic and genotypic effects was investigated to identify a total of 30 loci not overlapping with known major flowering genes across the four PhEnv, including four loci in common between landraces and modern lines.
Q.ICD.Eps-03 on chr. 2A was identified in PhEnv2 and PhEnv3 in modern germplasm for DTH and in PhEnv4 across environments in landraces for all the traits (DTH, CGDD and CDL). Because its role was significant across contrasting environments, it is likely that its impact on flowering time is independent from climatic conditions (Eps), explaining up to 10% of the phenotypic variation. Giunta et al. [51] also reported a QTL on chr. 2A in durum wheat explaining 16.9% of the phenotypic variation under long day, which falls within our QTL confidence interval (CI) when aligned to the Svevo genome [52].
Q.ICD.Eps-11 was identified on chr. 4A and accounted 2 to 18 % of the phenotypic variation in PhEnv2 (Lebanon) in both germplasm types. This PhEnv differs from all others because both Vrn and Ppd requirements were presumably met. It is therefore likely that this locus is also acting as an earliness per se (Eps), with an effect visible only when Vrn and Ppd requirements are fulfilled. Kamran et al. [54] also detected an overlapping QTL (QFlt.dms-4A1), which induced early flowering in spring wheat population.
Q.ICD.Ppd-12 on chr. 4B controlled significant variation (3 to 30%) for both germplasm groups only in PhEnv4 (Lebanon off-season), and for modern lines only in the combined analysis across environments. PhEnv3 has similar climatic conditions to PhEnv4 (no vrn and no ppd requirements), but in PhEnv3 the photoperiod is decreasing, rather than remaining constant throughout the season. Hence, this unique locus might be specifically linked to the promotion of earliness in durum germplasm cultivated at high latitudes during the summer cycle. Giunta et al. [51], Sanna et al. [55], and Milner et al. [56] reported QTLs linked to flowering time that overlap with Q.ICD.Ppd-12 at 26.9 Mbp of chr 4B and suggested it might correspond to the effect of the Rht-B1 gene.
Q.ICD.Vrn-17 chr. 6B was effective in modern germplasm in all environments except Lebanon (PhEnv2), while for landraces it was identified only in the combined analysis. PhEnv2 is the only environment for which complete vernalization requirements can be met. As such, this locus is likely involved in controlling vernalization requirements in durum wheat with a more refined mode of action compared to Vrn1. Giunta et al. [51] and Würschum et al. [30] also identified QTLs on chr. 6B at 590.8 Mbp with a significant effect on the determination of heading time. Hence, given its stronger sensitivity to temperature changes, it might be of value to consider this locus in addition to Vrn1 to breed for earliness.

Define usable alleles for earliness via haplotype analysis of multiple QTLs
Four major QTLs (Q.ICD.Eps-03, Q.ICD.Eps-11, Q.ICD.Ppd-12 and Q.ICD.Vrn-17) were selected because identified in both germplasm types and with some consistency across PhEnv. Haplotype analysis was conducted for these QTLs to define their additive nature across climatic conditions ( Figure 4). In PhEnv1 the haplotype TATC resulted in significantly earlier modern lines than all other haplotypes, except CATC that matched. Comparison to the other haplotypes suggests a major role for the A allele of Q.ICD.Eps-11 to promote earliness. In PhEnv2 again TATC resulted in the earliest flowering haplotype, matching also CATC and CCTC, this time promoting the role of the TC combination for Q.ICD.Ppd-12 and Q.ICD.Vrn-17. In PhEnv3, TATC was the earliest flowering type, matching four other haplotypes. Only the nucleotide T in Q.ICD.Ppd-12 appears as shared among all haplotypes and hence playing a pivotal role. Finally, in PhEnv4 the haplotype TATC resulted in significantly earlier flowering genotypes than all other haplotypes, except CATC. Comparison to the other haplotypes suggests a major role for the A allele of Q.ICD.Eps-11 in promoting earliness. Considering these results together, the TATC haplotype appears as the best combination to promote earliness across at all tested conditions, and that good additive effect could be confirmed for at least three of the main QTLs (Q.ICD.Eps-11, Q.ICD.Ppd-12, and Q.ICD.Vrn-17), while it was not possible to discern the contribution of Q.ICD.Eps-03.
Across environments, two ICARDA's elites 'Icavicre' and 'Ouassara', having the TATC haplotype, were the earliest flowering overall. In addition to TATC haplotype, early flowering genotypes at individual PhEnv could be identified with different haplotypes, probably due to additional effect of other loci. For instance, 'CaMdoH25' and 'Icambel' in PhEnv1 with CCTC haplotype, 'Bradano' and 'Massara1' in PhEnv2 with CCTT and CATC, 'IDON37-039' and 'Moulsabil2' in PhEnv3 with CATC and CCTC, and 'IDON37-053' and 'Waha' in PhEnv4 with CCTC, were the earliest flowering entries at individual sites.

Conclusions
The present study took advantage of contrasting climatic conditions to assess the effect of different loci on the control of flowering time. Testing of known major flowering loci explained only partially the large variation observed. The use of GWAS allowed to define additional loci involved in the control of flowering for durum wheat, and the literature showed that other groups also identified the same. Hence, it appears that for durum wheat the flowering transition is controlled by more loci than common wheat, and that breeders will have to rely on these additional loci to promote earliness. In particular, haplotype analysis revealed an important additive role for three of the identified QTLs with a significant effect in controlling flowering time, beside the known major loci. The conversion of the markers underlying these loci in ready-to-use assays will promote their deployment by durum wheat breeders.

Plant material
A durum wheat core collection comprised of 96 landraces from 24 countries and 288 cultivars and elite breeding lines from eight countries, International Center for Agricultural Research in the Dry Areas (ICARDA), and International Maize and Wheat Improvement Center (CIMMYT) were used for this study. Detailed information regarding plant material is described in an earlier publication [57].  Figure 1. The experiments were conducted according to an augmented design with 19 blocks and four repeated checks. Days to heading (DTH) was recorded as the number of days elapsed from the date of sowing to the onset of flowering determined at 50% of the plot with the tip of the spike emerging from the flag leaf (Zadoks scale stage 51). Daily records at each environment were minimum and maximum temperatures, and length of daylight in minutes. To estimate the cumulative growing degree days (CGDD) needed for flowering, the average daily temperature from planting to flowering were summed for each site following procedures described by Klepper et al. [58]. In case of wheat, a range of 0 to 32°C temperature is considered optimal for growth, therefore, those values below and above these temperatures were converted to 0 and 32° C, respectively. Cumulative day length (CDL) was instead measured as the total sum of minutes of sunlight needed from planting to heading at each environment. In the cases of Terbol off-season and four environments along the Senegal River, the high temperatures and short photoperiod prevented few modern lines and several landraces from flowering before harvest time. Nevertheless, these results were deemed of great interest for the analysis so were converted to the maximum value for DTH, CDL, and CGDD recorded at each environment and used to run markertrait association studies.

Genotyping
All accessions were profiled by 35K Affymetrix Axiom wheat breeders' array (www.affymetrix.com) at Trait Genetics (Gatersleben, Germany) following the manufacturer instructions. In addition, three KASP assays were also run to characterize the variation for Ppd-A1, Ppd-B1, and Vrn-A1. A total of 10 sub-populations were identified based on genetic diversity as explained by Kabbaj et al. [57] and Sall et al. [59]. As described in earlier publication [3,60], 7,652 highfidelity polymorphic single nucleotide polymorphism (SNPs) were obtained, showing less than 1% missing data, minor allele frequency (MAF) higher than 5%, and heterozygosity less than 5%. The sequences of these markers were aligned with a cutoff of 98% identity to the durum wheat reference genome [52] (available at: http://www.interomics.eu/durum-wheat-genome), to reveal their physical position. A sub-set of 500 highly polymorphic SNPs were selected on the basis of even spread along the genome, and used to asses population sub-structure, which reveal the existence of 10 main subgroups [57]. To avoid discovery bias these 500 markers were then removed from all downstream analysis. Linkage disequilibrium was calculated as squared allele frequency correlations (r 2 ) in TASSEL V 5.0 software [61] using the Mbp position of the markers along the bread wheat reference genome and plotted using the "Neanderthal" method. The linkage disequilibrium decay was measured at 51.3 Mb as reported in Bassi et al. [60].

Statistical analysis
Grouping of environments was conducted on statistical software RX 64 (3.3.3) through stats package [62] via hierarchical clustering based on Euclidean distance calculated by principal component analysis (PCA), using as input the climatic variables measured in each environment: maximum temperature, average maximum temperature, minimum temperature, average minimum temperature, average cumulative growing degree days (CGDD), average days to heading (DTH), and average cumulative day length (CDL). The resulting clusters were defined as phenological environments (PhEnv). Combined analyses of variance were conducted across environments and PhEnv for DTH, CGDD, and CDL, assuming environment, PhEnv, and genotypes as fixed effects.
Each year x location combination was considered as one environment. Best linear unbiased estimates (BLUEs) were derived for the individual environment, PhEnv, and across environments based on a linear mixed model. All analyses were carried out with GENSTAT (version 2010) and free statistical package RX64 version 3.3.3.
The broad sense heritability was estimated as indicated by Falconer et al. [63]: Where σg2 is the genotypic variance and σp2 is the phenotypic variance. The genotypic and phenotypic variance components were estimated based on the method suggested by Burton and Devane [64]: Where MSg and MSge are the mean square due to genotype and GxE interaction, MSe is the error mean square, and r is the number of replicates.
The association analysis was performed with TASSEL version 5.2.38. The marker-trait association test was carried out using mixed linear model (MLM) based on the kinship matrix estimated by Kabbaj et al. [57]. The analysis was performed using the BLUE across environments and for each of the four PhEnv. In each case, landraces and modern lines were analyzed separately. The significance of marker-trait association (MTA) was assessed by Boneforroni corrected equation as suggested by Duggal et al. [65], assuming 288 marker trait hypothesis (12,000 Mbp of durum genome divided by LD decay of 51.3 Mb), which resulted in a LOD threshold of 3.0 and 3.4 for p<0.05 and p<0.01, respectively. Significant MTAs located at less than twice the LD decay distance were merged into one QTL. These QTLs were then further assessed by factorial regression to determine the true marker effect for each of the four PhEnv. Pearson's critical value [66] for correlation was squared to obtain a critical r 2 = 0.0225 for p<0.05, which was used as threshold to determine significant effect on phenotypic variation.