Mapping Agronomic and Quality Traits in Elite Durum Wheat Lines under Di ﬀ ering Water Regimes

: Final grain production and quality in durum wheat are a ﬀ ected by biotic and abiotic stresses. The association mapping (AM) approach is useful for dissecting the genetic control of quantitative traits, with the aim of increasing ﬁnal wheat production under stress conditions. In this study, we used AM analyses to detect quantitative trait loci (QTL) underlying agronomic and quality traits in a collection of 294 elite durum wheat lines from CIMMYT (International Maize and Wheat Improvement Center), grown under di ﬀ erent water regimes over four growing seasons. Thirty-seven signiﬁcant marker-trait associations (MTAs) were detected for sedimentation volume (SV) and thousand kernel weight (TKW), located on chromosomes 1B and 2A, respectively. The QTL loci found were then conﬁrmed with several AM analyses, which revealed 12 sedimentation index (SDS) MTAs and two additional loci for SV (4A) and yellow rust (1B). A candidate gene analysis of the identiﬁed genomic regions detected a cluster of 25 genes encoding blue copper proteins in chromosome 1B, with homoeologs in the two durum wheat subgenomes, and an ubiquinone biosynthesis O-methyltransferase gene. On chromosome 2A, several genes related to photosynthetic processes and metabolic pathways were found in proximity to the markers associated with TKW. These results are of potential use for subsequent application in marker-assisted durum wheat-breeding programs. A.R.B. analysed the genotypic and phenotypic data. R.M.-G. isolated the DNA. R.M.-G., A.R.B. and P.H. carried out the AM analyses. R.M.-G., A.R.B. and P.H. drafted the manuscript. All authors have and agreed to the published version of the manuscript.


Introduction
Wheat is one of the most widely grown crops worldwide (FAO, 2015), and is essential for the human diet [1]. Its importance and worldwide dominance are due, in part, to its agronomic adaptability. Durum wheat (Triticum durum) is a tetraploid wheat species (AABB genomes) mainly grown in the Mediterranean basin, in the Northern Plains (between the USA and Canada), in the arid areas of South Western USA and in Northern Mexico [2]. Durum wheat is well-adapted to a broad range of climatic conditions (including dry environments) and marginal soils, and has low water requirements [3,4]. Climatic conditions, as temperature and water availability, together with biotic stresses, can strongly affect durum wheat development and production [3][4][5][6]. Crop adaptation field seasons is included in Supplementary Materials Figure S1. The agronomic and quality assessment of the panels over seasons is summarised in Table 1. Table 1. Agronomic and quality assessment of wheat field trials. The number of lines, year, location and water regime applied is shown.

Wheat Panel
No Plant material for genetic analysis was harvested for each line at the 4th leaf stage (growth stage 14 on the Zadoks scale [51]) and immediately frozen in dry-ice. Samples were stored at −80 • C until DNA extraction. Approximately 100 mg of frozen tissue was used for DNA isolation with a DNeasy Plant Mini Kit from Qiagen, following the manufacturer's protocol. DNA sample quality and concentration were assessed using electrophoresis on a 0.8% agarose gel and the restriction enzyme Tru1I (ThermoFisher) was used to check for the absence of nucleases in DNA prior to genotyping.
Samples were genotyped by Diversity Arrays Technology Pty Ltd. (Montana St, University of Camberra, Bruce ACT 2617, Australia) (DArT) using DartSeq TM . A total of 35,509 polymorphic dominant DArT markers and 9142 biallelic SNP markers were generated. Both datasets were thinned by removing one marker from each pair with a correlation coefficient of >0. 95. The final dataset consisted on 14,588 DArT markers (of which 8411 were mapped) and 5716 SNP markers (4142 mapped markers). DartSeq TM genotyping and mapping of the corresponding markers to the wheat reference genome sequence RefSeq v1 from the International Wheat Genome Sequencing Consortium (IWGSC, http://www.wheatgenome.org/) was performed by DArT (diversityarrays.com), as described by Sukumaran et al. [43]. The distribution of markers across the A and B subgenomes is given in Table 2.

Phenotypic Data
Ten agronomic and quality traits were assessed for three durum wheat elite line panels: days to heading (days, DTH); plant height (cm, PH); lodging (%, LOD); yellow rust (%, YR); yellow colour (YC); sedimentation index (cm 3 , SDS); sedimentation volume (cm 3 , SV); grain protein content (%, GPC); thousand kernel weight (g, TKW); and grain yield (Kg/ha, GY). Agronomic traits (DTH, PH, LOD, YR, TKW and GY) were assessed under both water treatments (FI and RI), and quality traits (YC, SV, SDS and GPC) were only evaluated under FI conditions. Visual disease evaluation and phenology assessments were made in the field, while quality parameters were evaluated on grain samples post-harvest. DTH, PH and LOD and YR were visually assessed at the field trials in Yaqui, while YR was assessed at Toluca. To assess DTH, heading date was recorded as the time when 50% of the spikes have emerged from the flag leaf sheath (stage 59 in Zadoks scale [51]); PH was recorded by measuring the distance between the stem's base and the top of the spike (awns not included); LOD was assessed as the percentage of lodging plot; and YR was assessed as the percentage of leaves with rust pustules. YC was assessed by a rapid colorimetric method with a Minolta color meter following CEN/ TS 15,465:2008 [52-54]; SDS was evaluated following UNE 34,903:2014 [55,56]; SV and GPC were assessed by Near-infrared spectroscopy (NIRs) [57]; TKW was measured by weighing 2 samples of 100 entire kernels randomly chosen previously dried at 70 • C for 48 h.
The correlation between the assessed traits was analysed using the 'cor' function in R [58][59][60]. Then, an analysis of variance (ANOVA) was undertaken, using the 'aov' function in R [61], to obtain the descriptive statistics for each trait.
Traits were analysed using a Q + K linear mixed-model [62,63] which follows the model equation: where y is a vector of observed phenotypes; X, S and Z are matrices related to β, α and µ, respectively; β is a vector of fixed effects; α is a vector of marker effects; Qv is a vector of population effect; µ is a vector of polygenic effects (with covariance proportional to a kindship or relationship matrix); and ε is a vector of residuals. These analyses were carried out using GenStat (14th Edition) to generate the best linear unbiased estimates (BLUEs) of variety performance in different ways: (i) across years and blocks; (ii) across years for each block (FI and RI); (iii) across a reduced dataset (years 2013 and 2014) and blocks; and (iv) across the reduced dataset for each block. The resulting datasets (available in Supplementary Material Table S2) were then used in different association mapping analyses.

Population Structure and Linkage Disequilibrium
Population structure was assessed using principal component analysis (PCA) based on the combined DArT and SNP genotyping datasets. Euclidean distances were calculated using the R package 'ggfortify' [64] and the PCA was visualised with 'ggplot2' [65].
The pattern of linkage disequilibrium (LD) was assessed between each pair of SNP markers on the same chromosome across the two constitutive genomes with the allele frequency correlation (r 2 ) using the 'popgen' package in R [66]. A heatmap was obtained with the D' and r 2 values for each chromosome and a scatterplot to determine LD decay (genetic distance in cM).

Association Mapping (AM)
The AM analyses were performed on the BLUEs obtained above using an additive model with 'rrBLUP' [67] and 'GWASpoly' [68] packages in R in different ways. Two marker-based kinship matrices (k-matrix), created from a subset of 14,588 DArT and 5716 SNP markers, respectively, were used for the adjustment based on relatedness of individuals (Supplementary Materials Tables S3 and S4). A minor allele frequency (maf) threshold of 0.05 was used. To establish a p-value detection threshold for statistical significance of associations, the Bonferroni correction, which employs a threshold of α/m to ensure the genome-wide type error I of 0.05, was applied with a total of 1000 permutations.
Associated DartSeq TM and SNP markers were blasted against the wheat reference assembly RefSeqv1 [69] with no indels or mismatches allowed, using an ad hoc Java program, to confirm their physical mapping location on the A or B genomes. The molecular markers were also mapped against the durum wheat genome (https://www.interomics.eu/durum-wheat-genome) to confirm their physical positions. In addition, to identify candidate genes, results were filtered, selecting the hits with best e-value for each molecular marker, and candidate genes were manually selected based on their annotations.

Phenotypic Assessment
Results from the ANOVA for all the traits across years and water treatments are shown in Table 3. The mean phenotypic values across years were calculated for each block and panel to evaluate the influence of water conditions on the assessed traits (Table 3). Days to heading during the field seasons assessed ranged from 63 to 94 days. In plots with lower water availability (RI block), the spike emergence from the flag leaf took place approximately 11 days earlier than in FI plots. Plant height ranged from 39 to 110 cm showing differences between water regimes, with a decrease of 25-30 cm under RI conditions. Likewise, and as result of the RI treatment, GY (ranging from 1.35 to 10.63 ton/ha) and TKW (from 29.6 to 63.2 g) also varied, being reduced by 4-5 tons/ha and 7-10 g, respectively, in the low water availability RI treatment. This strong RI treatment resulted in very low heritability values for DTH, PH and LOD. Several significant phenotypic correlations were observed between the analysed traits (

Population Structure and Linkage Disequilibrium
The PCA used a total of 14,588 DArT and 5716 SNP markers. The first and second principal components explained 3.91% of the genetic variation ( Figure 2). No underlying genetic structure was detected within or between the panels assessed. LD was estimated using the mapped SNP markers dataset. LD decay was determined within 20-30 cM for all the chromosomes (Figure 3). Using the classification defined by Maccaferri et al. [70], the markers presented loose linkage (Class 2), showing a distance value between 21 to 50 cM.

Population Structure and Linkage Disequilibrium
The PCA used a total of 14,588 DArT and 5716 SNP markers. The first and second principal components explained 3.91% of the genetic variation ( Figure 2). No underlying genetic structure was detected within or between the panels assessed. LD was estimated using the mapped SNP markers dataset. LD decay was determined within 20-30 cM for all the chromosomes (Figure 3). Using the classification defined by Maccaferri et al. [70], the markers presented loose linkage (Class 2), showing a distance value between 21 to 50 cM.

AM Analysis
Thirty-seven significant marker-trait associations (MTAs) were detected for TKW and SV across all years and water treatments with most of the significant markers located on chromosome 2A (Table 4). Twenty DArT and seven SNP markers were found in association with TKW on chromosome 2A (with additive effects ranging from −3.41 to 3.46). In addition, eight unmapped DArT and one SNP marker were also associated with TKW (additive effects ranged from −3.39 to 3.46 g). Most of these MTAs showed a negative additive effect, reducing the final weight value (ranging from −2.84 to −3.19 g), and only two MTAs were found to increase TKW (values of 2.97 and 3.09 g). Finally, a single SNP associated with SV was located on chromosome 1B (showing a positive effect increasing the final value by 1.26 mL). The resulting Manhattan and QQ-plots from this AM analysis are included in The AM analyses on partitioned subsets of the data consistently detected the QTLs for TKW and SV. Nevertheless, the individual assessment of the water treatments significantly reduced the number of MTAs found, due in part to less available data for the RI block (Supplementary Material Table S6). The initial dataset of 294 durum wheat elite lines was reduced to 200 lines (assessed during the 2013 and 2014 seasons) to give a dataset balanced across assessment years. Using this reduced dataset for AM analysis, the results confirmed the QTLs previously found for the full dataset (Supplementary Materials Table S6). The analysis also detected an additional locus for SV on chromosome 4A, and a locus for YR on chromosome 1B (with additive effects of −0.84 and 2.79, respectively).

Candidate Genes Analysis
The marker SNP620, located on chromosome 1B and detected in association with SV, was found included into a cluster of 12 genes encoding blue copper proteins (BCP), with homoeologs in the two durum wheat subgenomes ( Figure 4, Table 5 and Supplementary Material Figure S6). In the hexaploid wheat genome, this set of genes form a cluster of homeolog triads [72] with a total of 31 genes (Supplementary Materials Table S7 and Figure S7). Additionally, another interesting gene (TraesCS1B01G568400LC.1) was found closer this marker, coding for the ubiquinone biosynthesis O-methyltransferase.
There were several markers located in chromosome 2A, in close proximity to some interesting genes. Markers SNP1183, SNP1184 and DArT3165 were found next to several genes encoding reductase-1 ( Figure 4 and Table 5). In addition, the marker SNP8395, also located on the same chromosome, was found in proximity to the gene TraesCS2A01G309700.1, which encodes a type A response regulator 1 ( Figure 4 and Table 5).
Significant MTAs from the partitioned analysis also allowed the identification of potentially interesting genes. On chromosome 1B, marker DArT1744, previously described by Mérida-García et al. [73] related to high molecular-weight glutenin subunits, was found in proximity to genes encoding isocitrate dehydrogenase kinase/phosphatase G and leucine-rich repeat receptor-like protein kinase family protein. These genes participate in the carbohydrate metabolism during the Krebs cycle and play a crucial role in plant development and stress responses, respectively [74,75]. On this chromosome, another marker (SNP809) was found in proximity to some interesting genes encoding sugar transporter proteins. Additionally, some markers located on chromosome 2A were found in proximity to Acyl-CoA N-acyltransferase genes (SNP1206, SNP8395 and DArT3180) and chloroplastic zeaxanthin epoxidase (SNP1189) (Supplementary Materials Table S8).

Candidate Genes Analysis
The marker SNP620, located on chromosome 1B and detected in association with SV, was found included into a cluster of 12 genes encoding blue copper proteins (BCP), with homoeologs in the two durum wheat subgenomes (Figure 4, Table 5 and Supplementary Material Figure S6). In the hexaploid wheat genome, this set of genes form a cluster of homeolog triads [72] with a total of 31 genes (Supplementary Materials Table S7 and Figure S7). Additionally, another interesting gene (TraesCS1B01G568400LC.1) was found closer this marker, coding for the ubiquinone biosynthesis Omethyltransferase.
There were several markers located in chromosome 2A, in close proximity to some interesting genes. Markers SNP1183, SNP1184 and DArT3165 were found next to several genes encoding reductase-1 ( Figure 4 and Table 5). In addition, the marker SNP8395, also located on the same chromosome, was found in proximity to the gene TraesCS2A01G309700.1, which encodes a type A response regulator 1 (Figure 4 and Table 5). Significant MTAs from the partitioned analysis also allowed the identification of potentially interesting genes. On chromosome 1B, marker DArT1744, previously described by Mérida-García et al. [73] related to high molecular-weight glutenin subunits, was found in proximity to genes encoding isocitrate dehydrogenase kinase/phosphatase G and leucine-rich repeat receptor-like protein kinase family protein. These genes participate in the carbohydrate metabolism during the Krebs cycle and play a crucial role in plant development and stress responses, respectively [74,75]. On this chromosome, another marker (SNP809) was found in proximity to some interesting genes encoding sugar transporter proteins. Additionally, some markers located on chromosome 2A were found in proximity to Acyl-CoA N-acyltransferase genes (SNP1206, SNP8395 and DArT3180) and  Table 5. Candidate genes for markers with significant marker-traits associations. Genes located in proximity of markers found in association with TKW and SV across years and water treatments (within a ±50 kb window). Values for physical position and distance are indicated in base pairs (bp). Chr: chromosome position. Blue copper proteins are shown in blue colour; ubiquinone biosynthesis O-methyltransferase in purple colour; the regulator response gene in brown colour; and for reductase 1 genes in green colour. Physical positions and gene annotations are based on the wheat reference assembly RefSeqv1 [69].

Discussion
The maintenance of crop production is a current and pressing need given growing populations and reduced availability of arable land [76]. There is an increasing need for breeding programs to improve yield potential and the adaptation of new varieties to different biotic and abiotic stresses [77]. Abiotic stresses, including drought and heat, are impacting productivity on the million hectares of wheat grown worldwide each year [78]. Detailed molecular and phenotypic characterization are valuable tools in the dissection of complex traits [79], and especially those that are influenced by water availability [14].
The improvement of key traits is essential for better end-use production quantity and quality in wheat [80]. In this study, we analysed a set of 10 agronomic and quality traits under full irrigation conditions (FI), with an additional six traits also assessed under low water availability (RI) in order to understand trait variation under contrasting water regimes in the CIMMYT durum wheat breeding programme. Irrigation conditions influenced some important yield and yield-related traits such as GY and TKW, as well as adaptive traits including DTH and PH ( Table 3). The RI treatments had decreased final yields in line with previous observations [6,81]. Previous reports have also shown TKW to be reduced by high temperatures [17], most likely related to water availability. Groos et al. [82] assessed the genetic basis of grain yield and protein content in a segregating population of wheat RILs grown at six locations and also identified effects from GxE interactions involving protein content and yield. Our mean trait values corroborated this, with the highest values for GY recorded for FI blocks across panels (see Table 3). A similar trend was shown for DTH, PH and TKW, which decreased under low-water regimes.
Correlations between the assessed traits showed that GY was positively correlated with two different phenology traits (PH and DTH). This is in agreement with Maccaferri et al. [6], who showed important positive and negative correlations for GY and DTH, and also positive correlations for GY and PH in several environments with different water regimes. DTH and PH were also positively correlated (Supplementary Materials Table S5), with taller plants having a longer time period to the emergence of the tip of the spike stage.
Wheat TKW is an important yield component with a direct effect on grain yield [83,84]. In line with this, our results showed a significant and positive correlation between TKW and GY. However, the previously reported negative correlation between TKW and DTH [6] was not observed, potentially as result of temperatures and water availability from emergence to heading, and also from heading to harvest. Rharrabti et al. [17] previously highlighted a positive correlation between protein content and TKW, which is in agreement with the results obtained in the present study. They highlighted that warm temperatures during grain filling could negatively affect this correlation.
Significant associations between endosperm proteins (gliadin and glutenin subunits) and SV have been previously highlighted [85,86]. Here we found a positive correlation (r = 0.15) between SV and GPC. This correlation is thought to be the result of grain protein content influencing the sedimentation volume value [87], which depends on the degree of protein hydration and oxidation [88]. Finally, for sedimentation index (SDS) analysis, we observed a negative correlation with protein content, in agreement with results presented by Rharrabti et al. [17,45]. This is also in agreement with Oelofse et al. [89] who highlighted the significant influence of protein content on the SDS sedimentation test [90][91][92].
The SNP and DArT markers used to analyze patterns of genetic structure ( Figure 2) and LD (Figure 3) for the durum wheat lines revealed no detectable genetic structure and consistent patterns of LD across chromosomes (LD was estimated to extend~20 cM). These results support the suitability of durum elite line sets currently in use in breeding programmes for the practical application of GWAs analysis. The rate of unmapped markers was lower for SNP than for DArT markers (27.5% vs. 42.0%, Table 2), indicating higher precision in genetically mapping SNP markers, probably as a result of co-dominance.
In the assessment of MTAs for quality and yield-related traits, different AM analyses were performed on subsets of the dataset. Several MTAs for SV and TKW were detected across years and water regimes, located on chromosomes 1B and 2A, respectively. All GWAS analyses corroborated the major QTLs previously detected, and reported two new QTLs, one for YR in chromosome 1B, and another for SV in chromosome 4A.
Associations on chromosome 1B were significant for wheat quality. There are known loci including Gli-B1/Glu-B3 on this chromosome, which are of great importance for some gliadin and glutenin subunits [93]. In fact, the candidate gene analysis reported the presence of a high molecular-weight glutenin subunit (HMW-GS) gene in the proximity of marker DArT1744 (found to be significantly associated with SV and SDS), which was previously described by Mérida-García et al. [73] in association with specific weight. In line with this, Pogna et al. [93] highlighted the importance of Glu-B3 for determining protein quality with these endosperm proteins showing significant effects on SV, which also showed a high and positive correlation with SDS in our study (Figure 1 and Supplementary Materials Table S5). Likewise, Blanco et al. [86] reported three QTLs on chromosomes 1B, 6A and 7B (based on the analysis of 259 polymorphic markers) associated with SDS and SV in a recombinant inbred line population. In the present study, we found a SNP marker (SNP620) associated with SV, showing a positive additive effect of 1.26 (see Table 4) and also with SDS (marker effect of 0.11) (Supplementary Materials Table S6). This marker was previously placed on chromosome 1B, in the same location as MTAs for gluten index and sedimentation index [73]. Other previous studies, such as Reif et al. [94] and Fiedler et al. [95], also reported markers associated with SV on chromosome 1B, but with differing genetic positions. The additional locus for SV found on chromosome 4A (marker DArT9459) has not been previously reported.
Marker SNP620, significantly associated with SV, is located within a cluster of homoeolog gene triads coding for blue copper proteins (Table 5 and Supplementary Materials Figure S6). These proteins have been described containing a copper atom, and participate in redox processes [96], with a crucial role in the electron shuttle in plants [97]. In addition, Yao et al. [98] described the blue copper protein genes as the targets of miR408 in wheat, which is involved in the regulation of gene transcription required for heading time [99]. In our study SNP620 was also found in proximity to a gene coding for an ubiquinone biosynthesis O-methyltransferase. Liu et al. [100] highlighted its crucial role as an electron transporter in the electron transport chain of the aerobic respiratory chain. This ubiquinone gene is involved in plant growth and development, is implied in chemical compounds biosynthesis and metabolism which are involved in plant responses to stress, and also participates in gene expression regulation and cell signal transduction [100].
On chromosome 1B we also found a significant MTA for yellow rust, in agreement with previous studies in durum and bread wheat, which placed different markers significantly associated with this trait, but in differing genetic positions [101][102][103][104]. The candidate gene analysis revealed the proximity of this marker (SNP809) to sugar transporter protein genes. Sugars are formed during the photosynthetic process and are essential for plant nutrition. The sucrose transport has been considered of great importance for plant productivity [105]. In line with this, the sucrose is involved in the gene expression regulation of the supposedly sugar-sensing pathway [106,107].
The majority of MTAs for TKW were located on chromosome 2A, showing both positive and negative effects. Previous studies have reported different markers in association with this quality trait, including a number mapped on chromosome 2A [38,[108][109][110]. One of the markers found by Yao et al. [38] (SSR marker gwm445 on chromosome 2A (68.2 cM), belongs to the same QTL found in this study for the marker SNP1153 (chromosome 2A, 68.6 cM), and also found by Juliana et al. [111] in bread wheat lines from CIMMYT's first year-yield trials. Sukumaran et al. [43] analysed a durum wheat panel of 208 lines under yield potential, heat and drought stress conditions, and identified markers on chromosome 2A with a similar position to those detected in this study (4 markers at 70 cM and 6 markers at 69 cM) under heat stress conditions. They highlighted that several SNP markers were related to transmembranes or were uncharacterized proteins. We found several candidate genes for this important TKW QTL (Table 5) among which the most striking feature is the presence of four reductase 1 genes (NADPH-dependent 6 -deoxychalcone synthase) and a type A response regulator 1 (Figure 4). These genes are both related with photosynthesis. Hu et al. [112] highlighted that NADPH plays a crucial role in biological processes in plants, such as the regulation of the production of reactive oxygen species (ROS) for the stress tolerance [113,114]. Additional GWAS analyses using reduced datasets revealed other interesting genes for this QTL (chromosome 2A, Supplementary Materials Table S8), encoding for the Acyl-CoA N-acyltransferase and the chloroplastic zeaxanthin epoxidase. The first gene has several functions in signaling and metabolic pathways [115]. The zeaxanthin epoxidase plays an important role in the xanthophyll cycle and abscisic acid (ABA) biosynthesis. The xanthophyll cycle has a main function in the dissipation of light energy excess and also increasing the photosynthetic system stability [116].
The proposed approach has successfully detected genetic markers with significant associations with TKW, SV, SDS and YR. These are of potential use in durum wheat breeding programs, and can be further interrogated to the candidate gene level using the RefSeqv1 bread wheat genome reference [69] and the durum wheat genome reference [71].  Figure S4: Quantile quantile-plots from genome-wide association studies (GWAS) analysis for durum wheat DArT markers (mapped and unmapped). DTH: days to heading; PH: plant height; LOD: lodging; GY: grain yield; TKW: thousand kernel weight; YC: yellow color; SV: sedimentation volume; SDS: sedimentation index; and GPC: grain protein content; Figure S5: Quantile quantile-plots from GWAS analysis for durum wheat SNP markers (mapped and unmapped). DTH: days to heading; PH: plant height; LOD: lodging; GY: grain yield; TKW: thousand kernel weight; YC: yellow color; SV: sedimentation volume; SDS: sedimentation index; and GPC: grain protein content; Figure S6: Blue copper protein gene cluster on durum wheat chromosome 1B. High confidence genes are shown in green colour, low confidence genes are shown in yellow; Figure S7: Cluster tree of blue copper protein gene homoeologs in bread wheat (RefSeqv1 [69]). For chromosome 1A, high confidence (HC) and low confidence (LC) genes are shown in brown and orange colour, respectively; for chromosome 1B, HC and LC genes are shown in dark and light green colour, respectively; for chromosome 1D, HC and LC genes are shown in dark and light blue colour, respectively; Table S1: Durum wheat elite lines assessed; Table S2: Best Linear Unbiased Estimates (BLUEs) outputs for all assessed traits and the association mapping analyses performed in durum wheat: [i] across years and blocks; [ii] across years for each block (FI and RI); [iii] across years and blocks for a reduced dataset (years 2013 and 2014); and [iv] across the reduced dataset for each block. DTH: days to heading; GPC: grain protein content; GY: grain yield; PH: plant height; SDS: sedimentation index; SV: sedimentation volume; TKW: thousand kernel weight; and YC: yellow colour; YR: yellow rust; LOD: lodging; Table S3: Kinship matrix for durum wheat DArT markers; Table S4: Kinship matrix for durum wheat SNP markers; Table S5: Phenotypic correlations between the assessed traits in durum wheat and their corresponding p values. YR: yellow rust; DTH: days to heading; PH: plant height; LOD: lodging; GY: grain yield; TKW: thousand kernel weight; YC: yellow color; SV: sedimentation volume; SDS: sedimentation index; and GPC: grain protein content; Table S6: Marker-trait associations found for the association mapping analyses performed in durum wheat: [i] across years and blocks; [ii] across years for each block (FI and RI); [iii] across years and blocks for a reduced dataset (years 2013 and 2014); and [iv] across the reduced dataset for each block. SV: sedimentation volume; TKW: thousand kernel weight; SDS: sedimentation index; YR: yellow rust; "-": unmapped marker; Table S7: Homoeolog triads for blue copper protein genes mapped in the wheat reference assembly RefSeqv1 [69]; Table S8: Candidate genes for GWAS analyses performed in durum wheat: [i] across years for each block (FI and RI); [ii] across years and blocks for a reduced dataset (years 2013 and 2014); and [iii] across the reduced dataset for each block. Molecular markers mapping positions are shown both in the durum wheat genome [71] and the wheat reference assembly RefSeqv1 [69]; Supplementary Material S1. R script used to perform the GWAS analysis; Supplementary Material S2. R script used to perform the LD analysis. Funding: This research was funded by project P12-AGR-0482 from Junta de Andalucía (Andalusian Regional Government), Spain (Co-funded by FEDER). ARB is supported by the UK Biotechnology and Biological Sciences Research Council (BBSRC) 'Designing Future Wheat' cross-institute strategic programme (BB/P016855/1). PH is supported by the Spanish Ministry of Economy, Industry and Competitiveness (MINECO) project AGL2016-77149-C2-1-P.

Conflicts of Interest:
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.