followed the same agronomic practices, with the difference that at the time of booting (Z45) a 10 m2 and 1.5 m high plastic tunnel was placed over each block (Figure 1) and left there until early dough stage (Z83). An electronic thermometer (temperature data logger) was placed in the middle of each block (normal and heat stressed) to reveal that the temperatures were up to 16◦ C higher inside the plastic tunnels, to reach a maximum of 49 ◦C (Figure 1). Marchouch is a drought prone site, and no rainfall occurred after Z45 in any of the two field seasons.

**Figure 1.** Mean temperature difference of 18 days over two seasons between the plastic tunnel-mediated heat stress and normal field conditions between 8 a.m. and 8 p.m., and a picture of the plastic tunnel at 9 a.m.

The following traits were recorded: days to heading (DTH) measured at the moment when the awns became visible, plant height (PH) measured from the ground to the top of the highest spike excluding the awns, and the number of fertile spikes per meter square (Spkm2) was counted in a 0.25 m2 area. The whole plot was harvested by hand and the dry biomass (Biom) was weighed before threshing. Grain yield (GY) was weighed for each plot and expressed as kg ha−1. The weight of a thousand kernels (TKW) was expressed in grams. The harvest index (HI) was calculated as the ratio between GY and Biom. The grain number per spike (GNSpk) was derived from dividing grain number per meter square by Spkm2 as follows:

$$\text{Grain number} / \text{m}^2 = \frac{\text{Grain weight of the plot}}{1.5 \text{m}^2 \times \frac{\text{TKW}}{1000}} \tag{1}$$

$$\text{GNSpk} = \frac{\text{Grain number} / \text{m}^2}{\text{Spkm}^2} \tag{2}$$

The second and third sets were field tested in Kaedi, Mauritania (16◦14" N; 13◦46" W) during season 2014–2015 and 2015–2016 where the temperature reached a maximum of 41 ◦C and an average maximum daily temperature of 34 ◦C throughout the season. The trial was carried out under augmented design with a plot surface of 4.5 m2. Standard agronomic management practices were adopted. Full details for this experiment are published elsewhere [35].

#### *2.3. Data Analysis*

A mixed linear model was run using the lme4 package [36] in R [37] to obtain best linear unbiased estimates (BLUEs) of the normally distributed traits. For count traits (DTH, Spkm2, GNSpk), the generalized mixed linear model was used to get the BLUEs by Proc GLIMMIX in SAS. In both models, genotype, treatment, year, and replication were considered as fixed effects and block as random effect

nested in treatment and year. Broad-sense heritability was calculated based on variance components from random model using the method suggested by DeLacy et al. [38]:

$$\text{H}^2 = \frac{\sigma^2 \text{g}}{\sigma^2 \text{g} + \frac{\sigma^2 \text{G} \text{x} \text{T}}{t} + \frac{\sigma^2 \text{G} \text{x} \text{Y}}{y} + \frac{\sigma^2 \text{G} \text{x} \text{T} \text{x} \text{Y}}{ty} + \frac{\sigma^2 \text{c}}{ty \text{r}}} \tag{3}$$

where: σ2 *G*×*T* = genotype × treatment variance, σ2 *G*×*Y* = genotype × year variance, σ2 *G*×*Y*×*T* = genotype × treatment × year variance, σ2 *e* = residual variance, *r* is the number of replications per treatment, *t* is the number of treatments, and *y* is the number of years.

Box-and-whisker plots where constructed by ggplot2 package [39] using the BLUEs combined over year per each treatment. The relationship between the target trait GY and yield components (GNSpk, TKW, Biom, HI) was studied using the Pearson correlation coefficient and the additive regression model. The critical value of the correlation significance was determined at 0.30 for *p* < 0.05 and 0.39 for *p* < 0.01 for 40 df using the corrplot package [40]. The additive model incorporates flexible forms (i.e., splines) of the functions to account for non-linear relationship contrary to linear regression model estimated via ordinary least squares [41]. For the additive model, the effective degree of freedom term determines the nature of the relationship between the predictor and the response variables where EDF = 1 indicates linearity and EDF > 1 the non-linearity. The additive regression analysis was performed using the mgcv package [42].

Two stress tolerance indices were calculated to identify the heat tolerant genotypes. The stress susceptibility index (SSI) [43,44] was calculated as follows:

$$\text{SSI} = \frac{\left[1 - (\text{Ys}) / (\text{Yp})\right]}{\left[1 - (\text{Ys}) / (\text{Yp})\right]} \tag{4}$$

where Ys and Yp are yield values of the genotypes evaluated under heat stress and normal conditions, respectively, and Ys and ¯ Yp are the mean yields of the lines evaluated under heat stress and normal ¯ conditions, respectively.

The stress tolerance (TOL) [45] was calculated as follows:

$$\text{TOL} = \text{Yp} - \text{Ys.} \tag{5}$$

The classInt package [46] was used to identify the possible number of class intervals of the indices for the frequency distribution of the subset.

The cut-off value for tolerant vs. susceptible genotypes for SSI was equal to 1, with lines having SSI < 1 being stress tolerant. Regarding the TOL index, the smaller TOL values indicate the genotypes with low yield depression and hence more tolerant. The experiment-wide TOL mean (1608 kg ha−1) was identified as the cut-off value for tolerant vs. susceptible. The emmeans package [47] based on ANOVA model was used to discriminate among the grain yield means of haplotypes.

#### *2.4. Genotyping and Marker-Trait Associations*

Details of the genotyping step of the core set and panel have been previously discussed in Kabbaj et al. [34] and Sall et al. [35]. Briefly, 7652 high-fidelity 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 cut-off of 98% identity to the durum wheat reference genome [48] (available at: http://www.interomics.eu/durum-wheat-genome), to reveal their physical position. The average length of the Axiom probe is of 75 bp, hence the 2% allowed miss-match was set to account for the existence of 1 SNP within each sequence. A sub-set of 500 highly polymorphic SNPs were selected on the basis of even spread along the genome, and used to identify the existence of population sub-structure, which revealed the existence of 10 main sub-groups [34]. To avoid bias, these 500 markers were then removed

from all downstream association analysis. Linkage disequilibrium was calculated as squared allele frequency correlations (*r*2) in TASSEL V 5.0 software [49], using the Mb position of the markers along the bread wheat reference genome. Linkage disequilibrium (LD) decay was estimated and plotted using the "Neanderthal" method [50]. The LD decay was measured at 51.3 Mb for *r*2 < 0.2 as presented in Bassi et al. [51].

The genome wide association study (GWAS) was based on BLUEs of all the traits that displayed a significant treatment effect and the two stress tolerance indices. Two models were fitted and compared using two covariate parameters, Q (population structure) and K (Kinship). Q model was performed using a general linear model (GLM), and Q + K model using a mixed linear model (MLM). The best model for each trait was selected based on the quantile-quantile (Q-Q) plots [52]. Flowering time (DTH) was used as covariate in all analyses to remove the strong effects of flowering genes from the study. The value calculated for the LD decay of 51.3 Mb indicated that this association panel interrogated the 12,000 Mb of the durum wheat genome via 248 "loci hypothesis," and hence the Bonferroni correction for this panel was set to 3.1 LOD for *p* < 0.05 as suggested by Duggal et al. [53]. Local LD decay for *r*2 < 0.2 was calculated for a 100 Mbp window around the marker with highest LOD for all marker-trait associations (MTAs) identified at a distance inferior to 104 Mbp (twice the LD decay). The MTAs that occurred at a distance inferior to twice the local LD were considered to belong to the same QTL. QTL associated to flowering time were removed from all downstream analyses (Table S2). A regression analysis was performed between the haplotype of the peak marker of each QTL to determine possible duplicate or homeolog loci. In addition, all the MTAs analyses were performed using Tassel 5 software [49].

#### *2.5. Markers Conversion to KASP (Kompetitive Allele Specific PCR)*

The array sequences of 20 markers associated to traits (MTA) were submitted to LGC Genomics for in-silico design of KASP primers using their proprietary software. Those that passed the in-silico criteria were purchased and used to genotype the independent validation set. For each marker that amplified and showed polymorphism, the regression cut-off between phenotype and haplotype was imposed at *r* = 0.105 following Pearson's critical value [54]. KASP markers AX-95260810, AX-94432276, and AX-95182463 were tested for association with grain yield, while AX-94408589 for association with biomass. In addition, the top 20 and worst 20 lines were considered as the true positive and true negative for heat tolerance. Hence, the accuracy was calculated as the ratio of the correct allelic call among all, sensitivity as the ratio of the correct positive allelic among the top 20 yielding lines, and specificity as the ratio of the correct negative (wt) allelic calls among the 20 worst yielding lines. The sequence of the validated KASP markers is provided in Table S3, or the primers can be ordered directly at LGC Genomics indicating the Axiom code used in this article.

#### **3. Results**

#### *3.1. Agronomic Performance of the Genotypes and Sensitivity of Traits to Heat Stress*

The combined analysis of variance across four environments (two different temperature treatments over two crop seasons) revealed significant genotypic differences for all traits measured (Table 1). The yield performance of the genotypes across environments averaged 2171 kg ha−1 and ranged from 352 kg ha−1 obtained under heat stress conditions for the lowest yielding line DWAyT-0215, to 4658 kg ha−1 under normal conditions for the highest yielding line DWAyT-0217. The top yielding line under heat-stress was the ICARDA/Moroccan cultivar 'Faraj' with an average yield of 2249 kg ha−1 over the two seasons.

**Table 1.** Descriptive statistics, component of trait variation, and heritability (h2) among a set of 42 durum genotypes (G) tested under two treatments (T): normal and plastic tunnel-mediated heat stress during seasons 2015–2016 and 2016–2017.

\*, \*\* Significant at the 0.05 and 0.01 probability levels, respectively.

The treatment effect was significant only for Biom, GY, HI, and GNSpk, whereas DTH, PH, Spkm2, and TKW were not significantly affected by treatments (Figure 2). The yield components were all significantly reduced under heat stress except TKW that showed a slight increase for the genotypes exposed to heat. The genotypes tested under plastic-tunnels had 61%, 54%, 42%, and 17% lower average GNSpk, GY, HI and Biom, respectively, compared to control. Relatively high heritability was observed for all the phenological and agronomical traits except for HI that had the lowest heritability (h2 = 0.20).

**Figure 2.** Boxplot of the best linear unbiased estimates (BLUEs) for various traits under two different environmental conditions (Heat: plastic tunnel-mediated heat stress and Normal) across two years. \*\* indicate significant difference between the means of control and heat-stressed plants at *p* < 0.05.

#### *3.2. The Traits Interrelationship under Each Environmental Condition*

Correlation analysis (Figure 3; Tables S4 and S5) was first conducted to investigate the interrelationship among all agronomic traits. Under both treatments, GNSpk had the highest association with GY (*r* = 0.81 under heat, *r* = 0.67 under normal), while Spkm2 and TKW were the least correlated with GY. Biomass was also correlated with GY with *r* = 0.61 under heat and *r* = 0.67 under normal conditions. HI also showed a significant positive correlation with yield under both treatments, but its effect was stronger under heat stress (*r* = 0.72) than normal conditions (*r* = 0.54). DTH was not significantly correlated to any trait except HI (*r* = −0.44) under normal conditions.

**Figure 3.** Relationships between grain yield (GY) and yield components (grain number per spike (GNspk), harvest index (HI), dry biomass (Biom), number of fertile spikes per meter square (Spkm2), weight of a thousand kernels (TKW)) under plastic tunnel-mediated heat stress and normal conditions assessed by Pearson correlation and simple generalized additive model. The continuous grey line represents a linear relationship; the dashed grey line represents a non-linear relationship. The thickness of the line indicates the level of predictivity of the trait for GY. The length of the lines represents the correlation, the shorter the line the more the trait is correlated to GY.

Among yield components, the only significant and positive associations under the two environmental conditions were observed between Spkm2, TKW, and Biom and between HI and GNSpk. Under heat conditions, a positive and significant correlation was noticed between GNSpk and Biom while under normal conditions HI was positively associated to TKW (Figure 3; Table S4).

The additive model was then used to further determine the nature of the relationship between GY and each predictor variable under normal and heat conditions (Figure 3; Table S5). The similarities observed between the two treatments in terms of the nature of relationship between GY and each of the predictors were the constantly linear and non-linear relationship between Spkm2, TKW and the response variable GY, respectively.

GNSpk was considered the best predictor (deviance = 0.73%) with a complex relationship (EDF = 2.64) with GY under heat stress, whereas under normal conditions this trait was the second best predictor (deviance = 0.44%) with a linear relationship (EDF = 1). A similar trend was observed for HI in both treatments. Biom was found to be the best predictor (deviance = 0.52%) for GY with a non-linear relationship (EDF = 2.52) under normal conditions (Table S2; Figure S1).

#### *3.3. Stress Tolerance Indices*

Two different stress tolerance indices were calculated for GY: SSI and TOL (Figure 4). The genotypes showed wide variation for these indices. Seven SSI groups were identified with four having an SSI lower than 1 and the three remaining groups of genotypes having SSI > 1. The frequency distribution of the panel showed a wide variation and indicated the presence of susceptibility, with 45% of the genotypes falling in the very heat-susceptible class of SSI higher than 1, and only 7% of the lines showing high tolerance at SSI < 1. For TOL index, seven groups were also identified with 48% of the lines showing high yield depression and 5% of the genotypes presenting high stability. The smaller TOL values indicate the genotypes with low yield depression and hence more tolerant. However, good heat tolerance can also be reached by low yielding lines, but their value for breeding would be questionable. Hence, a scatterplot was devised to compare the GY under normal conditions and each of the heat indices (SSI and TOL). Five genotypes (four ICARDA lines, one Moroccan cultivar): Kunmiki, Berghouata1, Margherita2, IDON37-141, and Ourgh were found to have above average yield, low yield depression (low TOL values) and good heat tolerance (SSI < 1).

**Figure 4.** Two different stress tolerance indices SSI (stress susceptibility index) and TOL (tolerance index) of grain yield, comparing plastic tunnel-mediated heat stress with normal conditions for the 42 durum wheat genotypes. The bars plot shows the frequency distribution of SSI and TOL for the genotypes tested. The dashed red lines mark the separation between tolerant (left) and susceptible (right) genotypes. The scatter plot shows the yield performance of genotypes tested under normal conditions against each of SSI and TOL. The vertical dashed red lines indicate the average GY. The horizontal dashed red lines indicate the cut-off value for tolerant vs. susceptible genotypes for each index. Red dots indicate genotypes that were identified as superior by both bi-plots.

#### *3.4. Markers Associated to Heat Stress Tolerance*

A total of 204 MTAs were identified for four traits (GY, GNSpk, HI and Biom) under both stress and normal conditions and 49 MTAs were recorded for the two GY stress tolerance indices. Regression analysis and clustering based on local LD decay confirmed that these associations were distributed over 12 loci (Table 2 and Table S6). Chromosome 1A had the highest number of MTAs (27) while chromosome 4A had the lowest (6).

Under normal conditions, 56 MTAs were detected for three traits GY, GNSpk, and HI, with the third trait having the highest number of MTAs (48). No common region for these traits was identified under the non-stress environment. Under heat stress, a higher number of associations (148) were identified with trait variation (*r*2) ranging from 0.25 to 0.36. The highest number of MTAs were detected for GNSpk distributed over 10 different loci, followed by HI on six loci. A common region for GY, GNSpk, HI, and Biom was identified under the heat condition on chromosome 6BS. Loci associated with both GNSpk and HI were detected on 1AL, 1BL, 2AL, 3AL, and 3BL. For heat tolerance indices (SSI-GY and TOL-GY), 49 MTAs were identified. The common loci associated with the two indices were on chromosomes 2AL, 5AL, and 5BL, while the loci on chromosomes 1AL and 6BS were identified only for TOL-GY and SSI-GY, respectively.

A comparison of the significant loci under each treatment and including the heat tolerance indices indicated a locus on chromosome 2AL, which was consistently identified for the indices, and both treatments for GNSpk and HI. Two loci on chromosomes 3AL and 3BL were associated with GNSpk and HI under both control and stress conditions, but were not associated with any of the indices. Three significant loci on chromosomes 1AL, 5BL, and 6BS were shared among heat stress treatment and stress tolerance indices, but not under normal conditions, making of these the most interesting genomic regions that specifically respond to heat stress. Overall, a total of 12 unique significant loci were identified (numbered QTL.ICD.Heat.01–QTL.ICD.Heat.12) and can be consulted in Table 2. Local LD decay was estimated for the 100 Mbp genomic region surrounding the peak marker. It varied between 31.7 and 108.7 Mbp, or a −38% to 112% variation compared to the average LD decay calculated for the whole panel (51.3 Mbp). This variation was accounted for to determine the correct physical size in each genomic region to assign multiple MTAs to the same QTL.

have been converted into KASP markers and validated; GNspk—Grain number per spike; HI—Harvest index; TOL-GY—Tolerance index for grain yield; SSI-GY—Stress susceptibilityindexforgrainyield;GY—Grainyield;Biom—Biomass.

#### *3.5. E*ff*ect of Di*ff*erent Allele Combination on Yield Performance*

The loci identified on chromosomes 1AL, 5BL, and 6BS appeared as the most critical for heat tolerance and were then tested further. These regions were associated with the control of multiple traits under heat stress: GY, GNspk, HI, Biom and the two indices SSI-GY and TOL-GY. A set of 208 modern lines were investigated for haplotype diversity at these three loci. Five groups with different allelic combinations were identified (Figure 5). Their allelic effect on GY was then assessed when field tested under high temperatures along the Senegal River [35]. The haplotype class with positive alleles at all three loci had the highest GY average reaching 2381 kg ha−1 with a maximum value of 3856 kg ha−1. Genotypes of the haplotype classes with only two favorable alleles reached GY of 2199 and 2103 kg ha−1, while lines that only carried one positive allele 2103 and 2023 kg ha−1 (Figure 5). ANOVA confirmed that the haplotype group with all three positive alleles was significantly superior to the others.

**Figure 5.** Effect of different allele combinations of the significant loci on yield performance of 208 accessions tested under heat stressed conditions along the Senegal River. The circle indicates the average of each class over 2 years, and the whiskers show the standard error of the mean. The accessions were divided into five clusters based on their haplotype for three major QTLs: "+" mark the positive and "-" the wild-type alleles. Letters (a, b, ab) indicate significant differences between the clusters.

#### *3.6. Validation of Markers for Marker Assisted Selection*

To effectively deploy in breeding the most interesting QTLs via MAS, it is first required a step of validation using more affordable marker methodologies and in different genetic backgrounds and environments. A total of 20 MTA sequences linked to important agronomical and spike fertility traits were submitted for KASP primers design. Among these, only 14 could be successfully designed, and 11 identified a polymorphism within the validation set. Four showed significant (*p* < 0.05) correlation to the test phenotype (Figure 6). Three QTLs were represented by these four markers, AX-95260810 and AX-94432276 tagged QTL.ICD.Heat.08 on chromosome 5AL, AX-95182463 underlines QTL.ICD.Heat.09 on chromosome 5BL, and AX-94408589 tags QTL.ICD.Heat.10 on chromosome 6BS. The latter two QTLs are among the three main effect regions identified in this study (Figure 5). AX-95260810 reached 15% correlation to grain yield under heat, 74% accuracy, 43% sensitivity, and 100% specificity. Especially, its ability to identify 100% of non-heat tolerant entries is particularly remarkable. AX-95182463 and AX-94408589 also reached significant correlations of 14% and 32% for grain yield and biomass under severe heat, respectively, with sensitivities of 62% and 40%, accuracies of 30% and 65%, and specificities of 4% and 90%. Overall, AX-9526081 and AX-94408589 appeared as the most suitable for MAS application.

**Figure 6.** Kompetitive Allele Specific PCR (KASP) markers validation on an independent set of 94 elite lines of ICARDA tested under severe heat for grain yield and biomass. Correlation was measured between the BLUE for grain yield recorded along the Senegal River and the haplotype score. Accuracy, sensitivity, and specificity where determined using only the top 20 and worst 20 lines. AX-95260810 and AX-94432276 tag QTL.ICD.Heat.08, AX-95182463 tags QTL.ICD.Heat.09, AX-94408589 tags QTL.ICD.Heat.10.

#### **4. Discussion**

#### *4.1. Evaluation of the Phenotypic Performance of Yield and Yield Components under Normal and Heat Stress Conditions*

Several studies reported that wheat plants are very sensitive to elevated temperatures during flowering and grain filling phases [9,55,56], due to a reduction in seed development and fertility [56–58]. This study evaluated a set of durum wheat genotypes derived from a global collection for GY and yield components under heat and normal conditions. The genetic and phenotypic diversity shown by this set together with its relatively similar flowering time, promote it as an ideal panel to test heat tolerance. Further, the plastic tunnel method deployed here allowed to increase the temperatures well above 21 ◦C, the value that defines the absence of the stress [9]. A similar methodology was also successfully deployed by Corbellini et al. [54] to study the effect of heat shock proteins on technological quality characteristics. Compared to timely vs. delayed sowing experiments to simulate heat stress, the use of the plastic tunnel method avoids incurring false discovery due to changes in the phenological behavior of plants.

In the present study, a short and severe episode of heat stress was applied from the beginning of heading to the early dough stage, and resulted in 54% reduction in grain yield. This was in agreement with the study conducted by Ugarte et al. [59] that found a reduction of up to 52% when thermal treatment was applied via transparent chambers. Interestingly, our stress treatment caused an average temperature increase of 10 ◦C, which caused an average GY reduction of 5.4% for each 1 ◦C raise. This value is well within the 4.1% to 6.4% interval suggested by Liu et al. [10] for 1 ◦C raise in temperatures. GNSpk was the most affected trait (−61%) with the highest positive correlation to GY. This is in good agreement with previous studies that have shown that seed setting is the most sensitive parameter to heat stress, with a noticeable influence on yield [28,60–62]. Still, its non-linear relationship to yield confirms the complexity of the trait. Biom and HI were also found to have an influence on yield [63,64] with different relationships based on the occurrence of the stress. The presence of dissimilarities of the associations between the two treatments indicates clearly that there is a trade-off among the yield components as previously reported by Sukumaran et al. [65] for grain weight and grain number. Variation of one of the yield components affect the others positively or negatively. Compared to the simple regression, the additive model allowed to reveal the complexity of the relationship between GY and yield related traits.

The stress index SSI was developed by Fisher and Maurer [43] and modified by Nachit and Ouassou [44] as a useful indicator and a good parameter for selection. It measures the severity of the heat stress [66,67] and was also used in earlier studies in wheat to seek heat tolerant genotypes [23,68,69]. The TOL index is instead useful for selecting against yield depression, and it was used in several studies for heat or drought tolerance in wheat [27,44,67,70]. Improving heat tolerance should not be based on the use of these criterions alone as was suggested by Clarke et al. [71]. It is important to select simultaneously for good yield performance coupled with good adaptability (SSI < 1) and stability (low TOL) [44]. In that sense, the accessions Kunmiki, Berghouata1, Margherita2, and IDON37-141 originated from ICARDA durum wheat program, and Ourgh, a Moroccan cultivar, have been identified as high yielding genotypes that also show good heat stress tolerance based on the two indices.

#### *4.2. Dissection of Heat-Specific QTLs Associated with Yield-Related Traits and Stress Tolerance Indices*

The significant correlation identified between yield and its components were not linear in nature, and tend to change their mode of action based on the occurrence of the stress. Therefore, several physiological processes are simultaneously involved in protecting the wheat plant from the heat stress [72], and there is value in dissecting it into its genetic components. In this study GWAS was used to identify the genetic regions controlling the response of the various traits. To prevent the confounding effect that phenology-related loci might have [73], MTAs were identified for DTH and removed from downstream analysis. Additionally, flowering time was used as covariate in all analyses for the other traits. Very few MTAs for DTH were observed either in normal or stressed conditions due to the synchronized flowering of the entries used in this study. This indicated the absence of confounding effects between the two trials. i.e., almost all the accessions were exposed to the same conditions in each developmental phase [74] before imposing the stress.

Out of 12 QTLs identified, three occurred only when the heat stress was imposed, including indices. These three main genomic regions occurred on chromosomes 1AL, 5BL, and 6BS, and were considered as QTLs controlling heat tolerance. These three loci were confirmed by mean of haplotype analysis on a larger panel of modern lines (208 entries) field tested under severe heat along the Senegal River valley [35], to confirm that the presence of the positive alleles at all three loci provided a significant GY advantage of +182 kg ha−1 (+8%). The QTL on the long arm of chromosome 1A controlled GNSpk, HI, and TOL-GY, and it explained up to 27% of the phenotypic variation. In a study with double haploid population of bread wheat, Heidari et al. [75] identified a major QTL on the same chromosome (1A), influencing grain number per spike, grain weight per spike, and spikes/m2. However, their phenotypic assessment was not performed under heat stress, the marker systems used was different compared to our study and the locus was identified in the short arm of chromosome 1A. Therefore, it is quite difficult to align the results from that study to the current one. Another study had previously reported many MTAs on chromosome 1A detected for yield components under heat stress, but all were found to have a pleiotropic relationship with days to heading and were also located on the short arm of 1A [26],

instead of 1AL found here. A heat-specific QTL was also detected on the same chromosome in the short arm for spikelet compactness and leaf rolling in bread wheat [76]. An earlier study identified a QTL on 1AS for yield but associated with different stress conditions [77]. To the best of our knowledge, this is the first time that this region on 1AL is presented as associated to GNSpk, HI, and TOL-GY in durum wheat under heat stress conditions. The second major QTL region was detected on the long arm of chromosome 5B and found to be associated with GNSpk and the two indices SSI-GY and TOL-GY, contributing to 37% of the phenotypic variation. A region in the short arm of the same chromosome has been previously reported to be associated with grain number per square meter in bread wheat [76], and controlling thousand grain weight in durum wheat [27] under combined drought and heat stress. Shirdelmoghanloo et al. [25] and Acuna-Galindo et al. [78] reported loci for grain weight and other important traits on chromosome 5B under heat and non-heat conditions in hexaploid wheat. On the other hand, the same chromosome has been previously suggested to carry heat-specific QTLs for yield per se in bread wheat [26]. Sukumaran et al. [27] identified markers for heat susceptibility (HSI or SSI) and tolerance (TOL) indices for yield and grain number per square meter on the short arm of the chromosome 5B. Mason et al. [64] also detected QTL for HSI for kernel number on 5BL in bread wheat. The genomic region identified in this study on 5BL is likely to be a new QTL since no information has been reported earlier for this locus associated to GNSpk, SSI-GY, and TOL-GY in durum wheat and specific to heat stress, but we cannot exclude that it overlaps with previously reported QTLs. A third heat-responsive locus was identified on the short arm of chromosome 6B related to GY, SSI-GY, GNspk, HI, and Biom accounting for 36% of the phenotypic variance. An earlier study on bread wheat identified a locus on chromosome 6BS underpinning chlorophyll loss rates and heat susceptibility index for grain weight and chlorophyll loss rates under heat-stress conditions [25]. Under post-anthesis high temperatures stress, Vijayalakshmi et al. [20] reported a QTL on the short arm of chromosome 6B for senescence related traits in hexaploid wheat. McIntyre et al. [79] and Pinto et al. [21] reported QTLs on chromosome 6BL that were associated with many important traits (grain number per square meter and grain yield and water-soluble carbohydrate content) related to drought and heat tolerance. Ogbonnaya et al. [26] found a locus on the short arm of chromosome 6B for grain yield under heat stress in bread wheat. These previously reported QTLs in 6B could overlap with the one identified in this study, but they were either identified not in association with heat tolerance or detected in hexaploidy wheat. Therefore, this region is also assumed to have been reported for the first time here in relationship to heat tolerance for durum wheat. This locus affects multiple traits (GY, GNspk, HI, Biom, and two heat susceptibility indexes) and hence it is of good importance for deployment in breeding. The principal breeding objective is to develop varieties with high grain yield and stability when exposed to different stresses. However, grain yield is a complex trait controlled by many genes and strongly influenced by the environment [80–86]. Therefore, a good understanding of traits and underlying loci associated with tolerance to elevated temperatures is of a great importance for breeding new heat tolerant cultivars [87].

#### *4.3. Pyramiding Heat-Tolerant QTLS via MAS*

Three loci on chromosomes 1AL, 5BS, and 6BS showed an additive nature by means of haplotype analysis (Figure 5), revealing that only the combination of all three positive alleles generated a true yield advantage. Among the most heat tolerant elite lines identified here 'Kunmiki', 'Berghouata1', and 'Ourgh' confirmed to harbor the positive alleles for all three loci. This prompts their use in crossing schemes to pyramid the positive alleles, as well as the deployment of simple marker system to conduct MAS.

Axiom to KASP marker conversion and validation was attempted for 20 MTAs. Eleven KASP markers generated polymorphic haplotypes in an independent set of ICARDA elite lines. Four revealed a significant (*p* < 0.05) correlation to GY and biomass assessed under severe heat along the Senegal River Valley (Figure 6). In particular, AX-95182463 tags QTL.ICD.Heat.09 located on chromosome 5B and it revealed good correlation and sensitivity, but lacks in accuracy and specificity, and it is hence protected from Type II errors, but prone to Type I, with several elite lines wrongly identified as carrying the positive alleles. AX-95260810 tags QTL.ICD.Heat.08, linked to the two stress tolerance indices for GY (SSI-GY and TOL-GY) located on chromosome 5A. AX-94408589 tags QTL.ICD.Heat.10 located on chromosome 6B, and associated to several traits GNspk, HI, Biom, SSI-GY. In these two cases, the KASP markers explained 15% and 33% of the phenotypic variation of an independent validation set, with 100% and 90% specificity, and 74% and 65% accuracy, but medium sensitivity (43% and 40%). As such, these markers are protected against Type I errors (no false positive), but prone to Type II errors, with several elite lines identified as not carrying the positive allele while instead being tolerant to heat. Hence, while all converted KASP markers are prone to different types of errors, these three markers can be considered as validated and ready to be deployed in breeding. The combination of the three might represent a more stringent approach to protect against both types of errors. An additional nine QTLs were identified in this study, and their KASP conversion and validation are still ongoing and will require better targeted efforts to be achieved.

#### **5. Conclusions**

Heat stress causes a complex cascade of negative effects on the wheat plant, resulting in drastic reductions in grain yield. The deployment of heat tolerant varieties that will benefit greatly farmers requires first to enhance our understanding of this mechanism and loci governing it. Our study combined a discovery phase with a core set tested over two field seasons in Morocco under artificial heat-treatment with plastic tunnels, followed by a different confirmation set of germplasm grown for two seasons in Kaedi, Mauritania under severe natural heat, and completed with one final validation set tested one season in Kaedi. Our results confirmed that spike fertility (GNSpk) and maintenance of green leaves (Biom) are the most critical traits to drive tolerance to this stress, and hence should be the primary targets of durum wheat breeders. Further, the deployment of plastic tunnels proved to be a strategic methodology to study this stress and reveal its mechanisms without affecting the phenology of the plant. In addition, 12 loci were identified as responsible for controlling the main heat tolerance traits. Among these, three were activated only when the stress occurred and hence represent ideal targets for breeding. Two of these were validated into a KASP marker and are now ready for deployment via MAS, especially if associated with a third, also validated, KASP. Finally, three ICARDA elite lines and one Moroccan cultivar were confirmed as tolerant to heat, with high grain yield, and carrying positive alleles for three main QTLs. These are freely available and should be incorporated as crossing parents by other breeding programs. Altogether, this study has confirmed the key traits for heat tolerance as well as a new methodology to study it in durum wheat, it has revealed the main loci controlling these traits and proceeded to validate three of them for MAS, and it has also provided freely available elite lines to breed new cultivars better adapted to the stress.

**Supplementary Materials:** Table S1: List of durum wheat genotypes evaluated under plastic tunnel-mediated heat stress in the present study, Table S2: Markers associated with days to heading (DTH) under heat stress and normal conditions, Table S3: Sequence information of the KASP markers, Table S4: Pearson correlation matrix between all the measured traits under heat conditions (upper part) and normal (lower part) conditions. GY—Grain yield; Biom—Biomass; HI—Harvest index; Spkm2—Spikes per square meter: GNspk—Grain number per spike; TKW—Thousand kernel weight; DTH—Days to heading. \*, \*\* Significant at the 0.05 and 0.01 probability levels, respectively, Table S5: Correlation (r), linear regression estimated via ordinary least squares (OLS) and flexible regression estimated via regression additive model. (**a**) Under heat stress. (**b**) Under normal conditions, Table S6: Regression matrix between the haplotype of the peak markers for the 13 identified QTLs. \*, significant loci similarity at *p* < 0.05 consistent with homeologous relationship; \*\*, significant loci identity (*p* < 0.01) consistent with wrongly assigned genomic position, Figure S1: Plots of the additive regression model showing GNspk, biom, TKW, spkm2 and HI as the spline function of the target trait grain yield (GY). (**a**) Under heat stress. (**b**) Under normal conditions

**Author Contributions:** F.M.B., M.N. and K.E.H. conceived and designed the study. K.E.H. and F.M.B. performed the field experiment. A.T.S. performed the field experiment in the Senegal river. A.A. contributed in the genotyping. K.E.H. and F.M.B. analyzed the data. K.E.H. Wrote the original draft. K.E.H., B.B., A.F.M., A.A., M.N., and F.M.B. wrote or reviewed the manuscript. All authors read and approved the final manuscript.

**Funding:** This study was funded by the Australian Grains Research and Development Corporation (GRDC) project ICA00012: Focused improvement of ICARDA/Australia durum germplasm for abiotic tolerance, while the field work along the Senegal River was funded by the Swedish Research Council (Vetenskapsradet) U-Forsk2013 project 2013-6500, "Deployment of molecular durum breeding to the Senegal Basin: capacity building to face global warming" and U-Forsk2018 project 2017-05522, "Genomic prediction to deliver heat tolerant wheat to the Senegal River basin: phase II." The marker conversion work was covered by the International Treaty on Plant Genetic Resources for Food and Agriculture 2014-2015-2B-PR-02-Jordan: "An Integrated Approach to Identify and Characterize Climate Resilient Wheat for the West Asia and North Africa."

**Acknowledgments:** The authors wish to acknowledge the technical assistance provided by A. Rached and all ICARDA durum wheat program staff in handling field activities.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Screening of Diverse Ethiopian Durum Wheat Accessions for Aluminum Tolerance**

## **Edossa Fikiru Wayima 1, Ayalew Ligaba-Osena 2, Kifle Dagne 3, Kassahun Tesfaye 3, Eunice Magoma Machuka 4, Samuel Kilonzo Mutiga 4 and Emmanuel Delhaize 5,\***

Received: 15 July 2019; Accepted: 7 August 2019; Published: 9 August 2019

**Abstract:** Acid soils and associated Al3+ toxicity are prevalent in Ethiopia where normally Al3+-sensitive durum wheat (*Triticum turgidum* ssp *durum* Desf.) is an important crop. To identify a source of Al3+ tolerance, we screened diverse Ethiopian durum germplasm. As a center of diversity for durum wheat coupled with the strong selection pressure imposed by extensive acid soils, it was conceivable that Al3+ tolerance had evolved in Ethiopian germplasm. We used a rapid method on seedlings to rate Al3+ tolerance according to the length of seminal roots. From 595 accessions screened using the rapid method, we identified 21 tolerant, 180 intermediate, and 394 sensitive accessions. When assessed in the field the accessions had tolerance rankings consistent with the rapid screen. However, a molecular marker specific for the D-genome showed that all accessions rated as Al3+-tolerant or of intermediate tolerance were hexaploid wheat (*Triticum aestivum* L.) that had contaminated the durum grain stocks. The absence of Al3+ tolerance in durum has implications for how Al3+ tolerance evolved in bread wheat. There remains a need for a source of Al3+-tolerance genes for durum wheat and previous work that introgressed genes from bread wheat into durum wheat is discussed as a potential source for enhancing the Al3+ tolerance of durum germplasm.

**Keywords:** aluminum; resistance; Ethiopia; durum; hydroponics; soil acidity; *Triticum turgidum* ssp. *durum* Desf.

#### **1. Introduction**

Durum (*Triticum turgidum* ssp *durum* Desf.) wheat and bread wheat (*Triticum aestivum* L.) provide the Ethiopian population with a large proportion of its caloric intake [1,2]. Despite a general increase in production and productivity of wheat during the last two decades (1998–2016), Ethiopia imports a substantial amount of both durum and bread wheat. Durum production as a proportion of the total wheat produced in Ethiopia has declined over the past few decades from about 80% in the 1980s to an estimated current proportion of only about 20%. Despite a reduction in the proportion of grain produced in Ethiopia, durum constitutes 50% to 80% of the wheat grain imported in any given year [1], indicating that demand for durum remains strong. One of the abiotic stresses that may be contributing to the decline in durum production is the prevalence of acid soils in Ethiopia. Durum wheat is very sensitive of the toxic Al3+ found in acid soil and is the most sensitive of the small-grained crops to Al3+ toxicity [3]. By contrast, bread wheat has Al3+-tolerant alleles of *TaALMT1* [4] and *TaMATE1B* [5] along with other yet to be cloned Al3+-tolerance genes [6].

Worldwide, acidic soils are one of the most important limitations to agricultural production [7]. In Ethiopia, acid soils account for about 34% of agricultural land area that range from slight to strongly acidic soils [8]. Acid soils generally occupy the western part of the country extending from southwest to northwest, although strongly acidic soils occur mainly in the western part of the country including the lowlands. Acid soils are particularly prevalent in the highlands of Ethiopia [9] and the application of lime commonly results in improved yields of various crops [10,11]. The productivity of wheat over acidic areas of Ethiopia is low as compared to parts of the country where soil acidity is absent. For instance, in areas with strongly acidic soils, which occur widely in the western and southwestern parts of the country, the productivity of wheat is as low as 0.8–2.0 t/ha. These areas include West Wallaga, Illu-Ababora, Jimma, Gamo Gofa, Asosa, and Metekel zones. On the other hand, a relatively higher grain yield, ranging from 2.5 to 3.1 t/ha, is obtained in regions with near neutral soil pH including West Arsi, Arsi, Siltie, East Shewa, and Bale zones (Supplementary Materials file, Table S1).

Aluminum is solubilized in acid soils into the toxic Al3+ which in sensitive plants typically affects the viability of the root apex. Al3+ inhibits root growth resulting in reduced uptake of water and nutrients ultimately hindering plant growth and development [12]. In addition to its toxic effects on plant tissues, Al3+ affects nutrient availability within the soil. All of these effects significantly reduce crop yield. While management, primarily by application of lime, is important for neutralizing acid soils, this can be costly and it can take years to correct acidity at depth. A complementary strategy for improving crop production on acid soils is the use of Al3+-tolerant germplasm developed through breeding or genetic modification along with liming practices. In several crop species, variation in Al3+ tolerance has been identified and selective breeding programs can be implemented to increase production on acid soils [13–15]. However as noted above, durum wheat is one of the most Al3+ sensitive of the small-grained crops and shows little variation in Al3+ tolerance [16]. Bread wheat is hexaploid and possesses the A-, B-, and D-genomes whereas durum wheat is tetraploid possessing only the A- and B-genomes. Although, as noted below, a major gene for Al3+ tolerance is found on the D-genome of bread wheat, there are other genes for Al3+ tolerance in bread wheat found on the A- and B-genomes [6]. Early reports that identified Al3+-tolerant durum genotypes with a level of tolerance similar to that of bread wheat [17] can now be attributed to misidentification or contamination of durum grain stocks with bread wheat [18], highlighting the importance of verifying the genetic identity of germplasm. Bread wheat, in contrast to durum, shows a large variation in Al3+ tolerance and much of this is conditioned by alleles of the *TaALMT1* gene located on chromosome 4D [19] and *TaMATE1B* located on chromosome 4B. *TaALMT1* encodes an Al3+-activated membrane channel permeable to malate with the malate exuded by root apices binding toxic Al3+ to protect the developing roots. Another Al3+ tolerance gene in bread wheat is *TaMATE1B* that encodes a citrate transporter in the plasma membrane and in some bread wheat genotypes confers a lower level of tolerance than *TaALMT1* [20] whereas in other genotypes *TaMATE1B* appears to be the predominant Al3+ tolerance gene [21]. Both *TaALMT1* and *TaMATE1B* have been introgressed from bread wheat into durum wheat. This required use of the pairing homeologous (*ph1c*) mutant to enable *TaALMT1* to be transferred from chromosome 4D in bread wheat to chromosome 4B in durum wheat [20]. In contrast to bread wheat, for durum seedlings grown on acid soil, *TaMATE1B* appears to provide a greater level of Al3+ tolerance than *TaALMT1* for reasons that are not understood.

Ethiopia has been considered as a center of diversity for durum germplasm and landraces grown by farmers are a potential source of agronomically important genes [22,23]. A 90 K single nucleotide polymorphism (SNP) chip analysis supports the notion that durum wheat in Ethiopia is particularly diverse [24]. Furthermore, others undertaking similar phylogenetic analyses based on SNP markers have suggested that Ethiopia is a second center of origin for durum wheat [25]. We speculated that because of the diversity of the durum germplasm coupled with a strong selection pressure imposed by extensive areas of acid soils in Ethiopia useful levels of Al3+ tolerance could have evolved in Ethiopian

durum landraces. Cultivated durum genotypes sourced from various countries were found to be all Al3+ sensitive [16] whereas others [26] identified relatively tolerant lines in a similar population. The Ethiopian landraces are likely to be more diverse in their genetic makeup than the populations previously screened for tolerance and provide an opportunity to identify unique genes that confer tolerance of acid soils and Al3+ toxicity in particular. As an initial strategy it would be preferable to identify Al3+ tolerance in germplasm already adapted to Ethiopian conditions instead of introgressing genes from other sources. In this study we screened a diverse set of Ethiopian durum germplasm for Al3+ tolerance using a rapid hydroponic screen and then assessed selected lines more thoroughly in the hydroponic screen and finally in field trials on acid soil. Despite Ethiopia having conditions that would favor the evolution of Al3+ tolerance, the diverse durum germplasm was found to be Al3+ sensitive and highlighted the importance of verifying the identity of grain stocks. The absence of Al3+ tolerance genes in diverse durum germplasm suggests that genes encoding Al3+ tolerance found on the A- and B-genomes of bread wheat arose subsequent to the hybridization with the D-genome that produced hexaploid wheat.

#### **2. Materials and Methods**

#### *2.1. Germplasm*

A total of 595 durum wheat accessions obtained from the Ethiopian Biodiversity Institute (EBI: www.ibc.gov.et; Supplementary Materials file, Table S2), were screened in a series of non-replicated trials for their Al3+ tolerance using a rapid hydroponic screen. Based on their performance at this preliminary stage, accessions were classified as Al3+-tolerant, -intermediate, or -sensitive (see below). Subsequently all the tolerant as well as selected accessions of the intermediate and sensitive classes, totaling 150 accessions, were evaluated in a replicated hydroponic experiment and a field trial. For the experiments that screened previously identified Al3+-tolerant lines of durum the selected lines are described by Raman et al. [26]. A durum line (Langdon 4D (4B)) that has the 4B chromosome substituted by the 4D chromosome of bread wheat [27] was included as a positive control of a confirmed Al3+-tolerant durum line [18].

#### *2.2. Hydroponic Culture*

A rapid hydroponic screening was undertaken using an apparatus comprised of a plastic basin (for holding nutrient solution), a plate that held seed of different wheat accessions separately submerged in nutrient solution and an aeration system (Supplementary Materials file, Figure S1). The method is based on one described previously [16]. Seedlings grown submerged in the aerated nutrient solution remained viable and could be transplanted to soil when required to bulk up grain. Grain harvested from single plants of selected seedlings from this "preliminary" screen was then used in a replicated hydroponic experiment using the same growth conditions. For the replicated experiment, the average performance of 10 seedlings for each accession was determined and the experiment repeated three times.

The nutrient solution comprised of 500 μM KNO3, 500 μM CaCl2, 500 μM NH4NO3, 150 μM MgSO4, 10 μM KH2PO4, 2 μM FeCl3, and 5 μM of Al2SO4.18H2O as described previously [28]. The pH of the solution was adjusted to about 4.3 with 1 M HCl. Dry grain was immersed in the solution and seedlings allowed to grow submerged for 5 days in the solution culture with the nutrient solution changed every day to maintain pH and Al3+ concentration relatively constant. After 5 days seedlings had typically developed three seminal roots and the total root length was recorded to the nearest 0.1 cm.

We used conventional hydroponics with the same nutrient solution composition as described above to compare the Al3+ tolerance of a set of durum and bread wheat genotypes in a range of AlCl3 concentrations. Relative root length for the various genotypes was calculated after 3 days growth as root length in solution that contained Al3+ relative to root length in control solution that lacked Al3+. Errors associated with relative root length were calculated as described previously [29].

#### *2.3. Determining the Identity of the Al3*+*-Tolerant Accessions*

The morphological similarity between durum and bread wheat made it difficult to establish whether the identified Al3+-tolerant lines were durum or bread wheat by using phenotypic traits alone. Therefore, a polymerase chain reaction (PCR) assay that targeted the Dgas44 sequence was used to determine whether the lines were durum or bread wheat [30]. Dgas44 is a D-genome specific repetitive sequence that can be used to distinguish hexaploid wheat from tetraploid wheats that lack the D-genome [31]. The sequence of forward and reverse primers of Dgas44 marker, respectively, were 5 -CTTCTGACGGGTCAGGGGCAC-3 and 5 -CTGAATGCCCCTGCGGCTTAAG-3 .

Ten grain of bulked up samples used in the field trial along with one verified bread wheat cultivar (Enkoy) and one durum wheat accession (8317), distinguished by its reddish/pinkish seed color, were planted in pots. Enkoy variety was included as a positive control (i.e., possesses the D-genome since it is a bread wheat cultivar), while accession 8317 was used as a negative control since it is a known durum accession that lacks the D-genome. Young green leaves were collected separately from three individual plants for DNA extraction. Leaf samples were freeze-dried in liquid nitrogen and pulverized with a Geno/Grinder 2000 and genomic DNA was extracted with a ZymoResearch kit (Plant/Seed MiniPrep) following the manufacturer's protocol. The extracted genomic DNA was quantified with a NanoDrop 2000 UV-Vis Spectrophotometer (Thermo Scientific, Waltham, MA, USA) and PCR of samples undertaken as previously described [18].

#### *2.4. Field Experiment*

The field experiment was conducted at Bedi (38◦36 3" E, 9◦5 59" N), which is located in Watabicha Minjaro Kebele, Welmera District, West Shewa Zone, Oromia Regional State, Ethiopia (Supplementary Materials file, Figure S2). It is situated at approximately 35 km west of Addis Ababa and about 25 km away from Holeta town in the north-east. Crops are generally planted from around mid-June, though it varies with crop type and the time of onset of rains.

Samples of the reddish-brown soil to a depth of 20 cm were collected from the experimental plot with an auger in a regular pattern following a line transect. The soil samples were submitted the following day to JIJE LABOGLASS PLC (Addis Ababa, Ethiopia) for analysis of the major parameters using standard procedures as shown in Table 1. Characterized as a clay type, the pH of a water extract of the soils was found to be 4.92 indicating acidity and the likely presence of Al3+ toxicity. Note that a water extract will generally have a higher pH value than if a CaCl2 solution is used for the extraction. For the limed treatments 451 kg of CaCO3 was manually applied to the soil the year prior to the trial covering an area of 200 m2. The lime was obtained from Guder Lime Factory (Guder, West Shewa Zone, Ethiopia) through collaboration with Oromia Agriculture and Rural Development Bureau. The plot was then ploughed immediately and re-ploughed after a week to thoroughly mix the lime with the soil to a depth of about 20 cm so that Al3+ would be detoxified in the soil solution.

A total of 150 accessions were selected based on their response to Al3+ toxicity in the preliminary hydroponic screen described above. These accessions were planted during the major rainy season on limed and un-limed plots in the field using a randomized complete block design. The two treatments of the experiment (limed and un-limed blocks) were each replicated twice, and the replicated blocks were spaced 1 m apart such that the long edge of the blocks were arranged adjacent to one another. All 150 accessions were planted side by side in each block resulting in a total area of 30 × 2 m for each block. Individual accessions were planted within the blocks randomly in single rows 2 m long with 20 cm spacing between rows. Contamination of un-limed plots with lime through erosion and splashing was avoided by considering slope and spacing factors in the experimental layout. Specifically, plots at a higher slope were assigned as un-limed plots and were separated from the limed plots by 2 m. N and P fertilizers were applied after 3–4 weeks of planting as urea and di-ammonium phosphate at blanket recommendation rate (100 kg/ha) since there was no pre-determined site-specific fertilizer application rate. Hand-weeding was done at early seedling stage (about one month) and before the booting stage. Mature plants were harvested to measure total biomass and grain yield.

**Table 1.** Chemical and physical characteristics of the soils of the study area.

#### *2.5. Statistical Analysis*

Normality of distributions of the replicated hydroponic and field data were tested with SPSS version 24 (IBM SPSS statistics for Windows 2016, Version 24.0. Armonk, NY, USA: IBM Corp.) and RStudio (version 1.0.143: Integrated Development for R. RStudio, Inc., Boston, MA, USA. Retrieved from www.rstudio.com). The non-parametric (Kruskal–Wallis) and parametric ANOVA were computed either with RStudio or SPSS in order to determine whether there was significant difference between accessions for their performance in Al3+ treatments. The general linear model in which both accessions and replications were considered as fixed factors was used for analysis of variance. The median absolute deviation, a non-parametric statistic, was computed with RStudio to assess the level of variability attributed to the trait. Furthermore, post-hoc analysis was conducted using Tukey's honestly significant difference for parametric statistics, while a nonparametric post-hoc was done with Dunn's test, a Kruskal–Wallis test based post-hoc with "agricolae" package of RStudio. The coefficient of determination (r2) was calculated to examine the relationship between various variables.

#### **3. Results**

#### *3.1. Hydroponic Screen*

Durum wheat accessions (595) were initially screened for Al3+ tolerance using a rapid hydroponic method of submerged seedlings. These accessions sometimes comprised of two or more genotypes since there was a large variation in performance between individual plants of the accession (Figure 1). Furthermore, there were visually observable differences within an accession such as variation in grain color. To take this heterogeneity into account an accession was scored based on its best-performing seedling. The use of average performance of plants in representing an accession would have resulted in rejection of many accessions because of a poor average performance such that a single plant within the accession with an acceptable level of Al3+ tolerance would be lost.

Accessions were classified into three phenotypic classes based on their total root length as tolerant (≥3.1 cm), intermediate (≤3.0 cm but ≥2.1 cm) and susceptible (≤2.0 cm). Moreover, accession 6956 was included in the tolerant class because of its exceptionally long primary seminal root. Using this criterion 21 accessions were classed as tolerant, 180 as intermediate, and 394 as susceptible to Al3+ toxicity (Supplementary Materials file, Table S3).

**Figure 1.** Examples of seedlings grown for 5 days in Al3+-containing nutrient solution. Accessions 204399 and 214516 were scored as tolerant as was accession 234641 that was clearly segregating for tolerance. Accession 238131 was scored as sensitive. Individual seedlings identified within an accession to be tolerant were planted out to bulk up grain and the harvested grain used in a further hydroponic screen and a field trial.

A total of 150 accessions that included all 21 of the Al3+-tolerant class, 79 of the best performers of the intermediate class, and a random selection of 50 of the sensitive class were evaluated in three replications to more precisely characterize their Al3+ tolerance. In this experiment accessions were represented in each experiment by the average performance of ≥10 seedlings using grain harvested from individually selected seedlings. The performance of accessions in the replicated experiment using the same hydroponic method as the preliminary screen was well correlated with the preliminary screen (*r*2 = 0.72; Figure 2).

**Figure 2.** Relationship between total root length of the individual seedlings selected (total root length from the preliminary screen, *n* = 1 for each accession) that had the longest roots and the total root length of the resulting progeny of the selections showing the average of 10 or more seedlings (in three replications). Error bars indicate the standard error of the mean for the three replications.

#### *3.2. Field Experiment*

The set of 150 selections used in the replicated hydroponic screen was assessed in the field on an acid soil. As there was limited grain harvested from single plants of the selected seedlings, the accessions were assessed as single rows 2 m long on both an un-amended plot and an equivalent plot that had been amended with lime to neutralize the acidity. The grain yield and biomass means for the limed plots of all accessions combined was greater than the means of the acid plot indicating that

soil acidity was clearly present at this site. Grain yield on the acid site was reduced on average to only 18% of the limed site whereas biomass was reduced to 28% of the limed site (Figure 3).

**Figure 3.** Liming improves biomass and grain yield of Ethiopian accessions grown on an acid soil. Ethiopian accessions (150) were grown in single rows 2 m long on an acid soil and the same soil that had lime incorporated to a depth of 20 cm. Final biomass and grain yield were determined and data of all accessions were combined with error bars indicating the standard error of the mean (*n* = 150). Student's *t*-test indicated significant differences between treatments (\*\*\*; *p* < 0.001).

The performance of the Al3+-tolerant lines on the un-limed acid plot was remarkable and clearly distinct from that of the sensitive accessions. Figure 4 shows the relative (un-limed/limed) biomass and grain yields of the accessions selected for the field trial against the root length as determined in the preliminary screen.

**Figure 4.** Root length as determined with the rapid screen identifies Al3+-tolerant accessions based on both mature biomass (**A**) and grain yield (**B**). The preliminary root length is a measure of the total root length of the most tolerant seedling in the pool of seedlings assessed for each accession. The dashed

vertical lines denote the cut-off points for root lengths of seedlings initially classed as; tolerant: >3.0 cm; intermediate: 2.0–3.0 cm; sensitive: <2.0 cm. The single seedling was grown to seed and the resulting progeny grown in two replicates as 2 m rows on un-limed and limed plots. The mean of the replicates was used to calculate the relative biomass and relative grain yields where values for the un-limed plots were divided by the limed plots and multiplied by 100 to express data as a percentage. The solid line shows a linear regression with *r*2 = 0.56 for (A) and *r*2 = 0.52 for (B).

Expressing data in relative terms takes into account inherent differences in plant vigor when Al3+ is absent but a similar relationship was found when using unmodified data of only the acid plot (Supplementary Materials file, Figure S3). Biomass and final grain yield of genotypes grown in the acid trial were strongly correlated (Figure 5A). When biomass and grain yield data were combined within each grouping of lines selected from the preliminary screen (sensitive, intermediate, and tolerant), they showed rankings consistent with their grouping based on the screen (Figure 5B).

**Figure 5.** Relative biomass and relative grain yield are highly correlated and the classification of accessions into groups with the preliminary screen of Al3+-tolerant, -intermediate, and -sensitive is consistent with final biomass and grain yield. Biomass and grain yield of accessions grown on the un-limed plot are expressed as a percent of the limed plot using the procedure described in the legend of Figure 4. (**A**) Relationship between relative biomass and relative grain yield for all accessions grown in the field (*r*2 = 0.84). (**B**) Relative biomass and relative grain yield of all accessions allocated to the three classes were combined with error bars denoting the standard error and a one-way ANOVA of the data showed significant differences between all groups for each of biomass and grain yield (ANOVA on ranks *p* < 0.001 between classes; *n* = 47 for Sensitive, *n* = 83 for Intermediate, and *n* = 20 for Tolerant classes).

Since previous studies have found that stocks of durum wheat lines can be contaminated with bread wheat, we sought to verify that both Al3+-tolerant and -intermediate selections were indeed durum wheat. It can sometimes be difficult to distinguish the species based solely on the phenotypes so we used the Dgas molecular marker that is unique to the D-genome to distinguish the species. Using Dgas we found that all the Al3+-tolerant and -intermediate accessions were actually bread wheat (Figure 6), with only the most sensitive genotypes being durum wheat.

**Figure 6.** Examples of accessions screened by PCR for the Dgas sequence. Panel (**A**) is of Al3+-tolerant selections while panel (**B**) are all Al3+-sensitive selections except for the first lane after the DNA ladder. Most samples are shown as duplicates in lanes side by side except for the last lane of (A) (arrow) whose duplicate is the first lane of (B) (arrow). The lane on the left is a 100 bp DNA ladder with sizes of the smallest four markers shown on the right. A band at about 300 bp is indicative of the presence of the D-genome. Blk is a no DNA sample while C1 is a known durum wheat (accession 8317) and C2 is a known bread wheat variety (Enkoy). All lines (38 lines) that yielded more than 3.5 g per row in the field trial (Supplementary Materials file, Figure S3B) were analyzed and all possessed the 300 bp band indicating that they were hexaploid wheat.

#### *3.3. Other Potential Sources of Al3*+*-Tolerant Durum Wheat*

As discussed in the Introduction, the Al3+ tolerance of durum has been enhanced by introgression of genes from bread wheat although not into Ethiopian germplasm [20]. As an alternate source of genes, a previous report identified three out of 420 tetraploid genotypes screened that were comparatively Al3+ tolerant although they did not approach the tolerance shown by a tolerant bread wheat used as a check line [26]. The lines had been confirmed as tetraploid with a molecular marker specific for the D-genome, so these genotypes were a potential source of genes that could be used in direct crosses to Ethiopian lines. To establish if these three lines had a useful level of Al3+ tolerance, we assessed their performance in hydroponic culture against bread wheat lines that varied in tolerance as well as a 4D (4B) substitution line of durum. However, all three lines were rated as sensitive with the best performer having a similar level of tolerance as Al3+-sensitive bread wheat (Figure 7). By contrast, a durum line where the 4B chromosome was substituted by the 4D chromosome of hexaploid wheat had a level of Al3+ tolerance comparable to the tolerant hexaploids.

**Figure 7.** Relative Al3+ tolerance of durum wheat genotypes previously identified as being relatively Al3+ tolerant in comparison to a durum substitution line with chromosome 4D of bread wheat and bread wheat lines possessing Al3+-tolerant and -sensitive alleles of the *TaALMT1* gene. Bellaroi is an Australian durum wheat cultivar; lines 202, 249, and 267 are accessions previously classified as being Al3+-tolerant durum genotypes [26]; Langdon is the parental durum cultivar used to generate the 4D substitution line; 4D (4B) is a durum line in the Langdon background where chromosome 4B has been substituted with chromosome 4D; Egret is an Al3+-sensitive bread wheat cultivar; ET8 and Carazinho are Al3+-tolerant bread wheat lines. Root length is expressed as percent relative to a zero control and error bars indicate cumulative errors of three to five seedlings.

#### **4. Discussion**

#### *4.1. A Rapid Screen Proves Robust and Correlates with Grain Yield in Field Trials when Ranking Wheat Germplasm for Al3*+ *Tolerance*

In this work we screened Ethiopian durum accessions for Al3+ tolerance first using a high throughput method with hydroponics, then with a replicated hydroponic screen and finally in the field on acid soil. Ethiopian durum accessions that include landraces are reported to show great diversity for many traits and molecular analyses have verified the unique nature of this germplasm [24]. However, despite this diversity our work shows that there does not appear to be useful levels of Al3+ tolerance in Ethiopian durum wheat and that seedlings identified to be Al3+ tolerant were in fact contaminating bread wheat. It should be noted that many of the accessions classed as tolerant only had a few Al3+-tolerant seedlings indicating that the durum grain stocks had been contaminated with bread wheat grain. This was verified by the use of a molecular marker that confirmed Al3+-tolerant seedlings to be bread wheat. We conclude that the rapid screen was sufficiently robust as a preliminary screen to identify Al3+-tolerant bread wheat contaminants. The accessions ranked as having intermediate tolerance in the preliminary screen could still have conferred a useful level of Al3+ tolerance for durum if they ranked similarly in the field. However, in the field the best performing accessions ranked as having intermediate Al3+ tolerance from the preliminary screen were all identified as bread wheat (Figure 6 and Supplementary Materials Figure S3B). Nevertheless, there exist Al3+-sensitive genotypes of bread wheat so the rapid screen does not replace the use of a molecular marker in identifying bread wheat lines. Verifying the identity of durum wheat grain stocks is of particular importance for germplasm banks but should also be confirmed by researchers. Bread wheat genotypes in the past have been mistakenly identified as durum wheat and the most relevant to this study was the incorrect identification of Al3+-tolerant durum lines that were subsequently found to be bread wheat [17,18]. In some cases, bread wheat has morphological characteristics in the field that are similar to durum wheat emphasizing the need for molecular analysis to establish the species identity [32].

The absence of Al3+ tolerance in unmodified *Triticum turgidum* ssp *durum* Desf. is consistent with previous studies that screened durum cultivars with presumably a lower level of diversity than the Ethiopian germplasm and also showed no or only a comparatively low level of Al3+ tolerance [16]. The simple and rapid screen undertaken with minimal equipment over only 5 days growth with seedlings submerged in nutrient solution showed a remarkable consistency in classifying the relative tolerance of germplasm grown in the field and measured for mature biomass and even grain yield (Figure 5B). Although there was a good general agreement of the hydroponic screens with performance of accessions in the field, there were exceptions. For instance, several accessions rated as having intermediate tolerance with the preliminary screen were the best or amongst the best for grain yield in the field (Figure 4B). This finding emphasizes the importance of verifying selections on acid soil whether in the field or in pots and that while hydroponic screens allow for rapid assessment of germplasm, the relative tolerance of genotypes can differ when grown in soil. For instance, a durum line carrying the *TaMATE1B* gene introgressed from hexaploid wheat (see below) shows marginal Al3+-tolerance in hydroponics that is considerably less than a line carrying *TaALMT1*, whereas the situation is reversed when the lines are grown on acid soil, with *TaMATE1B* lines outperforming *TaALMT1* lines [20].

#### *4.2. The Evolution of Al3*+ *Tolerance in Bread Wheat Occurred Subsequent to the Hybridization of the D-Genome*

Here we show that durum lines previously reported to be relatively Al3+-tolerant [26] had a level of tolerance similar to sensitive wheat (Figure 7), confirming the absence of useful Al3+ tolerance in durum wheat germplasm (excluding lines where genes from hexaploid wheat have been introgressed as discussed below). This observation is puzzling given reports of the presence of Al3+ tolerance genes on the A- and B-genomes of bread wheat. Although *TaALMT1*, the major gene for Al3+ tolerance in bread wheat, is located on chromosome 4D, there are several loci on the A- and B-genomes reported to confer Al3+ tolerance [6]. For example, *TaMATE1B* the only Al3+ tolerance gene other than *TaALMT1* that has been cloned from bread wheat, is located on chromosome 4B [33]. We surmise that the absence of Al3+ tolerance genes in diverse durum germplasm suggests that the multiple Al3+ tolerance genes found in bread wheat arose subsequent to the hybridization of the D-genome with the A- and B-genomes some 10,000 years ago. Transposable elements have been shown to enhance the level of expression of genes encoding transport proteins such as those of the MATE family [34] that confer Al3+ tolerance. The event of polyploidization commonly results in genome instability including the activation of transposable elements [35]. It is conceivable that activation of transposons when the hexaploid was formed has been key in enabling bread wheat to evolve Al3+ tolerance.

#### *4.3. Strategies to Enhance the Al3*+ *Tolerance of Ethiopian Durum Wheat*

To date the only durum germplasm verified to be Al3+ tolerant are those that were developed by introgression of Al3+ tolerance genes from bread wheat into durum wheat [20]. In those cases, the *ph1c* mutant was used to introgress the *TaALMT1* gene located on chromosome 4D of bread wheat and steps taken to avoid the hybrid necrosis that can occur when crossing bread to durum wheat to introgress the *TaMATE1B* gene on chromosome 4B. This germplasm is not considered to be genetically modified since it was developed by so-called natural means and can therefore be used in the field without restrictions. The germplasm is a source of Al3+ tolerance genes in a durum background that could be used to introgress one or both genes derived from bread wheat into selected Ethiopian durum germplasm. Many of the farms that grow durum wheat in Ethiopia are small holdings and it has been found that landraces otherwise known as "farmer's varieties" perform better than "improved" germplasm under many situations [36]. With this in mind it may be useful to cross the Al3+-tolerant germplasm described above into selected landraces adapted to local regions. Using local landraces as recurrent parents in multiple backcrosses while tracking the presence of Al3+ tolerance genes using

molecular markers should maintain the germplasm with various valuable landrace traits that could otherwise be lost [37]. The *TaMATE1B* gene seems a preferred candidate at this stage since it appears to be more effective than *TaALMT1* in a durum background and there exists a co-dominant marker that can be used to track the tolerant allele [20]. A recent publication has shown that the *TaMATE1B* gene introgressed into a durum genetic background confers a marked ability of roots from mature plants to withstand Al3+ toxicity when grown in an acid soil [38]. A simple backcrossing program where a single gene is introgressed into landraces could be useful in establishing whether or not *TaMATE1B* can improve durum grain production of landraces on acid soils of Ethiopia. We speculate that introgressing both *TaALMT1* with *TaMATE1B* may provide the greatest level of Al3+ tolerance to durum wheat as has been found for some bread wheat genotypes.

#### **5. Conclusions**

Here we show that despite the diversity of Ethiopian durum germplasm and large regions of acid soils, conditions conducive for the evolution of Al3+ tolerance, a useful level of Al3+ tolerance was not detected. A high throughput screen identified Al3+-tolerant seedlings within accessions but all of them were shown to be contaminating bread wheat. This finding highlights the importance of correct identification of germplasm and confirmation of species identity with molecular markers particularly for closely-related species. The finding that Ethiopian durum germplasm is Al3+ sensitive is consistent with Al3+ tolerance having evolved in bread wheat subsequent to the hybridization of the D-genome with the A- and B-genomes. A field trial on acid soil showed biomass and grain yields that correlated with classifications based on a high throughput screen, confirming the utility of the rapid screen for preliminary assessment of germplasm.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4395/9/8/440/s1, Figure S1: Photographs of the components of the hydroponic system. (A) The plastic cups to which dry grain was added. The holes allowed aeration of grain while the cups themselves separated the accessions from one another. The slots in the plastic stands at the end of a series of cups enabled the long air stone to be inserted and held firm. (B) Top view showing air stones and pumps along with the plastic basin used to hold the nutrient solution. (C) Top view of assembled equipment showing plastic cups with a partitioning plate inserted within each cup so that two accessions could be placed into each cup and separated from one another. Nutrient solution was added to the basin to fill and cover the cups that held the grain. The holes in the cups allowed sufficient aeration for growth of seedlings. Nutrient solution was changed every day to maintain pH at about 4.3 and Al2SO4 at 5 μM. Figure S2: Location of the study area in Ethiopia. Figure S3: Relationships between Al3+ tolerance determined with a rapid screen and (A) biomass or (B) grain yield for durum accessions grown on an acid plot in the field. Table S1: Wheat production across zones and special districts of Ethiopia showing yields averaged over 5 years. Table S2: Passport data of durum accessions used in hydroponics screening for Al3+ tolerance. Table S3: Preliminary screen of Al3+ tolerance of 594 Ethiopian durum wheat accessions grown in Al3+-containing nutrient solution where total root length after 5 days growth was used to classify genotypes.

**Author Contributions:** Conceptualization, E.F.W.; methodology, E.F.W. and E.D.; formal analysis, E.F.W. and E.D.; investigation, E.F.W., E.D., A.L.O.; resources, E.F.W., A.L.O., K.T., E.M.M.; data curation, E.F.W.; writing—original draft preparation, E.F.W. and E.D.; writing—review and editing, A.L.O., and K.D.; supervision, E.D., A.L.O., K.D., K.T. and S.K.M.; project administration, K.T. and E.M.M.; funding acquisition, E.F.W., A.L.O. and K.T.

**Funding:** The hydroponic and field experiments were sponsored by Addis Ababa University and the Eastern Africa Agricultural Productivity Program, respectively, while the Dgas analysis was financed by the BecA-ILRI Hub through the Africa Biosciences Challenge Fund (ABCF) program. The ABCF Program is funded by the Australian Department for Foreign Affairs and Trade (DFAT) through the BecA-CSIRO partnership; the Syngenta Foundation for Sustainable Agriculture (SFSA); the Bill & Melinda Gates Foundation (BMGF); the UK Department for International Development (DFID) and; the Swedish International Development Cooperation Agency (SIDA).

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Genetic Advance of Durum Wheat Under High Yielding Conditions: The Case of Chile**

#### **Alejandro del Pozo 1,\*, Iván Matus 2, Kurt Ruf 1,2, Dalma Castillo 2, Ana María Méndez-Espinoza 1 and María Dolores Serret 3**

Received: 22 June 2019; Accepted: 13 August 2019; Published: 15 August 2019

**Abstract:** In Chile, durum wheat is cultivated in high-yielding Mediterranean environments, therefore breeding programs have selected cultivars with high yield potential in addition to grain quality. The genetic progress in grain yield (GY) between 1964 and 2010 was 72.8 kg ha−1 per year. GY showed a positive and significant correlation with days to heading, kernels per unit ground area and thousand kernel weight. The gluten and protein content tended to decrease with the year of cultivar release. The correlation between the δ13C of kernels and GY was negative and significant (−0.62, *p* < 0.05, for all cultivars; and −0.97, *p* < 0.001, excluding the two oldest cultivars). The yield progress (genetic plus agronomic improvements) of a set of 40–46 advanced lines evaluated between 2006 and 2015 was 569 kg ha−1 per year. Unlike other Mediterranean agro-environments, a longer growing cycle together with taller plants seems to be related to the increase in the GY of Chilean durum wheat during recent decades.

**Keywords:** agronomic traits; carbon isotope; days to heading; grain quality; yield components

#### **1. Introduction**

Durum wheat (*Triticum turgidum* L. ssp. *durum*) covers ~17 million hectares worldwide, which is less than 10% of the total wheat area. However, its importance for human consumption is very high because it is used for making pasta, couscous, burghul and firik [1]. According to the International Grain Council, the largest producers of wheat in the world are the European Union, Canada, the United States, Turkey and Algeria.

For the production of high-quality durum wheat, dry environments are necessary, with warm days and cold nights during the growing season so that large grains are obtained with yellow color, vitreous kernels (more than 95%), hard texture and high test weight (about 82 kg hL−1), alongside high protein content (greater than 10%) and strong gluten (greater than 30% wet gluten), which gives elasticity to dough for industrial use [2]. In Chile, durum wheat is grown in Mediterranean climate environments from the Valparaíso Region (32 ◦S) to the Biobio Region (37 ◦S), but mostly under irrigation conditions or in areas where rainfall is sufficient to satisfy most or all of the crop potential evapotranspiration. The sowing area has increased from 9600 ha in 2001 to 27,000 ha in 2015, and the average yield for 2011–2015 was 6.7 Mg ha−1 [3].

Wheat yields in different regions of the world have increased greatly since the 1960s as a result of genetic improvement and better agronomic practices [4]. With the Green Revolution, breeding programs have seen the introduction of semi-dwarfing genes that interfere with the action or production of gibberellin [5], leading to a reduction in plant size and an increase in the partitioning of the above-ground biomass towards spikes and grains [6,7]. In bread wheat, the genetic gain in grain yield (GY) was positively correlated with harvest index and the number of grains per spike and per m2 [6,8]. In durum wheat, a Spanish study conducted with 12 cultivars from Italy and 12 from Spain released between 1930 and 2000 showed that the changes in grain yield were also associated with increases in the harvest index and number of grains per m2 [9]. Another study carried out on 14 cultivars released in Italy between 1900 and 2000 indicated that the total aerial biomass had not changed and that the increase in yield was associated with a reduction in plant height and an increase in the harvest index and number of grains per m2 [10]. These studies in durum wheat have been conducted under rain fed conditions and yields were below 6 Mg ha−1, however, there is no information about the genetic progress in durum wheat in high yield potential environments (>10 Mg ha−1) and how grain quality traits have been affected. Moreover, these works have not focused on studying exclusively the trends in breeding advances of post-Green-Revolution (i.e., semi-dwarf) durum wheat cultivars during the last half-century. This is despite the importance of this issue in the context of climate change and the fact that at least for bread wheat, there are studies reporting a stagnation in yields (or at least a drastic decrease in genetic advance) during the last decade [11].

In bread wheat, grain protein content, sedimentation value and wet gluten have increased in modern cultivars [8,12–14]. In durum wheat, modern dwarf and semi-dwarf cultivars have a higher gluten index compared to landraces or traditional Mediterranean cultivars [15]. Subira et al. [16] also reported significant changes in grain quality traits in a historical series of 24 durum wheat cultivars released in Italy and Spain in different periods of the 20th Century, particularly in gluten strength, sedimentation index and yellow color index. High protein content and 'strong' gluten are necessary to process semolina into a suitable final pasta product.

Physiological changes associated with breeding advances have also been reported for bread and durum wheat. For instance, modern bread wheat presented higher stomatal conductance (on an area basis) and carbon isotope discrimination (Δ13C; or a lower carbon isotope composition, δ13C), and lower oxygen isotope composition (δ18O) than older varieties [8,17–19]. In durum wheat, modern varieties have higher Δ13C (or lower δ13C) compared to landraces [20–22], although no clear differences were found for δ18O [21]. However, no information exists for durum wheat growing in a high-yielding Mediterranean environment.

The aim of this work was to analyze a) the changes in agronomic traits, grain quality and isotope composition in a set of ten durum wheat cultivars released in Chile between 1964 and 2010; and b) the progress in grain yield, plant height and test weight in selected advanced lines from the Instituto de Investigaciones Agropecuarias (INIA)-Chile breeding program. The experiments were conducted in a high-yielding Mediterranean environment between 2006 and 2015.

#### **2. Materials and Methods**

#### *2.1. Experimental Site, Plant Material and Growing Conditions*

The experiments were conducted at the Santa Rosa experimental field station (36◦32 S, 71◦55 W; 220 m.a.s.l.) of the Centro Regional de Investigación (CRI)-Quilamapu, Instituto de Investigaciones Agropecuarias (INIA). The climate corresponds to a humid Mediterranean type. During the experimental period (2006–2015), the monthly minimum average temperature was 3.1 ◦C (July) and the maximum 29.6 ◦C (January), and the average annual precipitation was 903 mm (Supplemental Table S1). The soil was a sandy loam, humic haploxerands (Andisol). Soil chemical characteristics of the top 10 cm were: pH 6.0, 8.87 mg kg−1 of N-N03; 17.05 mg kg−1 of P (Olsen), 0.45% of N-total, 4.5% of C and 0.33, 5.75, 0.65 and 0.48, cmol kg−1 of available K, Ca, Mg and Na, respectively [23].

Two different experiments were conducted. In the first experiment, ten cultivars released by the INIA breeding wheat program from 1964 to 2010 (Table 1) were evaluated during three consecutive seasons (2010 to 2012). The INIA cultivars derive from germplasm introduced from the International Maize and Wheat Improvement Center (CIMMYT) and probably all of them have the *Rht-B1* gene.

**Table 1.** Cultivars of durum wheat released by the wheat breeding program of the Chilean INIA between 1964 and 2010.

1 Year of cultivar release; 2 experimental line.

The experimental design was a complete block with four replications. Each plot consisted of five rows of 2.5 m length and 0.2 m apart. Sowing dates were in August of each year and the sowing rate was 220 kg ha−1. Fertilization consisted of 1.5 t ha−1 of lime (88%–90% CaCO3) before sowing, 260 kg ha−1 of diammonium phosphate (46% P2O5, 18% N), 200 kg ha−1 of potassium magnesium sulfate (22% K2O, 18% MgO, 22% S), 90 kg ha−1 of potassium chloride (60% K2O), 10 kg ha−1 of boronatrocalcite and 3 kg ha−1 of zinc sulfate (35% Zn) at sowing. After sowing, an extra 133 kg ha−1 of urea (46% of N) was applied at tillering initiation (Zadoks 20; [24]) and 201 kg ha−1 at the first node (Zadoks 31). Plots were furrow irrigated according to the needs of the crop (3–4 irrigations of ~50 mm each, per season). Weeds were controlled using the pre-emergence herbicide Bacara Forte 360SC, Bayer Crop Science (800 mL ha−1; 12:12:12% w/v a.i. of flufenacet/flurtamone/diflufenican) and the post-emergence Ajax, Anasac, Chile (10 g ha−1; 50% w/w a.i. of metsulfuron-methyl) and MCPA 750 SL, Anasac, Chile (800 mL ha−1; 95% w/v a.i. of 2-methyl-4-chlorophenoxyacetic acid)). Since the oldest cultivars showed susceptibility to rust (*Puccinia striiformis* and *Puccinia triticina*), two applications were made of the foliar fungicide Juwel-Top, Basf (100 mL ha−1; 12.5:12.5:15% a.i. of kresoxim-methyl/epoxiconazole/phenopropimorph). These applications were made before symptoms appeared, to avoid any interference of these diseases in the development of the plants.

In the second experiment, a selection of 46 advanced lines (F6–F8) of durum wheat from the breeding program (Durum Yield Nursery) and four check cultivars (Llareta-INIA, Corcolén-INIA, Lleuque-INIA and Queule-INIA) were tested each year from 2006 to 2015. Two trials of 25 genotypes each, including check cultivars, were established each year in an α-lattice design with five incomplete blocks per replicate, each block containing five genotypes. There were four replicates per genotype. The plots consisted of five rows of 2 m length and 0.20 m between rows. The seed rate was the equivalent of 220 kg ha−1. The sowing date was August of each year. Crop fertilization and weed control were as recommended for each year. Plots were furrow irrigated according to crop need (3–4 irrigations of ~50 mm each, per season). These trials were regularly conducted by the breeding program in order to test the most promising advance lines in comparison with the commercial cultivars (check cultivars); those advance lines with outstanding performance were evaluated for more than one year, and the rest were replaced by new ones. As a consequence, the set of advance lines evaluated in each year was composed of different elite genotypes.

#### *2.2. Agronomic Traits*

In Experiment 1 the following traits were evaluated: (a) Days from emergence to heading (DH) through periodic observations (twice per week), when approximately half of the spikes in the plot had already extruded; (b) the number of spikes per m2 (SM2) by counting the spikes in a 1.0 m length

of a row; (c) the harvest index (HI), determined from a sample from the 1.0 m row at maturity and calculated as the ratio of grain dry weight to total above ground dry weight; (d) the number of kernels per spike (KS) and thousand kernel weight (TKW) from 25 spikes taken at random from each plot and (e) the number of kernels per m2 (KM2) calculated as SM2 × KS. In Experiments 1 and 2, plant height (PH) from the ground to the top of the spike, excluding awns, was measured at maturity, and GY was assessed by harvesting 2 m2 (five rows, 2 m long).

#### *2.3. Grain Quality*

The test weight was evaluated in Experiments 1 (2010) and 2 (2006 to 2015), in samples of wheat free of impurities (obtained from each genotype and replicate) using a 250 cc Schopper scale (Louis Schopper, Germany). In addition, grain samples obtained from the genotypes and replicates evaluated in Experiment 1 (in 2010) were ground in mill for wet gluten and protein content determination. Wet gluten content was determined according to the International Approved Methods of Analysis (AACCI Method 38–12.02) in 10 g of pure flour mixed with 5.5 mL of a 2% saline solution, which was homogenized and then placed in a gluten washer (Glutomatic® 2200, Perten Instruments, USA) for 5 min; then the wet gluten was weighed and expressed as a percentage of the amount of pure flour. Protein content (%) was also determined in ground grain samples placed in a quartz cuvette and the reflectance spectrum between 800 and 2500 nm was determined using near infrared reflectance spectroscopy (NIRS), Bruker, USA. Yellow berry incidence was assessed on 100 g of kernels, separating and weighing the affected grains and then expressed in percentage.

#### *2.4. Total N Content and C and N Isotope Analyses*

Measurements were performed in mature grains harvested in 2011 (Experiment 1). The total N content was analyzed using an elemental analyzer (Flash 1112 EA; ThermoFinnigan, Bremen, Germany). The stable carbon (13C/ 12C) and nitrogen (15N/ 14N) isotope ratios of the same mature grains were determined in the same elemental analyzer coupled with an isotope ratio mass spectrometer (Delta C IRMS, ThermoFinnigan, Bremen, Germany). Nitrogen was expressed as a concentration (g N per g of dry weight) and atropine was used as a system check in the elemental analyses of nitrogen. The 13C/ 12C ratios of plant material were expressed in δ notation: δ13C = ( 13C/ 12C) sample/( 13C/ 12C) standard – 1, where 'sample' refers to plant material and 'standard' of known 13C/ 12C ratios. The 15N/ 14N ratios were also expressed in δ notation (δ15N) using international secondary standards of known 15N/ 14N ratios. More details are described in del Pozo et al. [8]. Measurements were performed at the Scientific Facilities of the University of Barcelona.

#### *2.5. Data Analysis*

Complete block analysis of variance (ANOVA) were performed for the set of cultivars evaluated in Experiment 1 using IBM SPSS Statistics software (SPSS Inc, USA). In addition, correlation analyses were performed between the year of cultivar release and agronomic, grain quality and isotope composition traits, and among the different traits. Trends for grain yield, plant height and test weight of 46 advanced lines and cultivars evaluated from 2005 to 2015 in Experiment 2 are also presented.

#### **3. Results**

#### *3.1. Agronomic Traits in Cultivars Released During the Past Six Decades*

Days to heading differed significantly among cultivars and also the year × cultivar interaction was significant (Table 2); it reduced in the 1990s, but then increased in the 2000s (Figure 1A). Plant height was significantly (*p* < 0.001) reduced from 108 cm in the 1960s to 90 cm in the 1970s, with a slight increase in 2010 (Table 2; Figure 1B).

**Table 2.** Mean sum of squares of the analysis of variance (ANOVA) for agronomic traits of ten durum wheat cultivars cultivated during three growing seasons (2010–2012).

Level of significance is indicated in bold (*p* < 0.01) and cursive (*p* < 0.05). DH: Days to heading; GY: Grain yield; PH: Plant height; SM2: Number of spikes per m2; KS: Kernels per spike; KM2: Kernel number per m2; TKW: Thousand kernel weight; HI: Harvest index.

**Figure 1.** Relationships between the year of release of ten durum wheat cultivars and: Day to heading (**A**), plant height (**B**), grain yield (**C**), harvest index (**D**), number of spikes per square meter (**E**), number of kernels per spike (**F**), number of kernels per square meter (**G**) and thousand kernel weight (**H**). Values correspond to the average (±SE) of three growing seasons (2010–2012) except for HI, which was determined in 2010 and 2011. The oldest (1964) cultivar (open circle) was not considered in the regressions. Mean values of cultivars for each year of evaluation are shown in Supplemental Table S2.

GY exhibited a positive and linear relationship with the year of cultivar release (*R*2 = 0.66; *p* < 0.001), this analysis excluding the oldest cultivar (1964; Figure 1C). The rate of increase in GY after 1960 was 72.8 kg ha−1 per year, and excluding the oldest cultivar it was 63.4 kg ha−1 per year. The SM2 of the ten cultivars ranged between 450 and 595 and increased significantly (*R*2 = 0.43; *p* < 0.05) with the year of release (Figure 1E). The HI was 0.36 in the 1960s and increased to 0.45–0.53 in the 1970s and onwards, whereas TKW was 51.1 g in the 1960s and rose to 57–68 g after the 1970s, but neither trait was correlated with the year of cultivar release (Figure 1D,H). Similarly, KS was not correlated with the year of release (Figure 1F), but KM2 increased significantly (*R*2 = 0.40; *p* < 0.05) with the year of cultivar release (Figure 1G).

The correlation matrix among the agronomic traits of the 10 cultivars evaluated during three growing seasons indicated that days to heading exhibited a positive and significant correlation with GY (*p* < 0.05) and KS (*p* < 0.01), and GY showed a positive and significant correlation with KM2 (*p* < 0.05) and TKW (*p* < 0.001; Table 3). Plant height was not correlated with GY. However, plant height showed a negative and highly significant (*p* < 0.001) correlation with TKW and HI. SM2 had a positive correlation with KM2 but a negative correlation with KS.

**Table 3.** Correlation matrix among agronomic traits evaluated in ten cultivars during three growing seasons (2010–2012).

\*: *p* < 0.05; \*\*: *p* < 0.01; \*\*\*: *p* < 0.001 DH: Days to heading; GY: Grain yield; PH: Plant height; SM2: Spike number per m2; KS: Kernels per spike; KM2: Kernel number per m2; TKW: Thousand kernel weight; HI: Harvest index.

#### *3.2. Grain Quality and Kernel Isotope Composition in Cultivars Released During the Past Six Decades*

The test weight increased curvilinearly with the year of cultivar release (Figure 2A). The gluten and protein content tended to decrease with the year of cultivar release, although the correlations were not significant (Figure 2B,C). Yellow berry was higher in two cultivars, but there was no clear pattern with the year of cultivar release (Figure 2D).

The relationships between the year of cultivar release and N concentration or δ15N in kernels were not significant (Figure 3A,B). The δ13C of kernels tended to decrease with the year of cultivar release, although the correlation was not significant (Figure 3C). In addition, δ13C was negatively correlated (*r* = −0.62; *p* < 0.05) with GY, but δ15N was not correlated (*r* = 0.03; *p* > 0.05).

**Figure 2.** Relationships between the year of release of 10 durum wheat cultivars and kernel test weight (**A**), wet gluten content (**B**), protein content (**C**) and yellow berry (**D**), determined in 2010. Values correspond to the average (±SE) of four replicates.

**Figure 3.** Relationships between the year of release of 10 durum wheat cultivars and the (**A**) total nitrogen, (**B**) natural abundance of 15N (δ15N) and (**C**) carbon isotope composition (δ13C) in kernels, determined in 2011. Values correspond to the average (±SE) of four replicates.

#### *3.3. Agronomic and Grain Quality Traits in Advanced Lines During the Last Decade*

The GY and plant height of advanced lines increased from 2006 to 2015, reaching a maximum in 2011 with averages of 12.7 ± 0.8 Mg ha−1 and 96 ± 4.2 cm, respectively (Figure 4A,B). GY was highly correlated with plant height (*r* = 0.85; *p* < 0.001). The check cv. 'Corcolén' followed a similar trend to the advanced lines. The average GY of advanced lines and cultivars had a positive and significant (*R*2 = 0.50; *p* < 0.001) relationship with the year of evaluation; the regression analysis indicated that the rate of increase in GY between 2006 and 2015 was 569 kg ha−1 per year. The test weight did not increase during this period (Figure 4C). No significant (*p* < 0.05) correlation was found between GY of advanced lines and the average temperature (maximum, minimum or mean) for the wheat growing season (August–January) from 2006 to 2015.

**Figure 4.** Grain yield (**A**), plant height (**B**) and test weight (**C**) for 40–46 advanced lines and cultivars of durum wheat grown under full irrigation in Santa Rosa, from 2006 to 2015. Box and whisker plots show the population minimum, 25th percentile/median/75th percentile and maximum. The open symbols indicate outlier data and the closed symbols indicate the check cultivar 'Corcolén'.

#### **4. Discussion**

#### *4.1. Agronomic Traits*

Modern cultivars of spring durum wheat from Chile have a very high yield potential (~13 Mg ha−1) in a Mediterranean environment, under fully irrigated conditions. The yield potential achieved in Chile is clearly higher than values recorded in the Mediterranean basin. For example, high-yielding conditions in Spain usually do not surpass 8 Mg ha−1 [25,26], which is clearly lower than those achieved in the Mediterranean conditions of Chile. The high-yielding conditions in Spain usually imply several irrigations per season, particularly during the critical period from stem elongation to the middle grain filling, which alongside natural rainfall aims to balance the water lost due to accumulated evapotranspiration. Even so, the potential yields achieved in Spain are lower than in Chile due to a number of factors, such as Spain's shorter crop cycle duration, its higher night temperatures and the higher temperatures during the reproductive stage. The genetic advance in GY of spring durum wheat in the high-yielding environment of central Chile was 72.8 kg ha−1 per year (0.73% per year) for the period 1964–2010, and 63.4 kg ha−1 per year when the cultivar released in 1964 was excluded from the analysis (Figure 1). This is higher than the findings for spring bread wheat (43.5 kg ha−1 per year or 0.51% per year) for a similar period (1964–2008) and in the same Mediterranean environment [8]. It is also clearly higher than the increase reported for durum wheat in Spain (24 kg ha−1 y−1; 0.44% y−1) from 1980 to 2003, with no clear additional improvements occurring thereafter [26]. In northwest Mexico, under fully irrigated conditions, the genetic progress of spring durum and bread wheat varieties developed by CIMMYT was 0.49% and 0.41% per year, respectively, between 1966 and 2003 [27], and 0.88% per year when comparing eight bread wheat cultivars released between 1962 and 1988 [28]. A more recent study conducted at the same site in Mexico indicated that the GY progress was 30 kg ha−1 per year (0.59%) for spring bread wheat cultivars developed from 1966 to 2019 [29]. In Spain, under moderately irrigated conditions, the genetic progress of GY was 0.36% and 0.44% for Italian and Spanish cultivars of durum wheat, respectively, for cultivars released between 1920 and 2000 [9]. In South Australia, under rain fed conditions, the annual rate of increase in GY was 25 kg ha−1 for 13 cultivars released between 1958 and 2007 [30]. In North China, the annual genetic progress of spring bread wheat ranged from 0.48% (32.0 kg ha−1) for cultivars released between the 1960s and the 1990s [31], and in Henan Province values of 51.3 kg ha−1 per year have been reported for the last three decades [32].

The yield progress observed in advanced lines of the INIA-Chile breeding program (Experiment 2), which includes genetic and agronomic progress, has been much higher (569 kg ha−1 per year) than in all the studies discussed above. This large increase in GY is explained partly by the genetic progress, but overall the improvements have derived from better agronomic management of durum wheat in the central-south of Chile, and this has included modifications to irrigation and particularly adjustments in fertilization practices conducted during the first three years of the program. In winter bread wheat, the yield progress was 246 kg ha−1 per year (2.6%) between 1976 and 1998 in central Chile under fully irrigated conditions [14]. Clearly, fine tuning of crop management can have large impacts on GY in high-yielding environments when lines or cultivars of high yield potential are available.

Plant height was reduced from 107 cm in 1964 to an average of 90 cm for the period 1970–2010 (Figure 1A), and this was the consequence of the introduction of semi-dwarfing genes in Chile in the late 1950s [33]. Plant height was negatively correlated with the year of release in Australia, in cultivars developed between 1958 and 1973, but not in cultivars released after 1973 [30], and in China, in cultivars released between 1960 and 2000 [34]. A negative correlation between plant height and GY was also reported in the study of Zhou et al. [34]. However, the comparison of advanced lines produced during the last decade (Figure 2) showed a positive correlation between plant height and GY. These results suggest that plant height of semi-dwarf wheat below 70–80 cm may limit light interception and thus canopy photosynthesis and yield potential in high-yielding environments.

The HI increased between 1964 and 1970, but after that there were no changes. Furthermore, the correlation between GY and HI was not significant. The maximum values of harvest index (0.53) found in the current work were higher than those reported by Royo et al. [9] in a set of Italian and Spanish cultivars of durum wheat released between 1920 and 2000 and tested in Spain. In studies where cultivars released before and after the green revolution were evaluated, HI and the year of cultivar release were positively correlated (e.g., [9] for durum; [6,8,31] for bread wheat), but there was no correlation in cultivars released after 1970 (Figure 1C; see also [29]).

The increase in GY was positively associated with days to heading and KM2 and TKW (Table 3). The increase in the crop cycle in an irrigated Mediterranean environment contrasts with the breeding trend observed in rain fed Mediterranean areas, where early flowering, shorter duration cultivars are selected to escape post anthesis drought [9,35,36].

TKW increased significantly from 1964 to 1970, but the correlation with the year of cultivar release was not significant for the period 1970–2000 (Figure 1F). Genetic progress in TKW can be positive, negative or null depending on whether kernel weight has been a selection target for breeders and whether there have been changes in the number of grains per year (the trade-off between seed size and number in crops; [37]). For instance, in durum wheat growing in Mediterranean environments, kernel weight was superior in modern cultivars in Turkey [38], but remained unchanged in Italian and Spanish cultivars from the 20th century [9]. In bread wheat, kernel weight has been reduced [8,14,39] or has not changed [31] with genetic improvement.

#### *4.2. Grain Quality and Kernel Isotope Composition*

The test weight increased in modern cultivars and was positively correlated with TKW (*r* = 0.44; *p* < 0.05). The values of test weight obtained in this work are higher than those found in durum wheat genotypes grown under rain fed conditions in different zones of Spain [40,41]. Unfortunately, gluten and protein content did not improve between 1964 and 2010. Other studies comparing older cultivars or landraces to modern cultivars of durum wheat from Mediterranean countries have revealed lower grain nitrogen or protein content in the more modern cultivars [15,16,20,21,42]. In addition, the presence of *Rht* dwarfing genes in bread and durum wheat seems to reduce the concentration of Zn, Fe, Mn and Mg in kernels [43].

A number of studies have reported a negative correlation between grain protein concentration and GY in durum wheat [44] and in bread wheat [45,46]. It is probable that the lack of genetic progress in protein content is related to the strong increase in GY of the Chilean cultivars. However, this negative relationship should not be a limitation for genetic improvement in quality traits in grains of durum wheat because protein composition seems to be more important than the concentration [16,47].

The relationship between kernel δ13C and GY was negative, suggesting that genotypes exhibiting higher water use are the most productive [21,48]. In bread wheat under fully irrigated conditions, modern and more productive cultivars showed lower δ18O and δ13C, and higher stomatal conductance [8,18,21,25]. This negative relationship between δ13C and GY (or positive relationship between Δ13C and GY) has also been found in rain fed Mediterranean conditions ([21] for durum wheat; [49,50] for bread wheat), suggesting that the most productive lines are those able to maintain higher stomatal conductance and use more water [51]. In addition, the stomatal conductance of post green revolution wheat cultivars in Australia seem to show a lower sensitivity to vapor pressure deficit above 2 kPa compared to older cultivars [52], and this can be associated with lower (more negative) δ13C values.

In summary, changes in a number of traits have occurred in durum wheat cultivars selected for high-yielding environments in Chile. The large genetic progress in grain yield was associated with increases in days to heading, KM2 and TKW. The test weight has also increased with the year of cultivar release, but the gluten and protein content have not improved between 1964 and 2010. Interestingly, the increase in yield potential seems related to longer duration and somewhat taller plants that are able to use more water.

#### **5. Conclusions**

This study provided evidence that a high genetic advance in GY for durum wheat is feasible under high yielding conditions. The increase in GY was a consequence of a greater number of kernels per m2 and higher kernel weight in the more modern cultivars. The test weight was lower in the 1960s and increased curvilinearly with year of cultivar release. The gluten and protein content did not improve between 1964 and 2010. GY was negatively correlated with kernel δ13C, suggesting that genotypes exhibiting higher water use are the most productive. The yield progress of a set of advanced lines evaluated between 2006 and 2015 was very high, due to genetic progress, but this was also due to management improvements, particularly adjustment of fertilization practices conducted during the first three years. Unlike other Mediterranean agro-environments, a longer growing cycle together with taller plants seems to be related to the increase in the GY of Chilean durum wheat during recent decades.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4395/9/8/454/s1, Table S1: Monthly minimum (T min) and maximum (T max) temperatures and precipitation (PP) at Santa Rosa, Table S2: Mean values of cultivar traits according to the year of evaluation.

**Author Contributions:** I.M. designed the experiment and selected the germplasm. D.C. and K.R. were in charge of the management of the experiment and evaluation of agronomic traits. M.D.S. contributed to the isotope analysis. A.M.M.-E. contributed to data analysis. A.d.P. performed the data analysis and was in charge of writing the text but all the authors contributed to the manuscript.

**Funding:** This work was supported by the research CONICYT grant FONDECYT N◦ 1180252, Chile, and the contribution of Maria Dolores Serret was supported in part by the AGL2016-76527-R project from MINECO, Spain.

**Acknowledgments:** We thank Alejandro Castro for technical assistance with the field experiments.

**Conflicts of Interest:** The autors declare no conflict of interest.

#### **References**

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Durum Wheat Seminal Root Traits within Modern and Landrace Germplasm in Algeria**

## **Ridha Boudiar 1,2,\*, Juan M. González 3, Abdelhamid Mekhlouf 1, Ana M. Casas 4 and Ernesto Igartua 4,\***

Received: 11 March 2020; Accepted: 14 May 2020; Published: 16 May 2020

**Abstract:** Seminal roots are known to play an important role in crop performance, particularly under drought conditions. A set of 37 durum wheat cultivars and local landraces was screened for variation in architecture and size of seminal roots using a laboratory setting, with a filter paper method combined with image processing by SmartRoot software. Significant genetic variability was detected for all root and shoot traits assessed. Four rooting patterns were identified, with landraces showing overall steeper angle and higher root length, in comparison with cultivars, which presented a wider root angle and shorter root length. Some traits revealed trends dependent on the genotypes' year of release, like increased seminal root angle and reduced root size (length, surface, and volume) over time. We confirm the presence of a remarkable diversity of root traits in durum wheat whose relationship with adult root features and agronomic performance should be explored.

**Keywords:** proxy traits; genetic resources; root screening; root architecture

#### **1. Introduction**

The root system of wheat includes two main types, seminal (embryonic) and nodal roots, also known as the crown or adventitious roots [1,2]. Both types of roots play a crucial role in plant growth and are active throughout the whole plant life. Seminal roots, however, could be more important under specific circumstances, like drought conditions, as they penetrate deeper into the soil layers than nodal roots, making water in deep layers accessible to the plant [3–5]. Seminal roots also play a capital role during crop establishment, as they are the only roots existing before the emergence of the fourth leaf. Seminal roots include one primary root, two pairs of symmetric roots at each side, and, at times, a sixth central root [6].

The main features of root systems are encompassed under two categories, root system architecture (RSA) and morphology. RSA is related to the whole, or a large subset, of the root system, and may be described as topological or geometric measures of the root shape. Root morphology, as defined by J. Lynch, refers to "the surface features of a single root axis as an organ, including characteristics of the epidermis such as root hairs, root diameter, the root cap, the pattern of appearance of daughter roots, undulations of the root axis, and cortical senescence" [7]. The traits often used to describe wheat roots are total root length, root surface area, root volume, root angle, number of roots, and root diameter [8–10].

Roots are difficult to measure readily in natural conditions. Root trait determination has become accessible through the development of phenotyping methods in artificial systems, for instance, gel chambers [11], rolled germination paper [12], clear pots and growth pouches [13], "Termita" chamber and Whatman paper system [14], or growth pouches system [15]. Seminal roots can be phenotyped early and easily compared to the root system of mature plants [16,17], and for this reason, they have been proposed as good candidates to act as proxy traits in wheat [18] and maize [19,20]. Nevertheless, phenotyping these traits could be of interest only if they are useful to predict root growth and functioning in adult plants [21,22]. Indeed, several studies have found useful associations with traits in adult plants of wheat species [23–26]. For instance, the seminal root angle was correlated with nodal root angle [5,27], and with grain yield under drought conditions [28]. The seminal root number was correlated with thousand kernel weight (TKW) under stress, while the primary root length at the seedling stage was correlated with TKW under wetter conditions [25]. A steeper angle between the outermost roots and a higher root number in wheat seedlings have been linked to a more compact root system with more roots at depth in wheat [11,24,29].

Genotypic variation in root architecture has been reported within genotypes of different crop species [30–32], including wheat [13,21,25,33]. The presence of variation for the trait of interest is an essential requirement to improve the adaptability of crops under changing environmental conditions [34].

Local landraces are considered well adapted to the region where they were grown and contain large genetic diversity useful to improve crops like durum wheat [35]. These landraces were replaced by high yielding but more uniform semi-dwarf cultivars, better adapted to modern agriculture. However, scientists are convinced that local landraces still constitute a genetic resource useful to improve commercially valuable traits [36]. It is assumed that root traits enhance response to drought stress [37], but the realization of their contribution to superior grain yield depends on the type of drought and the agro-ecological conditions [38]. A deep rooting ideotype ("steep, cheap, and deep") was proposed by [39] to optimize water and N acquisition, building on the assumption that deeper rooting genotypes will use water that is beyond reach for shallower rooting genotypes. Modern breeding has caused some shifts in the root system architecture of durum wheat, from shallower and densely rooted systems in landraces of Mediterranean origin to deeper and more evenly distributed systems throughout the soil depth in cultivars worldwide [5].

The current study aims at evaluating the diversity of seminal root traits, including root angle and depth, during early growth of a set of durum wheat genotypes, consisting of modern cultivars and local landraces which are representative of the germplasm adapted to the mostly semi-arid conditions of Algerian cereal-growing regions before and after the advent of modern breeding. The study aims to reveal morphological diversity that could have agronomic relevance and, therefore, interest breeders.

#### **2. Materials and Methods**

#### *2.1. Plant Material*

We studied thirty-seven genotypes (landraces and modern cultivars), representative of durum wheat (*Triticum turgidum* ssp. *durum* Desf.) grown in Algeria. Geographical origins were varied (Algeria, France, Italy, Spain, Tunisia), and included genotypes produced at international breeding programs addressing semi-arid areas, namely the International Maize and Wheat Improvement Center (CYMMIT), the International Center for Agricultural Research in the Dry Area (ICARDA) and the Arab Center for the Studies of Arid zones and Drylands (ACSAD). These genotypes are representative of different periods of agriculture in Algeria, before and after the Green Revolution (Table 1).

**Table 1.** Name, type of cultivar, origin, and year of release of 37 genotypes of durum wheat used in the experiment.

\*: Backcross.

#### *2.2. Root Phenotyping*

#### 2.2.1. Preparation of Seeds

Twelve seeds of uniform size and healthy aspects were visually selected from each genotype and surface sterilized in a sodium hypochlorite solution (1.25% + one detergent drop, Mistol Henkel Iberica®). Seeds of each genotype were soaked and shaken in the solution for 15–20 min. Then, they were rinsed four times with sterile deionized water, in sterile conditions. Twelve seeds of each genotype were placed in Petri dishes, each with two filter papers soaked with 4 mL of sterile water. Then the Petri dishes were placed in a dark room at 4 ◦C for four days, and then at 22 ◦C/18 ◦C in a growth chamber with a 12 h light/darkness photoperiod for about 16 h.

Finally, the pre-germinated grains were transferred to the rhizo-slide system, described in detail in the next section and Figure S1. The experiment was carried out at the Laboratory of Cellular Biology and Genetics, Department of Biomedicine and Biotechnology of the University of Alcalá, Spain.

#### 2.2.2.  Sheets of A4-size black cardboard (180 g/m2, www.liderpapel.com) and filter papers were previously sterilized in an autoclave and then soaked in the nutritive solution Aniol [40]. The nutritive solution was prepared by dissolving 0.5550 g of Ca Cl2, 0.8215 g of KNO3, 0.6352 g of MgCl2·6H2O, 0.0165 g of (NH4)2SO, 0.0400 g of NH4NO3 in 100 mL of distilled water, to which 500 μL/L of Plant Preservative Mixture (PPMTM, Plant Cell Technology) at pH 5.8 was added. Each 8 mL was used to prepare 1 L of nutritive solution. Black cardboard with a nick made at the top center was placed on a glass plate with the same dimensions; then the pre-germinated grain (with embryonic part downward) was positioned just below the nick and covered by a filter paper. A black plastic sheet was used to cover the filter paper to ensure obscurity for roots, shifted ~2 cm upwards to allow better contact of the cardboard, and filter paper sheets with the nutritive solution. Two rhizo-slides were confronted to each other by the glass plate side, and the set was placed vertically in a glass box (internal dimensions of 32.2, 22, and 16 cm, length, width, and height) with two liters of the nutritive solution at the bottom, and then secured with two paper clips. Each glass box held 6 glass plates with two rhizoslides each, for a total of 12 seedlings, consisting of two genotypes, 6 seedlings for each (Figure S1). In total, each genotype was replicated 12 times. More details on the system are found in Ruiz et al. (2018) [25].

Once placed in the rhizoslides and the glass boxes, the seedlings were grown in a growth chamber for 7 days at 22/18 ◦C and 12/12 h photoperiod, day/night. The 37 genotypes were processed in batches of 6. Pre-germinated seeds of each 6 genotypes were placed into six glass boxes, each holding 6 seeds of two different genotypes. A complete batch comprised six boxes, three glass boxes prepared each Monday, and three each Thursday, every week. In total, 7 batches (14 runs) were performed until the experiment was completed (accounting for some seedlings that had to be replicated for various reasons). The set of genotypes for each run was selected randomly.

During the experiment, the boxes were replenished with distilled water every two days, to refill to the initial solution level. At the same time, to minimize seedling failure, each single seedling received 10 mL of the nutritive solution, applied with a pipette, near each seed. On the eighth day, the rhizo-slides were opened and shoots were immediately collected. The fresh roots were scanned using a Canon "LiDE210" scanner at 300 ppi to capture the first image then overlapped roots were manually separated and a second scan was done. The individual plant shoot dry weight (SDW) was obtained after oven-drying at 80 ◦C for six hours.

#### 2.2.3. Image Analysis

The two images of a rhizo-slide were analyzed using SmartRoot software v.3.32 [41] plugin for ImageJ1.46R (http://imagej.nih.gov/ij/download.html). The first image was used to measure only root angles and the second one to assess the other root traits using manual and semi-automatic SmartRoot procedures. Each root of the seedling was traced, semi-automatically, and then SmartRoot automatically generated the corresponding traits. In total, ten variables from the Smartroot output were recorded for each seedling: total root length (TRL), primary root length (PRL), mean length of the other seminal roots (MRL), total root surface area (Surface), mean root diameter (Diameter), total root volume (Volume), root number (RN), and shoot dry weight (SDW). The root angle was determined

for each root with respect to the vertical (90◦). From this determination, we extracted the maximum vertical angle (MVA) represented by the root growing with the steepest angle, the least vertical angle (LVA) represented by the root growing with the widest angle, and mean vertical angle (MRA) of all the roots, for each seedling.

#### *2.3. Statistical Analysis*

The experiment was considered a completely randomized design, with 12 replicates per genotype. Statistical analyses were performed using the REML (Restricted maximum likelihood) procedure with Genstat 18 [42]. Genotypes were considered as fixed factors and replications were considered as a random factor. The "Genotype" factor (n-1 degrees of freedom) was broken down into a single degree of freedom comparison of landraces vs. cultivars (named "Type" effect), and a "within type" factor (n-2) which corresponds to the variation of genotypes within each type. Multiple means separation was carried out using LSD at 0.05 level, for variables in which the F-value for "Genotypes" was significant. A principal component analysis (PCA) and a hierarchical cluster analysis (HC) were performed using the R package FactoMineR [43]. The hcut function was used for tree cutting levels truncation. The R package Factoextra [44] was employed for extracting and visualizing the results. Broad-sense heritability (h2) was calculated on an entry mean basis using the REML procedure, as follows:

$$\mathbf{h}^2 = \sigma^2 \sqrt{\left(\sigma^2 + \left(\sigma^2 \sqrt{\mathbf{r}}\right)\right)} \tag{1}$$

where σ2 g is the genotypic variance, σ2 e is the error variance and r is the number of replications.

#### **3. Results**

#### *3.1. Genotypic Variability*