Investigating Stability Parameters for Agronomic and Quality Traits of Durum Wheat Grown under Mediterranean Conditions

: Durum wheat in the Mediterranean grows under rainfed conditions, where unpredictable climatic conditions result in substantial variation in grain yield and quality. Climate change intensiﬁes Genotype × Environment interactions and urges breeders to escalate their efforts to breed cultivars combining high performance and stability. The current study aimed to appraise the relations between twelve stability parameters derived by different statistical models for yield, yield-related and quality traits of durum wheat grown under Mediterranean conditions. Stability parameters were estimated in two experiments of twenty and sixteen cultivars, respectively. The parameters were categorized into three groups. Group A included Additive Main Effect and Multiplicative Interaction (AMMI)-derived parameters (ASV and AWAI), Wrickle’s ecovalence (Wi), Shukla’s stability variance ( σ 2 ), and the nonparametric parameters Si (1) and Si (2) . Group B included regression parameters (bi, Bi_A), Coefﬁcient of Variance (CV), and Superiority measure (Pi). Group C encompassed deviation from regression parameters (s 2 di-DJi) when the heterogeneity of the slope was signiﬁcant. Correlations between stability parameters for different traits and the between stability parameters and the traits per se were modest. Stability parameters of Group B had higher repeatability for grain yield. The results of the present study contribute to the adjustment of durum wheat breeding strategies.


Introduction
Durum wheat is a major crop for the Mediterranean basin [1]. This area is also the most significant import market and the largest consumer of durum wheat products [2,3]. The crop has been closely linked with regional diet and tradition. Further to pasta production, durum wheat is widely used for the preparation of bread, couscous, bourghul, and other traditional products [4].
The primary goal of all breeding programs has been the increase of grain yield. Studies have shown that genetic gains in grain yield are mainly achieved by increasing harvest index and number of seeds/m 2 , the latter trait being determined by number of spikes/m 2 and seeds/spike [5,6]. Adjusting the major components of plant phenology, that is plant height and days to anthesis (or heading), according to the targeting environment has been also crucial [2,7]. Moreover, breeding efforts target quality traits that are desirable for the industry, such as yellow pigment concentration, protein content, gluten quality and volume weight [8][9][10].
Durum wheat in the Mediterranean basin usually grows under rainfed conditions. Sowing occurs on November-December and the plants complete their cycle by the end of spring. Plants often experience drought and heat stress events during their development resulting in yield loss and, to lesser extent, affecting quality [9][10][11]. The Mediterranean basin is severely affected by climate change [12], increasing yield uncertainty, deteriorating quality, and impairing genetic gains [1]. On the other hand, breeders need to intensify their efforts to meet the demand arising from the growing population and the quality requirements of the industry. Thus, the ultimate goal of plant breeders is to breed cultivars combining high performance and stability for agronomic and quality traits. However, major challenges arise from the significant Genotype × Environment interactions (G × E) resulting from the differential response of genotypes at different environments [13,14]. The large proportion of G × E, which is usually associated with cross over interactions, impairs breeders' goal for high yielding and stable cultivars across environments [15,16]. This goal is especially cumbersome in the Mediterranean area due to high annual variation in grain yield and quality traits resulting from the unpredictable climate conditions [1,9,10].
A number of statistical models have been developed for analyzing G × E interactions [13]. The genotypic specific response in different environments (norm of reaction) can be captured by stability parameters [7,16,17]). The most widely used model for analyzing G × E is the linear model with the regression of individual genotypes' mean values over environmental means [18]. The slope of the regression of individual genotypes and their deviation from regression can be used as stability parameters [19,20]. Later, bilinear models have been introduced which are more efficient to overcome the limitations of the linear regression models [21] because they can capture a higher percentage of the variability in G × E interactions [13,16]. The most widely used bilinear model is the Additive Main Effect and Multiplicative Interaction (AMMI) model [22]. Stability parameters derived from AMMI models have been proposed [23], among them AMMI Stability Value (ASV) [24] and AMMI wide adaptation Index (AWAI) [25]. Several other parametric stability parameters have been introduced including the Wrickle's ecovalence (Wi) [26], Shukla's stability variance (σ 2 ) [27], Coefficient of Variance (CV) [28], and Superiority measure (Pi) [29]. Stability parameters based on non-parametric procedures have also been proposed [30,31] having advantages over parametric procedures, including that no assumptions about the distribution of the data are needed and that they are less sensitive to bias caused by outliers [13,17].
Lin et al. [32] categorized stability parameters into three types. Type I includes genotypes with small variance among environments and corresponds to the static (biological) stability. Type II refers to genotypes of which their response to environments is parallel to the mean response of all genotypes and corresponds to the dynamic (agronomic) stability. Type III is related to genotypes' residual mean square from the linear regression models. Stability parameters related to the regression coefficients on linear models can correspond to the biological or the agronomic concept depending on how a stable genotype is defined. Agronomic stability is more relevant to breeders' objectives for yield, nevertheless, biological stability may deserve consideration for quality traits [15].
Most of the published literature focuses on stability for grain yield. However, stability for yield, yield components, phenotypic (heading data and plant height), and quality traits may be related and may interact with the traits per se [33][34][35][36][37]. Significant correlations between stability parameters related to different stability concepts or between stability parameters and the traits per se implies the existence of trade-offs that need to be considered by plant breeders [38]. Therefore, the relations between agronomic and quality traits and their stability deserves further investigation.
The calculation of stability parameters depends on the set of genotypes and environments included in the analysis [21,32,39]. Furthermore, stability for wheat yield, heading date, and plant height might have more complex architecture than the traits per se [40]. The complex genetic structure of traits' stability implies that their expression can be substantially affected by the environment and might exhibit lower repeatability and heritability than the traits per se. Such attributes can limit the usefulness of stability parameters in plant breeding programs [41]. Despite its importance for plant breeding, the repeatability of the stability parameters has not been exhaustively studied [39,42,43]. The current study aimed to appraise the relations between twelve stability parameters derived by different statistical models for agronomic and quality traits of durum wheat grown under Mediterranean conditions. Stability parameters were estimated in two experiments of twenty and sixteen modern cultivars, respectively, with diverse origin. The correlations between stability parameters for different traits and the traits per se were investigated to shade light on plausible implications in their use and existence of trade-offs that need to be considered by plant breeders. The repeatability of the stability parameters has been marginally tackled by the literature. Therefore, the relations between stability parameters derived from eight cultivars that were common to the two experiments were investigated as a measure of their repeatability.

Plant Material
In Experiment One, twenty durum wheat cultivars (Triticum turgidum subsp. durum) of diverse origin, released for cultivation over the last 40 years, were selected (Table 1). These cultivars were found to be genetically diverse in previous studies applying molecular markers [44,45]. Six cultivars were bred by the Cypriot National Breeding Program and represent the main commercial cultivars grown in Cyprus over the last forty years. The other fourteen cultivars were released by other breeding programs, targeting areas with similar climatic conditions. 'Ourania' and 'Hekabe' are the cultivars currently grown in Cyprus combining high yield and good quality traits [46]. 'Omrabi5', 'Waha', and 'Korifla' are among the most widely grown cultivars released by ICARDA and are usually used as checks to international nurseries [25]. 'Simeto', 'Claudio', 'Iride', and 'Svevo' are appreciated for their quality traits, e.g., kernel weight, volume weight, gluten quality and protein content, respectively [47]. 'Matt', 'Mexicali81 , and 'Pisti' were selected because of their early heading which is one of the most important adaptation traits under Mediterranean environments [2].  The calculation of stability indices depends on the set of genotypes and environments included in the analysis. Therefore, a different set comprised of sixteen durum wheat cultivars were selected in Experiment Two, out of which eight were common to the two experiments (Table 1). Seven cultivars were bred by the Cypriot National Breeding Program including the two most widely grown cultivars in Cyprus 'Ourania', and 'Hekabe' and four promising lines. Five widely grown cultivars were kindly provided by the ICARDA. Based on the results of Experiment One, 'Pisti' and 'Matt' were included as early heading while 'Simeto' and 'Claudio' as late heading cultivars. The later ones were also included for their high kernel and volume weight, respectively.

Field Experimental Conditions and Design
The first field experiment (Experiment One) was conducted in three locations (Achelia-34 •  . The three experimental locations exhibit diverse edaphoclimatic conditions (Annex S1). Achelia has higher temperatures during winter and lower temperatures during spring, higher precipitation and deep clay soil favoring high yields. Athalassa has shallow sandy clay loam soil and lower precipitation during crop cycle resulting in drought stress during heading and grain filling. In addition, the higher day temperatures in spring and the frequent occurrence of extreme high temperatures during grain filling, very often result in heat stress. Dromolaxia also has sandy clay loam soil and it depicts intermediate climatic conditions and intermediate yield. Durum wheat is not cultivated in Athallasa because of the high probability of crop failure due to the adverse conditions. However, it was selected to represent the climatic conditions that would likely prevail in the current durum wheat Mediterranean zone in the future due to climate change [12].
The design of Experiment One was a randomized complete block with four replications. Each plot was 8 m long, and comprised of six rows, spaced apart 0.175 m. Seed rate was adjusted to 226 germinating seeds/m 2 . Experiments were set up at the end of November in Athalassa, and on December in Dromolaxia and Achelia. Planting always occurred after fallow. The fields were fertilized before sowing with 60 Kg/ha of N 2 and 60 Kg/ha of P 2 O 5 , respectively. In the 1st growing season, experiments in Achelia and Dromolaxia were additionally top dressed with 40 Kg/ha of N 2 at tillering. Weeds were chemically controlled at tillering (Atlantis ® Bayer, Illoxan ® Bayer, Granstar ® DuPont). During the 1st growing season, Valiant ® (Agriphar) was sprayed to control high infestation of cereal tortricid (Cnephasia pumicana) in Achelia. Infestation from diseases was negligible. Irrigation was applied in Athalassa during booting (30 mm) in the 1st growing season, and during tillering (50 mm) and booting (50 mm) in the 2nd growing season. Irrigation was also applied in Achelia during anthesis (60 mm) in the 2nd growing season.
The second experiment (Experiment Two) was conducted at three locations, Achelia, Dromolaxia and Polis (35 • 01 N, 32 • 25 E) for three consecutive growing seasons (2014/15, 2015/16, 2016/17). The three experimental sites represent the main durum wheat-producing areas of Cyprus. Polis and Achelia favor higher yields compared to Dromolaxia because of their higher precipitation and the lower temperatures during grain filling (Annex S1).
The design of Experiment Two was as described in Experiment One with two modifications. Plots were 6 m long and the number of replications was six. Experiments were set up in December, except from Polis at the growing season 2016/17 and Dromolaxia at the growing season 2014/15 where experiments were set up in January and November, respectively. In Achelia and Polis, the fields were fertilized before sowing with 60  ically controlled at tillering (Atlantis ® Bayer, Biopower ® Bayer, Denok amine ® Nufarm S.A.S) according to the manufacturer's instructions. Infestation from pest and diseases were negligible.

Data Collection for Agronomic and Quality Traits
For both experiments, heading date was recorded when the ears of 50% of the tillers had emerged from the flag leaf sheaths for approximately half their length and was expressed as growing degree days from emergence to heading (GDD). Plant height (PH) was measured at physiological maturity as an average of three measurements/plot, excluding awns. The plots were mechanically harvested and grain yield (GRYLD) was adjusted to 12% moisture level. Thousand kernel weight (TKW) was calculated as the mean weight of two samples of 200 seeds/plot. Volume weight (VW) was measured with a 0.5 L chondrometer (Seedburo) and expressed as Kg hl −1 .
In Experiment One, the number of fertile tillers/m 2 (NTLSM) was counted from the second and fifth rows in order to eliminate the border effect. Plants from these rows were hand harvested to estimate the number of seeds/spike (SPS). Sampling area for determining NTLSM and SPPS was 2 m long at each row. Number of seeds/m 2 (NSSM) was calculated using the formula "NSSM = NTLSM × SPS". In Experiment Two, NSSM was calculated according to the formula NSSM = (GRYLD (Kg/ha) × 100)/TKW. The two different methods applied for the calculation of NSSM provided similar results (data not shown).
Protein content (PRO), yellow pigment content (CAR), and gluten index (GI) were assessed on whole grain flour samples, only in Experiment One. Grain nitrogen content was determined according to Kjeldahl method. PRO was calculated multiplying the n value by 5.7 and expressed as a percentage on a dry weight base. CAR was determined as described by [48] and expressed in ppm. GI was measured according to [49], using a Perten Glutomatic (Perten Hägersten Sweden) with a minimum of two repeats per sample. PRO, CAR, and GI were assessed on two samples per cultivar. Each sample was collected from the bulk of two replications.

Statistical Analysis
Combined analysis of variance was conducted for all traits considering cultivars and environments as fixed factors using META-R software [50]. Means were compared by applying the LSD test at significance level 0.05. Broad sense heritability was calculated according to the formula H 2 = σ 2 g/[σ 2 g + (σ 2 ge/nLoc) + σ 2 e/(nLoc × nRep)], where σ 2 g, σ 2 e and σ 2 ge are the genotype, error, and G × E interaction variance components, respectively, nRep is the number of replicates and nLoc is the number of environments. Box plots were constructed to depict the trait variation in each environment. Additive Main Effect and Multiplicative Interactions (AMMI) analysis was performed by GEA-R [51]. The F-test proposed by [52] was used for determining the number of terms to be retained in a multiplicative model. The following stability parameters were calculated using GEA-R; Coefficient of Variance (CV) [28], Wrickle's ecovalence (Wi) [26], Shukla's stability variance (σ 2 ) [27], the slope of the regression of individual genotypes, and their deviation from regression bi and s 2 di according to Eberhart and Russell [19] and Bi_A and DJi according to Perkins and Jinks [20], Superiority measure (Pi) [29] and nonparametric coefficients Si (1) and Si (2) [30,31]. The former parameter is based on the mean of the absolute rank differences of a genotype over all tested environments and the latter on the variance among the ranks over all tested environments. In addition, two stability parameters derived from AMMI analysis were calculated; AMMI Stability Value (ASV) [24] and AMMI wide adaptation Index (AWAI) [25]. ASV is genotypes' stability considering its distance from IPCA1 and IPCA2 axes. AWAI is genotypes' stability as determined by its distance from each significant IPCs axis. Hierarchical cluster analysis was conducted to assess the relations between stability parameters for each trait. Squared Euclidean distances were calculated on standardized Z values, with a mean of 0 and a standard deviation of 1. Clustering was performed using the 'WARD' method. Spearman correlations were computed to depict relations between stability parameters and between stability parameters and the traits per se. For the assessment of parameters' repeatability, Spearman correlations between parameters derived from eight cultivars that were common to the two field experiments were calculated. Box plots construction, hierarchical cluster analysis and estimation of Spearman correlations were carried out using SPSS (IBM, SPSS ver. 26). Heat maps on Spearman correlations were constructed by Morpheus online software (https://software.broadinstitute.org/morpheus) [53].

Results
Average maximum and minimum temperatures were higher than normal for most months during growing seasons with the exception of the growing season 2011/12 (Annex S1). In some cases, temperatures were up to 2 • C higher than normal. Precipitation was exceptionally high during the growing season 2011/12, substantially lower than normal during the growing seasons 2015/16 and 2016/17, and close to average the other seasons. There were months with extremely higher or lower precipitation than normal. For example, February precipitation was negligible during the growing season 2016/17.
There were significant differences between cultivars for all traits for both experiments (Tables 2 and 3). The differences were significant in all environments with the exception of GRYLD in Athalassa and for PH in Dromolaxia during growing seasons 2013/14 and 2015/16, respectively. The percentage of variance explained by cultivar was low, particularly for GRYLD for both experiments. The yield component with the highest contribution to the genotypic variance was SPS, which had the highest heritability. The lower heritability was recorded for GRYLD for both experiments. The highest heritability and percentage of variance explained by genotypic effect was recorded for CAR followed by GI.    The differences between environments were significant for all traits for both experiments (Tables 2 and 3). The high percentage of variance explained by environments, particularly for GRYLD reflects the diverse climatic conditions during the experimentation (Annex S1). The favorable climatic conditions in Achelia during season 2011/12 resulted in high yield (7192 Kg/ha) while yield was below 2000 Kg/ha at the location Athalassa ( Figure 1). The adverse climatic conditions in Athalasa were particularly evident during grain filling. Two of the major characteristics appreciated by the industry, TKW and VW, were substantially lower in Athalassa. TKW was 23,56 g in 2012/13 and 25,46 g in 2014/15 and the respective VW was 66.55 and 69.31 (Kg/hl). On the other hand, CAR and PRO were substantially higher in Athalassa.
For Experiment Two, the highest yield was obtained in Achelia during the season 2015/16 (5982 Kg/ha) and the lower in Dromolaxia during seasons 2015/16 and 2016/17 (689 and 997 Kg/ha, respectively). Figures 1 and 2 show the genetic variation and the variation between environments for all traits. Some cultivars manage to sustain high or low performance for particular traits and appear on the box plots as outliers. For example, CAR was consistently low to all environments for cultivar 'Aronas' contrary to cultivar 'Matt'. Cultivars 'Iride' and 'Adnan2' retained their high SPS in most environments. The superiority of cultivar 'Simeto' for thousand kernel weight was more evident under the adverse conditions during grain filling at the location Athalassa.  Cultivar × Environment interactions were significant for all traits (Tables 2 and 3) with the percentage of variance explained ranging from 1.94% for GRYLD to 14.1% for GI in Experiment One and from 1.34% for GRYLD to 17.1% for VW in Experiment Two. Among yield components, G × E interactions were larger for SPS (7.08%) in Experiment One and for TKW (5.48%) in Experiment Two. In most cases, two or more PC axes from AMMI analysis were significant according to the Gollob test. The biplots of AMMI analysis are shown in Annex S2. The first two axes explained 74.81% and 53.26% of the variability for GRYLD, in Experiments One and Two, respectively. The percentage of variance explained by the heterogeneity of the slopes from the linear regression model was lower than that of the AMMI model. The heterogeneity of the slopes explained 11.23% and 21.17% of the variability for GRYLD in Experiments One and Two, respectively. The linear model was more efficient for explaining higher percentage of variability for quality traits (54.07% for CAR and 55.99% for GI). Significant differences between the slopes were detected to all traits except GRYLD, NTLSM, NSSM and SPS in Experiment One, and except TKW and VW in Experiment Two.
Figures 3-5 depict the relations between stability parameters for each trait. There were consistent high relations between Bi_A and bi, DJi and s 2 di, and Wi and σ 2 . The two stability parameters derived by the AMMI model (ASV and AWAI) consistently clustered together. This was also the case for the two nonparametric parameters (Si (1) and Si (2) ). With few exceptions, Wi-σ 2 and s 2 di-DJi grouped together in cluster analysis. Pi was the parameter which consistently had strong negative correlations with trait mean values. Significant negative correlations between CV and traits were observed in some cases and to a lesser extent positive correlations between bi-Bi_A and the traits, with the exception of VW where the correlation was negative. The correlations between trait mean values and the other stability parameters were non-significant. Pi and CV cluster together in the majority of cases. ASV and AWAI consistently clustered with σ 2 and Wi, in most cases with DJi, s 2 di and to lesser extent with Si (1) and Si (2) . Mean value, bi and Bi_A formed a distinct group for GRYLD in Experiment One, whereas in Experiment Two, CV and Pi had the greater distances from the other parameters (Figure 3a,b). CV, bi and Bi_A formed a distinct group for TKW and VW (both experiments- Figure 4) and for CAR and PRO (Figure 5c,d).  parameter bi was negatively correlated with bi SPS (p < 0.05), while ASV NTLSM was negatively correlated with ASV SPS (p < 0.05). TKW parameters CV and bi were positively correlated with CV and bi VW in both experiments (p < 0.01 and p < 0.05 in Experiment One and Two, respectively). TKW parameters bi and CV were also positively correlated with the respective GDD parameter (p < 0.05) in Experiment One, while in Experiment Two, TKW ASV was positively correlated with ASV GDD (p < 0.05). Positive correlations were observed between bi, CV and ASV GDD parameters and their respective CAR parameters (p < 0.05).       GRYLD was positively correlated to NSSM in both experiments (Figures 6 and 7) and the correlations between these traits were significant to almost all environments and stronger in lower yield environments (data not shown). Significant negative correlations between NSSM and TKW were observed in both experiments, contrary to the strong positive correlations between NSSM and SPS. This trend was consistent in almost all environments. GRYLD was not related with GDD, however, significant negative correlations were detected in few environments (data not shown). On the contrary, correlations between stability parameters for different traits were in general low and non-significant in most cases and for both experiments. Higher correlations were recorded between traits bi and CV stability parameters than the others. GRYLD stability parameter bi was significantly correlated with bi TKW (p < 0.05) and bi PH (p < 0.05) in Experiment One and with bi NSSM (p < 0.01) and bi TKW (p < 0.01) in Experiment Two. GRYLD stability parameter CV was significantly correlated with CV TKW (p < 0.01) and CV PH (p < 0.01) in Experiment Two. In some cases, significant correlations were observed between yield components stability parameters, nevertheless, the significant level was inconsistent among stability parameters and between the two experiments. For example, NSSM parameter CV was significantly correlated with CV NTLSM (p < 0.01) in Experiment One, however this trend was detected neither in Experiment Two nor in the other parameters. NTLSM parameter bi was negatively correlated with bi SPS (p < 0.05), while ASV NTLSM was negatively correlated with ASV SPS (p < 0.05). TKW parameters CV and bi were positively correlated with CV and bi VW in both experiments (p < 0.01 and p < 0.05 in Experiment One and Two, respectively). TKW parameters bi and CV were also positively correlated with the respective GDD parameter (p < 0.05) in Experiment One, while in Experiment Two, TKW ASV was positively correlated with ASV GDD (p < 0.05). Positive correlations were observed between bi, CV and ASV GDD parameters and their respective CAR parameters (p < 0.05).        Traits' stability parameters were not significantly correlated with the traits per se in most cases and for both experiments (Annex S3). In general, the correlations were higher for CV and bi parameters than the others. GRYLD was significantly correlated only with bi NSSM in Experiment One. In a few cases, yield components were significantly correlated with yield components' stability parameters, although this trend was inconsistent between parameters and the two experiments.
The repeatability of the trait means values and the stability parameters derived from the eight cultivars that were common to the two experiments are shown in Table 4. Mean values were highly and significantly correlated for GDD, PH, and GRYLD. Pi was the stability parameter with highly significant correlations for all traits. CV and bi GRYLD parameters were significantly correlated, while AMMI derived parameters (ASV and AWAI), S 2 di and σ 2 had higher correlations for PH and GDD.

Discussion
Genetic variation between cultivars was observed for all traits, although the environmental effect explained higher percentage of variability for almost all traits reflecting the substantial variation in climatic conditions between locations and years. GRYLD varied from 689 Kg/ha up to 7192 Kg/ha, representing the expected range in durum wheat yield under Mediterranean conditions [3]. In agreement with previous studies, SPS showed the largest genetic control, and G × E interactions among yield components [54], CAR was largely under genetic control [10], and GI was the trait equally affected by the genotypes and the environment [4].
G × E interactions were significant in all cases, and the percentage of variance explained by the interactions was higher than the genetic effect for GRYLD. The superiority of the AMMI model to explain higher percentage of G × E interactions compared to linear regression models [8,21] was also evident in the present study. The usefulness of the linear regression models are doubtful if heterogeneity of the slopes do not reach significance [15] and if it explains a small part of the G × E interactions [55]. The heterogeneity of the slopes was not significant for GRYLD and yield components except TKW in Experiment One and the percentage of variance explained was generally low in both experiments. The linear regression model, however, was more efficient to explain G × E variation for phenotypic and quality traits, suggesting that it can be used with less caution for these traits.
Regression coefficient according to Finlay and Wilkinson [18] and Perkins and Jinks [20] are similar notions except that, in the latter case, observed values are adjusted for location effects before regression [32]. The MS of the deviation from regression in both cases is the measure of stability in the above two linear regression models. Wrickle's and Shukla's stability parameters are equivalent because their second terms are constant and their first terms are proportional to each other [32]. Therefore, the strong and consistent relations between bi and Bi_A, between s 2 di and DJi, and between σ 2 and Wi observed in this study are frequently reported [14,56].
Becker and Leon [15] stated that the information given by Wi is rather similar to s 2 di, confirmed by the high correlations reported between these parameters [14,57,58]. Lin et al. [32] recommended that when data do not fit to the linear regression models or if the residual MS from the regression are heterogenous, Wi and σ 2 parameters should be used. In the present study, the relations between Wi-σ 2 with s 2 di-DJi were stronger in cases where the percentage of variance explained by the heterogeneity of the slope of the linear regression models was low. This trend was also observed in the relations between s 2 di-DJi and AMMI-derived parameters (ASV and AWAI) while the strong relations between AMMI-derived parameters and Wi-σ 2 were consistent. Consistent correlations between σ 2 , s 2 di, and AMMI-derived parameters were previously reported [43]. The relation between CV, bi, and Pi and the correlation between these parameters with grain yield is in agreement with the results of Flores et al. [57], Mohammadi and Amri [58], and Vaezi et al. [14]. Correlations between nonparametric parameters Si (1) and Si (2) and Wi, s 2 di, and ASV were reported for durum wheat grain yield [56,58].
Based on the above, the stability parameters in the present study can be categorized into three groups. Group A includes AMMI-derived parameters (ASV and AWAI) which provide very similar information to Wi and σ 2 . Deviation from regression parameters (s 2 di-DJi) also fall in this group when the heterogeneity of the slopes was not significant. Nonparametric parameters Si (1) and Si (2) were moderately related with the former parameters. Parameters in group A were not related with the trait mean values indicating that the stability calculated based on these parameters is an independent trait and accord to the static (biological) concept of stability. Group B includes parameters bi, Bi_A, CV, and Pi. These parameters were significantly related with the trait mean values suggesting that breeders can employ them for simultaneous selection for stability and the traits per se. The attributes of parameters in group B were related to the dynamic (agronomic) concept. It should be stressed though that, contrary to bi, Bi_A, and CV, Pi had very high and consistent correlations with the trait mean values (over 0.95). Pi was associated with durum wheat yield improvement under Mediterranean conditions [54]. However, the very high correlation of Pi with the traits' mean values implies that this parameter may provide to large extent similar information with the trait per se. Group C encompass s 2 di-DJi when the heterogeneity of the slope was significant. The deviation from regression parameters have received some criticism concerning their employment as stability statistics. Lin et al. [32] argued about the efficiency of s 2 di to estimate stability when s 2 di is large or MS deviation is heterogenous. Flores and Cubero [57], on the other hand, stated that when regression coefficients are not significantly different, s 2 di becomes an important statistic in estimating stability. Our results signify that the adequacy of the statistical model employed needs to be considered before interpreting the genotypes' stability parameters derived by the model, information that is frequently missing from the literature. Such information was also reported to be crucial for assessing the repeatability of stability parameters over seasons [39].
The strong relations between yield, yield components, and phenology, particularly heading date and the existence of trade-offs between these traits, are well established [59] and they were evident in the present study. Sadras et al. [34] and Grogan et al. [37] showed associations between stabilities (plasticities) of yield and phenology. Associations between stability (plasticity) for yield, yield components, and physiological traits were also presented [35,36]. Relations between stability for quality traits were reported by Knapp et al. [38] and Peterson et al. [33]. In the present study, correlations between traits' stability parameters and the traits per se were, in general, low and non-significant in most cases and for both experiments indicating that contrary to the traits, stability parameters are likely to be independent and can be employed with less caution regarding possible trade-offs. This was particularly evident for the stability parameters of groups A and C (ASV, AWAI, Wi, σ 2 , s 2 di-DJi, Si (1) , Si (2) ). Stability parameters of group B (CV, Bi_A, bi, and Pi) showed a higher tendency to corelate and the same trend was observed in the correlation of these parameters with the traits per se. Therefore, selection for stability based on parameters of group B could result to the simultaneous selection for stability for other trait(s) or trade-offs, thus breeders need to consider possible implications when using them.
The correlations between traits' mean values derived from the eight cultivars that were common to the two experiments were high and significant for GDD, PH, and GRYLD. The respective correlations between stability parameters were less consistent and non-significant in most cases which signifies the well-established notion that genotypes' stability depends on the set of genotypes and environments included in the analysis [21,32,38,39]. These results are in agreement with previous studies reporting higher and consistent repeatability of grain yield than the repeatability of stability parameters [39,42]. Exceptionally, Pi was the only stability parameter showing high repeatability for almost all traits, although this parameter was highly and consistently correlated with the respective trait mean values, raising consideration about the independence of the information that it provides.
Regression coefficient (bi) and CV were the parameters that showed repeatability for GRYLD. Coefficient of variation (CV) had the highest heritability for grain yield among other parameters in wheat [60]. Jalaluddin and Harrison [42] also reported high repeatability of these two parameters for wheat while repeatability of bi was less consistent in the work by Annicchiarico [39]. Repeatability of regression coefficients for soybean yield were moderate and higher than σ 2 , s 2 di, Si (1), and AMMI-derived parameters except SIPC f [43].
It has been suggested that the linear model can be used for the prediction of untested environments as long as it is in the range of the characteristics of the environments used for testing [16]. Despite their limitations, Flores et al. [57] stated that linear regression models remain useful when the number of genotypes and environments considered in the analyses are sufficiently large and when there are no extreme environments that bias the regression slope. The significant correlation of the regressions coefficient (bi) for GRYLD denotes that the inclusion of the extreme location of Athalassa in Experiment One did not affect the repeatability of this parameter. This can be attributed to the extremely low yield obtained in location Dromolaxia during the growing seasons 2015/16 and 2016/17 in Experiment Two.
The significant repeatability of bi for GRYLD in the present study implies that regression models can provide repeatable regression coefficients when specific assumptions are met. However, repeatability of the correlation coefficient was lower for the other agronomic traits, particularly for the phenotypic traits PH and GDD, for which stability parameters of groups A and C (ASV, AWAI, Wi, σ 2 , s 2 di-DJi, Si (1) , Si (2) ) showed higher repeatability. Further to a well-established notion that several stability parameters should be employed in a plant breeding program because unrelated parameters provide independent information associated with different contexts of stability [14,56], our results imply that unrelated stability parameters can show different repeatability depending on the trait. The later needs to be further investigated in future studies.
Climate change is severely affecting agricultural production in the Mediterranean, including Cyprus. According to recent models, higher temperatures and lower precipitation are expected and this trend is likely to be more evident during the hottest periods of the year [12]. These climatic scenarios were represented in the current study by the location Athalassa where the adverse climatic conditions, particularly during grain filling, resulted in lower grain yield and deteriorated quality in terms of kernel and volume weight. These findings argue for new breeding strategies where emphasis should be given to the breeding of cultivars that can sustain high kernel and volume weight under adverse climatic conditions. The cultivar 'Simeto' managed to sustain high kernel weight in Athalassa despite its lower yield potential due to the moderate tillering capacity and lower number of SPS. On the other hand, cultivars 'Iride' and 'Adnan' managed to retain high NSSM in Athalassa due to their high number of SPS, although they also retained very low kernel and volume weight which is a point of concern. These findings denote that, cultivars which are characterized by high and stable kernel weight might represent an alternative yielding strategy in Mediterranean environments, in the effort to sustain both yield and quality. Accordingly, Royo et al. [61] found that cultivars with high kernel weight are more stable in Mediterranean environments.

Conclusions
Twelve stability parameters were evaluated for agronomic and quality traits in two experiments of twenty and sixteen cultivars, respectively, of diverse origin. Genetic variation between cultivars was significant for all traits. However, the environmental effect explained higher percentage of variability for almost all traits, reflecting the substantial variation in climatic conditions between locations and years. G × E interactions were significant in all traits. The AMMI model was more effective to explain higher percentage of G × E interactions compared to linear regression models. The linear regression model was more efficient to explain G × E variation for phenotypic and quality traits than for grain yield and yield components. Stability parameters were categorized into three groups. Group A included AMMI-derived parameters (ASV and AWAI), Wi, σ 2 , Si (1), and Si (2) . Deviation from regression parameters (s 2 di-DJi) also fell in group A when the heterogeneity of the slopes was not significant. Group B included parameters bi, Bi_A, CV, and Pi. Group C encompassed s 2 di-DJi when the heterogeneity of the slope was significant. Stability parameters in groups A and C were not related with the trait mean values contrary to the stability parameters in group B which showed significant correlations. In general, correlations between stability parameters for different traits and between traits stability parameters and the traits per se were low and non-significant, although a tendency for higher correlation was observed for stability parameters in group B. Stability parameters in group B had higher repeatability for grain yield while repeatability for growing degree days to heading and plant height was higher for stability parameters in groups A and C. The results of the present study allowed the detection of cultivars combining high performance and stability for agronomic and quality traits. The accumulation of such knowledge is crucial for adjusting plant breeding strategies and variety recommendations to farmers.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/agronomy12081774/s1, Annex S1: Normal climatic conditions based on historical records and climatic conditions during the experimentation years. Annex S2: Biplots of AMMI analysis. Annex S3: Spearman correlations between stability parameters and the traits perse.