Contribution of Wild Relatives to Durum Wheat ( Triticum turgidum subsp. durum ) Yield Stability across Contrasted Environments

: Durum wheat ( Triticum turgidum subsp. durum ) is mostly grown in Mediterranean type environments, characterized by unpredictable rainfall amounts and distribution, heat stress, and prevalence of major diseases and pests, all to be exacerbated with climate change. Pre-breeding efforts transgressing adaptive genes from wild relatives need to be strengthened to overcome these abiotic and biotic challenges. In this study, we evaluated the yield stability of 67 lines issued from interspeciﬁc crosses of Cham5 and Haurani with Triticum dicoccoides , T. agilopoides , T. urartu , and Aegilops speltoides , grown under 15 contrasting rainfed and irrigated environments in Morocco, and heat-prone conditions in Sudan. Yield stability was assessed using parametric (univariate (e.g., Bi, S 2 di, Pi etc) and multivariate (ASV, SIPC)) and non-parametric (Si1, Si2, Si3 and Si6) approaches. The combined analysis of variance showed the highly signiﬁcant effects of genotypes, environments, and genotype-by-environment interaction (GEI). The environments varied in yield (1370–6468 kg/ha), heritability (0.08–0.9), and in their contribution to the GEI. Several lines derived from the four wild parents combined productivity and stability, making them suitable for unpredictable climatic conditions. A signiﬁcant advantage in yield and stability was observed in Haurani derivatives compared to their recurrent parent. Furthermore, no yield penalty was observed in many of Cham5 derivatives; they had improved yield under unfavorable environments while maintaining the high yield potential from the recurrent parent (e.g., 142,026 and 142,074). It was found that a limited number of backcrosses can produce high yielding/stable germplasm while increasing diversity in a breeding pipeline. Comparing different stability approaches showed that some of them can be used interchangeably; others can be complementary to combine broad adaption with higher yield.


Introduction
Durum wheat (Triticum turgidum subsp. durum (Desf.)) is an important cereal cultivated worldwide with an annual production of 40 million tones [1]. Its importance worldwide is the result of its grain characteristics, which make it suitable to develop various products namely, pasta, couscous, and burghul among others [2]. Most of the area cultivated with durum wheat is in the Mediterranean region, accounting for 60% of global production [3].
Durum wheat is generally grown under the rainfed conditions of the semi-arid regions, where it is exposed to several biotic and abiotic stresses [1,4]. For instance, Hessian fly (Mayetiola destructor), a major pest for wheat in North America and the temperate Mediterranean drylands, can cause significant yield losses of more than 30% in Morocco [5,6]. Dis-variation across environments. In general, these genotypes will not respond to improvement in the environmental conditions, nor increase the yield in favorable environments. The "agronomic concept" defines a stable genotype as one with the minimum contribution to the GEI. According to the agronomic concept, stable genotypes will respond to change in the environments [48].
This research was conducted to study the contribution of durum wheat wild relatives to yield stability under different environments characterized by drought, heat, and disease pressure, and under optimal conditions. Different stability approaches were used to characterize both the germplasm stability and the relationships between the testing environments.

Plant Material
The germplasm tested here is composed of 67 lines of backcrossing populations derived from interspecific crosses of two durum wheat cultivars (Haurani and Cham5) with four wild wheat progenitors. 29 lines were derived from hybridization with the tetraploid progenitor (Triticum turgidum subsp. dicoccoides (syn. Triticum dicoccoides), and 47 lines from crosses with the three diploid ancestors Triticum monococcum subsp. aegilopoides (syn. Triticum aegilopoides), Triticum urartu and Aegilops speltoides. The choice of the recurrent parent was based on the local adaptation and drought tolerance of Haurani. Cham5, on the other hand, is a high yielding variety released in several countries from the ICARDA breeding program. The wild parents were selected based on their origins and the available information on their resistance to disease (mainly leaf rust). Table 1 provides a summary of the number of lines derived from each cross, with the number of backcrosses and a detailed list with DOIs given in Supplementary Table S1. The two recurrent parents and eight checks, including the released varieties and the ICARDA elite lines, were included in the trials and represented 13% of the total nursery. The IG refers to the accession number of the wild parent at the ICARDA genebank. The number following the Asterix (*) refers to the number of backcrosses.

Testing Environments and Experimental Design
The trials were conducted in 15 environments representing six locations during different seasons, between 2015 and 2018. Five locations were in Morocco, representing the Mediterranean hot and temperate environments, while Wad Medani in Sudan represented the hot and irrigated environment ( Table 2). At Tessaout and Melk Zher, two trials were Agronomy 2021, 11,1992 4 of 21 planted in the same season with different water regimes, one under full irrigation (FIR) and the second under rainfed (RFD) or supplemental irrigation (SIR). The purpose was to assess yield losses and the effects of late drought by comparing the two treatments/environments. At Tessaout, the fully irrigated trials received six irrigations, the first one at sowing, and the rest supplemented at different growth stages. The drought-stressed trials at Tessaout were irrigated only at sowing to ensure simultaneous germination with the irrigated trials. At Melk Zher, drip irrigation was used to supply a total of 411 mm for MZIR-16. MZRF-16 (the stressed environment) received 127 mm between rainfall and irrigation. The trials at Wad Medani were irrigated at an interval of seven days. All other trials were conducted under rainfed conditions, with the exception of Marchouch during the 2016 season, which received one supplementary irrigation during the vegetative stage. The trials were randomized in an incomplete block design (alpha-lattice) with two replications. Each replication was composed of eleven incomplete blocks, with seven plots in each block. Each plot consisted of four rows of two meters length, with a distance of 0.25-0.30 m between rows, and a sowing density of 300 seeds/m 2 . The recommended agronomic practices (land preparation, fertilizers, weeding, etc.) for each environment were applied. At maturity, the grain yield (GY) was estimated by harvesting and weighing the two internal rows, avoiding the borders, and then converting to kg/ha.

Analysis of Variance and Genotype by Environment Interaction
In order to investigate the genotype-by-environment interaction (GEI) and estimate the variance components, a linear mixed model and an Additive Main Effects and Multiplicative Interaction Model (AMMI) [54,55] were used for the analysis of variance. The mixed model was fitted using the sommer R package [60] in R version 4.0.4 [61] with the environment as fixed effects, and the genotypes and GEI as random effects. Diagonal variance structures were used to account for the heterogenous residual variance between environments, and therefore to estimate the residual variance in each environment. Similarly, the replication and block effects were considered heterogenous and estimated in each environment as random effects. The genotypic effects allowed the computation of the best linear unbiased prediction (BLUP) across environments for each genotype. The Metan R package [62] was used to perform the AMMI model and the results were used to compute stability parameters based on the interaction principal components from AMMI. The best linear unbiased estimations (BLUEs) in each environment were computed using Meta-R Agronomy 2021, 11,1992 5 of 21 software [63]. The BLUEs were used to run a genotype and genotype-by-environment interaction model (GGE) [64] using the GGEBiplots R package [65]. The GGE model was used to assess the representativeness and discrimination of each environment. Meta-R was also used to estimate the genetic correlation (ρ g ) between environments for grain yield, following Equation (1) from Cooper and DeLacy [66].
where, ρ gij is the genetic correlation, ρ pij is the phenotypic correlation between the environments i and j, h i and h j are the broad sense heritabilities in environments i and j, respectively.

Analysis of Stability
The BLUPs computed from the linear mixed model were the first stability parameter used for the ranking of genotypes. The stability indices derived from the AMMI model were sums of the absolute values of the IPCA scores [56] and AMMI stability values (ASV) [57]. The SIPC and ASV were computed according to Equations (2) and (3), respectively.
where P is the number of IPCs retained via the F-test, λ k is the eigenvalue of the k th IPC and α ik is the genotype principal component score.
where SSIPC1 and SSIPC2 are the sum of squares of the first and second IPC, respectively. IPC1 and IPC2 are the scores of the genotypes in the first and second IPC, respectively. The weighed average of absolute scores (WAAS) [59] was estimated following Equation (4). A superiority index (WAASY) was derived by rescaling the yield and WAAS to balance productivity and stability [59]. In the present study, yield and stability were given the same weight (50/50) (Equation (5)).
where, PCA ik is the score of the i th genotype in the k th IPCA, and EP k is the amount of variance explained by the k th IPCA where rG i and rW i are the rescaled values for GY and WAAS, θ Y and θ S are the weights for grain yield and stability assumed to be 50 for each in this study. The rest of the stability parameters were computed using the BLUEs from each environment. The Agrostab package [67] was used for the estimation of the regression coefficient (Bi), squared deviation from the regression (S 2 di) [45], environment variance (EV) [50] and the coefficient of variation [49]. The Metan package was used to compute Shukla stability variance (σ i 2 ) [47], geometric adaptability index (GAI) [68] and the superiority index (Pi) [51].
Four non-parametric stability indices [52,53] were also estimated using the Metan R package: S i1 , which is the mean of absolute rank difference over environment, S i2 , which is the variance of the ranks, S i3 , which is the sum of absolute deviations and S i6 , which is the relative sum of squares of rank for genotype. The Pearson correlation coefficients between the stability indices were computed using the Hmisc R package [69] and plotted using the corrplot R package [70].

Analysis of Variance
The analysis of variance from the linear mixed model and AMMI showed highly significant effects of the environment, genotypes, and their interaction (p < 0.001) (Tables 3 and 4). The highest proportion of variance was explained by the environment (66.93%), followed by the GEI (18.74%), while the genotypes accounted for 8.39% of the variance ( Table 4). The diagonal structure of the error in the mixed model validated the assumption of residual heterogeneity between the environments (Table 3). Therefore, accounting for the heterogeneity of error variance would increase the precision of genotypic variance estimation and thus the BLUPs across environments and their precision. The first seven interaction principal components (IPCs) from the AMMI model were significant, and explained 80.6% of the GEI ( Table 4). The variance explained by the first two IPCs was relatively low, accounting for only 38.8% of the GEI. The first IPC captured 22.6% of the variance while the second (IPC2) accounted for 16.2%, which highlights the complexity of the interaction patterns.

Climatic Data
The testing locations in Morocco represent typical Mediterranean semi-arid and temperate environments, while Wad Medani in Sudan represents dry hot irrigated environments. The maximum temperature at Wad Medani was consistently above 30 • C and no rainfall was registered during both cropping seasons. In Morocco, the rainfall distribution seemed to be as important as the total amount of rainfall in determining the type of environment. For instance, the Marchouch and Annoceur locations received almost double the amount of rainfall (514 and 611 mm, respectively) in 2018 compared to 2017 (Table 5). In terms of rainfall distribution, MCH-17 was exposed to a severe drought during the vegetative stage, while drought was more intense during the reproductive stage at MCH-17 and TSRF-17. Melk Zher was characterized by a dry season where the total rainfall registered was 85.5 mm, and drip irrigation was applied to differentiate between fully irrigated and drought-stressed environments. The trials in Melk Zher were irrigated, where MZIR-16 and MZRF-16 received, between irrigation and rainfall, a total of 411 mm and 127 mm, respectively. In terms of temperature, Annoceur and Marchouch were characterized by a cooler winter in comparison to other locations, while higher temperatures were observed during the reproductive stages in all locations.

Environment Characterization for Yield
The average yield across all environments was 3387 kg/ha. This yield is in the range of the national average yield in the Mediterranean region, which varies between 1500 kg/ha under rainfed conditions, and 4500 kg/ha under irrigated or favorable rainfed conditions. Six environments had yields above average, while nine were unfavorable and had lower yields ( Figure 1). The highest mean yield was registered at MCH-18 (6468 kg/ha), followed by TSIR-17 (5955 kg/ha) and TSIR-18 (5251 kg/ha). The lowest yield was observed under heat stress at WMD-17 (1370 kg/ha), followed by yields at AT-17 (1751 kg/ha), AN-17 (2149 kg/ha), and WMD-18 (2183 kg/ha). The yield at MCH-18 reflected the favorable season in comparison to MCH-17 where the average yield was 2533 kg. The supplemental irrigation at MCH-16 resulted in favorable conditions and a yield of 3993 kg/ha. The contrasts between fully irrigated and rainfed trials at Tessaout and Melk Zher were high, as shown by the yield reduction between irrigated and rainfed conditions ( Figure 1). On average, the yield at MZIR-16 was 73% higher than MZRF-16, while TSIR-17 had an average yield 164% higher than TSRF-17. The broad sense heritability ranged between 0.08 at AN-17 and 0.90 at WMD-18, and was moderate in the rest of the environments (Figure 1). Significant genotypic effects were observed in 12 environments and only three had non-significant genotypic effects (AN-17, AN-18, and WMD-17).
kg/ha under rainfed conditions, and 4500 kg/ha under irrigated or favorable rainfed conditions. Six environments had yields above average, while nine were unfavorable and had lower yields (Figure 1). The highest mean yield was registered at MCH-18 (6468 kg/ha), followed by TSIR-17 (5955 kg/ha) and TSIR-18 (5251 kg/ha). The lowest yield was observed under heat stress at WMD-17 (1370 kg/ha), followed by yields at AT-17 (1751 kg/ha), AN-17 (2149 kg/ha), and WMD-18 (2183 kg/ha). The yield at MCH-18 reflected the favorable season in comparison to MCH-17 where the average yield was 2533 kg. The supplemental irrigation at MCH-16 resulted in favorable conditions and a yield of 3993 kg/ha. The contrasts between fully irrigated and rainfed trials at Tessaout and Melk Zher were high, as shown by the yield reduction between irrigated and rainfed conditions (Figure 1). On average, the yield at MZIR-16 was 73% higher than MZRF-16, while TSIR-17 had an average yield 164% higher than TSRF-17. The broad sense heritability ranged between 0.08 at AN-17 and 0.90 at WMD-18, and was moderate in the rest of the environments (Figure 1). Significant genotypic effects were observed in 12 environments and only three had non-significant genotypic effects (AN-17, AN-18, and WMD-17).  lation results were corroborated by the GGE biplot for the association between environments ( Figure 3). The GGE biplot confirmed the high discrimination ability of MCH-18, TSIR-18, and MZIR-16. MCH-18 had the advantage of efficiently representing most of the other environments compared to TSIR-18 and MZIR-16. Low yielding environments such as WMD-17, AN-17 and AN-18 showed low discrimination ability for the genotypes (Figure 3).

Parametric Stability Indices
More than 50% of the tested lines had BLUPs for yield above average. The check Marzak had the highest yield (2854 kg/ha), while line 142003 had the lowest yield (1274 kg/ha). Four other checks (129080, Icarachaz, Louiza and Faraj), the recurrent parent Cham5, and four of its derivatives, were also among the highest yielding lines ( Table 6). The recurrent parent Haurani had the second lowest BLUP (1290 kg/ha), and the highest The color and size of the circles represent the direction and strength of the correlation between environments.

Parametric Stability Indices
More than 50% of the tested lines had BLUPs for yield above average. The check Marzak had the highest yield (2854 kg/ha), while line 142003 had the lowest yield (1274 kg/ha). Four other checks (129080, Icarachaz, Louiza and Faraj), the recurrent parent Cham5, and four of its derivatives, were also among the highest yielding lines ( Table 6). The recurrent parent Haurani had the second lowest BLUP (1290 kg/ha), and the highest yielding of its derivatives was 142001 (Haurani*2/T. urartu), with a yield 58% higher than its recurrent parent. The range of the BLUPs suggests that the germplasm tested here have high variation in yield. Therefore, the BLUPs should be taken into consideration for the interpretation of the other stability parameters in order to balance productivity and stability. Table 6. Ranking of the best ten stable and the least ten stable genotypes of durum wheat genotypes using parametric stability indices. According to the joint regression, the ideal genotype would have a regression coefficient (Bi~1) combined with low squared deviation (S 2 di) and high BLUPs. However, these three conditions could be met only for a few lines having average yield performance ( Figure 4). Some lines such as Icarachaz, 129080, and 142074 (Cham5*3/T. dicoccoides) combined higher yields with high Bi. In addition, Icarachaz was regarded as unstable according to S 2 di, as it had high variance (Table 6). Plotting Bi versus BLUPs identified an interesting group of lines which combined high productivity and high stability (Figure 4). This group included three checks (Marzak, Louiza and Faraj), the recurrent parent Cham5, and three of its derivatives. Marzak combined the highest yield (2952 kg/ha) with a regression coefficient of 0.96. The two derivatives 142009 and 142061 had yields of 2699 and 2597 kg/ha paired with regression coefficients of 1.09 and 1.05, respectively. They are derived from the same three backcrosses of Cham5 with T. aegilopoides. The derivative line 142005 (Cham5*4/Ae. speltoides) combined a yield of 2677 kg/ha with a Bi of 1.13. The landrace Haurani and its derivatives showed specific adaptation to unfavorable environments; they had low yields paired with low regression coefficients. The line 141972 (Haurani*2/T. urartu) had the lowest Bi (0.41) and was ranked 72 according to the BLUPs. Two Haurani derivative lines showed a significant yield improvement compared to their recurrent parent. One of these two derivatives (142001) was ranked the most stable according to its S 2 di. Since a lower S 2 di is associated with more predictable performance, adding the yield suggested that the lines 141995, 141966, 142071 and 142070 are desirable (Table 6, Figure 4). The S 2 di did not succeed in selecting genotypes combining high yields with broad adaptation. For instance, the lines 142003 and 142062 were ranked high according to the S 2 di, but their BLUPs and Bi suggested they were poorly adapted to all environments.

Desirable
Agronomy 2021, 11, x FOR PEER REVIEW 11 of 23 they had low yields paired with low regression coefficients. The line 141972 (Haurani*2/T. urartu) had the lowest Bi (0.41) and was ranked 72 according to the BLUPs. Two Haurani derivative lines showed a significant yield improvement compared to their recurrent parent. One of these two derivatives (142001) was ranked the most stable according to its S²di. Since a lower S²di is associated with more predictable performance, adding the yield suggested that the lines 141995, 141966, 142071 and 142070 are desirable ( Table 6, Figure 4). The S²di did not succeed in selecting genotypes combining high yields with broad adaptation. For instance, the lines 142003 and 142062 were ranked high according to the S²di, but their BLUPs and Bi suggested they were poorly adapted to all environments. Each point represents a single breeding line. Vertical and horizontal lines represents the grand mean ± 1 SD of BLUPs and Bi, respectively. Blue color represents the first ranked breeding lines according to Bi, red color represents the first ranked breeding lines according to the squared deviation from the regression (S²di). Some labels were repelled using black arrows to avoid overlapping.
The CV showed high efficiency to be used for negative selection to discard the low yielding unstable genotypes ( Figure 5). Some lines, such as 141979, 142062 and 142039, were among the lowest yielding, and had high CVs (67%, 65% and 60%, respectively). Marzak was ranked fifth with a CV of 42%, while the first ranked line 142015 (Cham5*2/T. urartu) combined a CV of 38% with a BLUP of 2316 kg/ha. The two lines 142001 and 141792 had yields below average, coupled with respective CVs of 40% and 41%. These two lines, derived from the same cross of Haurani and T. urartu with two backcrosses, were also selected as stable using the EV. However, the EV showed more affinity to select for the biological concept of stability than the CV. Eight of the best ten lines selected using EV had yields below average and only two were above average (142015 and 142055). High yielding lines like Icarachaz (3rd BLUP) and 129080 (2nd BLUP) were considered unstable, and were ranked 70 and 72 for EV, respectively (Table 6). Each point represents a single breeding line. Vertical and horizontal lines represents the grand mean ± 1 SD of BLUPs and Bi, respectively. Blue color represents the first ranked breeding lines according to Bi, red color represents the first ranked breeding lines according to the squared deviation from the regression (S 2 di). Some labels were repelled using black arrows to avoid overlapping.
The CV showed high efficiency to be used for negative selection to discard the low yielding unstable genotypes ( Figure 5). Some lines, such as 141979, 142062 and 142039, were among the lowest yielding, and had high CVs (67%, 65% and 60%, respectively). Marzak was ranked fifth with a CV of 42%, while the first ranked line 142015 (Cham5*2/T. urartu) combined a CV of 38% with a BLUP of 2316 kg/ha. The two lines 142001 and 141792 had yields below average, coupled with respective CVs of 40% and 41%. These two lines, derived from the same cross of Haurani and T. urartu with two backcrosses, were also selected as stable using the EV. However, the EV showed more affinity to select for the biological concept of stability than the CV. Eight of the best ten lines selected using EV had yields below average and only two were above average (142015 and 142055). High yielding lines like Icarachaz (3rd BLUP) and 129080 (2nd BLUP) were considered unstable, and were ranked 70 and 72 for EV, respectively (Table 6). Agronomy 2021, 11, x FOR PEER REVIEW 13 of 24 The selection intensity of Shukla stability variance (σi 2 ) was centered around the genotypes having an average yield performance coupled with low variance. The accession 141995 (Cham5*4/Ae. speltoides) was ranked first, followed by the lines 141966 and 142000, which are derived from Cham5 crossed to T. dicoccoides with two and three backcrosses, respectively. Shukla variance attributed lower stability to the genotypes with high and low yields such as Icarachaz and 141972 ( Table 6).
The superiority index (Pi) and the geometric adaptation index (GAI) selected similar genotypes for stability. Nine genotypes were included in the best ten lines using both GAI and Pi, among which were five checks, the recurrent parent Cham5 and four of its derivatives. The stable derivatives were issued from crosses with T. aegilopoides (142009 and 142060), T. dicoccoides (142074) and Aegilops speltoides (142005). Pi and GAI tended to select genotypes with high yield potential, as most of the low yielding genotypes were regarded as unstable (Table 6).
By giving the same weight to productivity and stability, the superiority index from the weighed average of absolute scores (WAASY) balanced productivity and general adaptation for the ranking of genotypes ( Table 6). The lines 141986 (Cham5*3/T. dicoccoides) and 142045 (Cham5*3/T. urartu) were the most stable, followed by Louiza and line 141966. WAASY also included the two highest yielding checks, Marzak and 129080, within the most stable genotypes. Interestingly, derivatives from crosses/backcrosses of Cham5 with the four wild species (Ae. speltoides, T. urartu, T. aegilopoides and T. dicoccoides) had superior stability/productivity compared to the recurrent parent. It was also noticeable that the same crosses resulted in genotypes with contrasting performance in yield and stability. This was the case for lines 141984 and 142074, both derived from Cham5*3/T. dicoccoides (Table 6).  Francis and Kannenberg (1978) coefficient of variation (CV). The red area highlights the undesirable genotypes while the green area highlights the stable and desirable genotypes. Vertical and horizontal dashed lines represent the mean of BLUPs and CV, respectively. Some labels were repelled using black arrows to avoid overlapping.
The selection intensity of Shukla stability variance (σi 2 ) was centered around the genotypes having an average yield performance coupled with low variance. The accession 141995 (Cham5*4/Ae. speltoides) was ranked first, followed by the lines 141966 and 142000, which are derived from Cham5 crossed to T. dicoccoides with two and three backcrosses, respectively. Shukla variance attributed lower stability to the genotypes with high and low yields such as Icarachaz and 141972 ( Table 6).
The superiority index (Pi) and the geometric adaptation index (GAI) selected similar genotypes for stability. Nine genotypes were included in the best ten lines using both GAI and Pi, among which were five checks, the recurrent parent Cham5 and four of its derivatives. The stable derivatives were issued from crosses with T. aegilopoides (142009 and 142060), T. dicoccoides (142074) and Aegilops speltoides (142005). Pi and GAI tended to select genotypes with high yield potential, as most of the low yielding genotypes were regarded as unstable (Table 6).
By giving the same weight to productivity and stability, the superiority index from the weighed average of absolute scores (WAASY) balanced productivity and general adaptation for the ranking of genotypes ( Table 6). The lines 141986 (Cham5*3/T. dicoccoides) and 142045 (Cham5*3/T. urartu) were the most stable, followed by Louiza and line 141966. WAASY also included the two highest yielding checks, Marzak and 129080, within the most stable genotypes. Interestingly, derivatives from crosses/backcrosses of Cham5 with the four wild species (Ae. speltoides, T. urartu, T. aegilopoides and T. dicoccoides) had superior stability/productivity compared to the recurrent parent. It was also noticeable that the same crosses resulted in genotypes with contrasting performance in yield and stability. This was the case for lines 141984 and 142074, both derived from Cham5*3/T. dicoccoides (Table 6).

Non-Parametric Stability Indices
Despite the differences in the ranking of genotypes between the four non-parametric stability parameters, some lines were identified as stable by all of them. S i1 , S i2 , S i3 and S i6 ranked the line 142053 (Cham5*2/T. dicoccoides) among the most stable, and the lines 141995 and 141966 were selected based on three of these parameters (S i2 , S i3 and S i6 ) ( Table 7). S i3 and S i6 selected similar lines for stability while S i6 had more affinity for the selection of high yielding genotypes (129080, 142053, Marzak). The lines selected by S i1 and S i2 could have low and high yields; this resulted in lower ability to differentiate superior from poorly adapted genotypes based on these indices (Table 7). Table 7. Ranking of the best ten stable and the least ten stable durum wheat genotypes identified using four non-parametric stability indices.

AMMI Derived Stability
The two stability parameters derived from the AMMI model selected average performance lines with low contribution to GEI as stable ( Figure 6). The most stable lines with low SIPC were 142008 (Cham5*3/T. dicoccoides), and 142071 and 142001, both derived from T. urartu crossed with Cham5 and Haurani. The lines 142045 (Cham5*3/T. urartu), 142063 (Cham5*3/T. dicoccoides) and 142032 (Cham5*3/T. aegilopoides) had the lowest contribution to GEI according to the ASV and were therefore the most stable. The inclusion of BLUPs with SIPC and ASV resulted in the identification of desirable lines with high yield potential ( Figure 6). The checks Louiza and the derivative lines 142074, 142009, 142060 can be recommended for combining stability with high productivity. ASV and SIPC showed some difference for the highest yielding lines, with Marzak and 129080 ranked differently by the two parameters.

Association of the Stability Indices
The correlation between stability parameters formed two major groups of positively correlated parameters and ranged from r 2 = −0.95 (BLUPs and Pi) to r 2 = 0.94 (BLUPs and GAI; EV and Bi). The first group was formed by BLUPs, WAASY, GAI, Bi and EV, with all indices positively correlated with the BLUPs. However, considering the interpretation of different parameters, the lines selected by EV will have low yields. The second group included the variance parameters (S 2 di and σ i 2 ), the two AMMI-derived stability indices, and the non-parametric index S i2 (Figure 7). ASV showed moderate correlation with the other parameters of the second cluster, ranging from 0.46 with S i2 to 0.76 with σ i 2 . As they have the same interpretation, the parameters in this group would select similar lines for stability. A strong positive correlation was found between Pi, S i3 and S i6 . These parameters were negatively associated with the second cluster and showed low to moderate correlation with the first group of parameters (Figure 7). WAASY was significantly correlated with all the other parameters except with S i1 , which itself did not correlate with any other parameter.
Pi, Superiority index; EV, Environment variance; GAI, Geometric adaptability index; σi 2 , Shukla variance; BLUPs, Best linear unbiased predictions; Bi, Regression coefficients; S 2 di, Squared deviation from the regression; CV, coefficient of variation; WAASY, superiority index from the weighed average of absolute scores; ASV, AMMI stability value; SIPC, sums of the absolute values of the IPCA scores, S i1 , Mean of absolute rank difference over environment; S i2 , Variance of the ranks; S i3 , Sum of absolute deviations; S i6 , Relative sum of squares of rank for genotype. . The green area represents the most desirable lines. Some labels were repelled using black arrows to avoid overlapping.

Association of the Stability Indices
The correlation between stability parameters formed two major groups of positively correlated parameters and ranged from r 2 = −0.95 (BLUPs and Pi) to r 2 = 0.94 (BLUPs and GAI; EV and Bi). The first group was formed by BLUPs, WAASY, GAI, Bi and EV, with all indices positively correlated with the BLUPs. However, considering the interpretation of different parameters, the lines selected by EV will have low yields. The second group included the variance parameters (S²di and σi²), the two AMMI-derived stability indices, and the non-parametric index Si2 (Figure 7). ASV showed moderate correlation with the were negatively associated with the second cluster and showed low to moderate correlation with the first group of parameters (Figure 7). WAASY was significantly correlated with all the other parameters except with Si1, which itself did not correlate with any other parameter.

Dissection of the Genotype by Environment Interaction
The analysis of variance showed the complexity of the GEI and the crossover interaction, meaning that genotypes responded differently to the changes in environment. High GEI is common in multi-environment trials [68,71,72], which reduces selection accuracy and genetic gains [18]. The dissection of the GEI and proper characterization of the environments and germplasm is therefore essential to improve genetic gains [73]. In this study, the genetic correlation allowed us to understand the relationships between pairs of environments and their interaction. The small contribution of Marchouch to the GEI showed its suitability for efficient selection of superior genotypes for other environments. Despite the climatic differences between years, the same superior genotypes can be se-

Dissection of the Genotype by Environment Interaction
The analysis of variance showed the complexity of the GEI and the crossover interaction, meaning that genotypes responded differently to the changes in environment. High GEI is common in multi-environment trials [68,71,72], which reduces selection accuracy and genetic gains [18]. The dissection of the GEI and proper characterization of the environments and germplasm is therefore essential to improve genetic gains [73]. In this study, the genetic correlation allowed us to understand the relationships between pairs of environments and their interaction. The small contribution of Marchouch to the GEI showed its suitability for efficient selection of superior genotypes for other environments. Despite the climatic differences between years, the same superior genotypes can be selected for the three seasons at Marchouch, and also for other locations. The crossover interaction due to year effects was more pronounced in the other locations; these effects were mainly associated with the total and distribution of rainfall and temperature for rainfed trials. These findings confirm previously reported effects of weather conditions on the GEI in durum wheat in the Mediterranean basin [3,74,75]. Annoceur was characterized by low heritability during both seasons, and it was among the locations where low genetic gain for durum wheat was achieved [76]. The testing environments exposed durum wheat to favorable conditions, drought, and heat stresses, in addition to disease pressure, mainly leaf rust at Allal Tazi and tan spot at Annoceur. In addition, the contrasting environments at Tessaout and Melk Zher showed their effectiveness to select for drought tolerance [77]. The genetic correlation showed that some genotypes can be selected simultaneously for drought and heat tolerance. This approach of using different environments was reported to be useful to select for productivity and tolerance to major biotic and abiotic stresses simultaneously [42].

Considerations for the Use of Durum Wheat Wild Relatives
Much of the progress in durum wheat breeding and genetic gains was achieved at the cost of a reduction in genetic diversity [76]. Maintaining the level of genetic gains requires a diverse gene pool for the breeders to select for adaptive traits [78]. The use of crop wild relatives (CWR), therefore, is strategic as it allows several traits to be increase simultaneously, while recovering lost diversity [79]. In this context, the present study investigated the potential contribution of wild relatives to yield stability in durum wheat. The choice of wild parents from the primary gene pool was of high importance, as several barriers are associated with the use of species from the secondary and tertiary gene pool [80,81]. The primary gene pool allows a higher frequency of recombination, which can be useful for quantitative traits [82]. Interestingly, the performance of the derivative lines from the four wild parents was not affected by the number of backcrosses, which suggests the possibility of maintaining a certain level of diversity while using these species by reducing the backcrossing. Instead, more care should be given to the selection scheme and intensity to reduce the frequency of undesirable linkages. Our results showed that the same crosses with T. dicoccoides and T. aegilopoides produced lines outperforming the recurrent parents, but also lines with very low yields. El Haddad et al. [83] reported that no linkage drag on agronomic performance was associated with the use of wild relatives in durum wheat. In the case of bread wheat, the adoption of an appropriate selection strategy resulted in the elimination of unfavorable traits and the release of high-yielding and stable varieties from synthetic hexaploid-derived lines [84]. This is where pre-breeding becomes crucial in the process of gene introgression from CWR. Pre-breeding should retain favorable alleles while returning the background of the elite parent through reasonable top crossing [78].

Impact of Wild Relatives on Yield Potential
The use of two contrasting recurrent parents, the landrace Haurani with low yield potential and high tolerance to drought, and Cham 5, a high yielding cultivar released in several countries, was useful. Haurani derivatives showed an important advantage in both yield and stability, and they overpassed their recurrent parent in most of the environments. For instance, line 142001, derived from T. Urartu, had significantly higher BLUPs associated with higher stability, according to both Bi and S 2 di. The same line showed a yield increase under drought stress, combined with earliness and cooler canopy under heat stress [77,85]. Under favorable conditions, some Haurani derivatives (142064 for example) can yield twice as much as the recurrent parent, which indicates a high contribution to yield potential. In the case of Cham5, which is high yielding, the expected contribution of wild relatives can be more pronounced under stressed environments. The top yielding lines under drought stress at TSRF-17 were Cham5 derivatives crossed with the four wild parents. In fact, line 142026 (Cham5*3/T. urartu) outyielded Cham5 under optimal conditions (TSIR-17) and was subject to lower yield losses due to drought. Similarly to our results, a significant increase in yield and stability was also reported using T. dicoccoides and T. monococcum under drought and terminal stress [42]. The value of CWR in Cham5 derivatives was even greater under heat stress during both seasons. The contribution to heat tolerance was observed at the level of yield and its components (T. dicoccoides and T. urartu), phenology (T. urartu) and physiological response (T. aegilopoides, T. urartu, and T. dicoccoides) [77,[85][86][87]. The results from MCH-18, which was the most favorable environment, were interesting. Cham5 had the highest yield, followed by four of its derivatives, with yields above 9 tones/ha. This finding showed that the use of CWR from the primary gene pool does not always come with a penalty on yield. These findings confirm that the increase in yield from crosses with wild relatives is mainly attributed, but not restricted, to the improvement of resistance/tolerance to biotic and abiotic stresses [25,26].

Yield Stability of the Durum Wheat Derivatives
Our study used different stability approaches to select stable genotypes for both biological (static) and agronomic (dynamic) stability. It was not surprising that most of the checks had high agronomic stability, especially that Cham1 (syn. 129080), Cham5 (syn. 129081), Louiza and Marzak had wide adaptation combined with moderate to high yield potential. When balancing productivity and stability, several lines derived from crosses with T. aegilopoides, T. dicoccoides, T. urartu and Ae. speltoides were identified. For example, lines 141995 (Cham5*4/Ae. speltoides) and 141966 (Cham5*2/T. dicoccoides) had yields above average, lower variance (σ i 2 ) and deviation from regression (S 2 di). In addition, they were highly stable according to SIPC, WAASY and the non-parametric indices (S i2 , S i3 and S i6 ). These lines are desirable under unpredictable weather conditions, as they can maintain average performance during poor seasons and respond positively to favorable conditions. The line 142074 (Cham5*3/T. dicoccoides) can also be recommended for its dynamic stability, as it was highly ranked for Pi, BLUPs, GAI, S i6 and had average stability for WAASY. These findings are in line with the reported contribution of T. dicoccoides and Ae. speltoides to yield stability while conserving high yield potential [83,88]. Simmonds et al. [89] reported the release of a bread wheat variety derived from a cross with T. dicoccoides, which exhibited high stability under different environments and cultural practices. The contribution of CWR to yield stability can be the result of the simultaneous improvement of yield and its components under unfavorable environments [90][91][92][93][94][95][96]. Our results showed that the use of wild relatives can supply lines with wide adaptation, and lines specifically adapted to unfavorable or favorable environments.

Association between Different Stability Parameters
Selection from multi-environment trials is an important component for plant breeders, and the adoption of appropriate statistical analysis is crucial to improve the selection accuracy. The use of a linear mixed model provided a significant advantage as it accounted for the heterogeneity of the variance between environments. The advantage of mixed models in increasing the prediction accuracy has been reported in other studies [59,97]. Interestingly, BLUPs did not correlate significantly with the SIPC and ASV from the AMMI model. The reason behind this is that these indices select lines with low contribution to GEI, which will favor lines with average yield performance. This also explains why SIPC and ASV were clustered with the S 2 di and σ i 2 and confirms the positive association previously reported by Sneller et al. [56]. The negative correlation between CV and BLUPs was useful, as it allowed the elimination of the less stable lines combining low yield potential with high variation. The value of the regression coefficient to select for agronomic stability was confirmed by its correlation with Pi, BLUPs and GAI. Different and sometimes opposite correlations between the stability statistics have been reported by other studies [21,68,[98][99][100]. These correlations varied based on the population structures and environmental characteristics. In this view, flexible stability statistics such as the superiority index WAASY can be more useful. The weights of productivity and stability can be adjusted depending the yield, size of the population, and the environment to select superior genotypes.

Conclusions
The backcrossing population of durum wheat provided a view on the performance of different lines derived from similar crosses with regard to linkage drag. Many lines can be selected for the unpredictable climatic conditions of the Mediterranean region. These lines combine high productivity and stability. It was concluded that it is possible to develop diverse adapted high-yielding germplasm using CWR. This can be achieved by deploying the proper selection schemes at the pre-breeding level. WAASY can be recommended as a flexible stability index to select for stable and highly productive germplasm.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/agronomy11101992/s1, Table S1: List of accessions and checks used in the trials with accession number, DOIs and pedigree.