Adaptability and Stability of Faba Bean (Vicia faba L.) Accessions under Diverse Environments and Herbicide Treatments

The adaptability and stability of 37 faba bean (Vicia faba L.) accessions with different levels of tolerance to metribuzin or imazethapyr was assessed across 12 season–location–herbicide experiments. Significant Genotype x environment (GE) interaction was found for the days to flowering (DFLR), plant height (PLHT) and grain yield (GY). Performance and stability of the accessions regarding PLHT and GY were assessed using four different stability parameters: cultivar superiority, static stability, Wricke’s eco-valence and Finlay and Wilkinson’s regression model. The stability parameters ranked these genotypes differently suggesting that PLHT and GY stability should be assessed not only on a single or a few stability parameters but on a combination of them. GGE biplot analysis indicated that the environments representing metribuzin treatment at Marchouch 2014–2015 and the non-treated treatment at Terbol 2018–2019 are the ideal environments for evaluating faba bean genotypes. GGE biplots showed herbicide tolerant accession IG12983 with simultaneous average PLHT, GY and stability across the environments. The performance of other tolerant accessions, namely IG13945, IG13906, IG106453, FB2648, and FB1216 was less stable but superior under specific mega environments. Therefore, utilizing these accessions in faba bean breeding programs would help broaden the adaptability to diverse locations–season–herbicide treatments.


Introduction
Faba bean (Vicia faba L.) was domesticated 10,000 years Before Christ (BC) in the Near East where archeological findings of domestic [1] and wild specimens were discovered [2]. Today, faba bean is considered the fourth most important cool season food legume after chickpea, lentils, and peas as it is grown on 2.57 million ha area with a total production of 5.4 million tons in 2019 [3]. This crop plays critical role in supporting nutritional and food security and enhancing soil structure in many countries, including China, Egypt, Ethiopia, United Kingdom, Australia, France, Sudan, and Morocco [3]. Faba bean production has increased 2% annually over the past three decades while global area remains stagnant.
Faba bean is affected by several biotic and abiotic stresses including parasitic and non-parasitic weeds. Among annual weeds, broadleaved and grass species like Anthemis arvensis L., Chenopodium album L., Convolvulus arvensis L., Sinapis arvensis L. and Avena sterilis L. compete with faba bean crop [4,5]. Parasitic weeds like Orobanche crenata L. and Cuscuta sativa L. also severely affect faba bean in many production areas [6][7][8]. An integrated weed management that employs a variety of chemical and non-chemical methods is an effective way of minimizing the losses caused by weeds. The development of resistant cultivars to multiple herbicides with different modes of action [9][10][11] is still the most economical and environmentally friendly control strategy to reduce the cost of the weed control practices, and avoid herbicides injury to the crop and herbicide resistance of some weeds [9,12,13]. Research conducted at ICARDA [9,14] resulted in the identification of accessions that tolerate metribuzin and/or imazethapyr herbicides. However, demonstrating the adaptability of these accessions to a wide range of environments can increase their economic value as climate change is expected to reduce the production of faba bean in many regions. Therefore, there is a need to study the yield stability of these accessions under different environmental conditions. Grain yield is a very complex trait which is strongly influenced by genotype (G), environment (E) and genotype x environment (GE) interaction [15,16]. GE interaction is of major importance for breeders, given that it reduces the association between phenotypic and genotypic values across environments [15,17]. It also affects the identification of relevant test environments, the allocation of resources within a breeding program and the choice of germplasm and breeding strategy [18]. GE interaction is a challenge in the case of legume breeding as previous studies have suggested a high proportion of variance due to environment (E) and GE interaction on the expression of grain yield in pulse crops including faba bean [19,20].
Environmental variation has a major effect on the variation of yield (up to 80% or higher) [21] in developed pure lines with narrow genetic base, but genotypes and GE interaction are more relevant for germplasm evaluation and selection and they must be considered simultaneously when selecting a genotype; in other words, an ideal genotype should have both high mean yield performance and high stability across environments [22,23].
The major objective of this research was to assess yield stability of selected faba bean accessions for herbicide tolerance under different environments with combined effect of herbicide treatment, location, and season. The second objective was to identify megaenvironments and the best environments to screen faba bean for herbicide tolerance.

Phenological Traits
Combined analysis of variance showed significant differences among the 37 accessions across the 12 environments for days to flowering (DFLR) and days to maturity (DMAT) reflecting the presence of genotypic variability for both traits. In addition, the interaction between Genotype x Herbicide treatment was highly significant (Table 1) indicating that the studied genotypes behaved differently under different herbicide treatments for DFLR and DMAT in the different environments. Very high values of narrow sense heritability (0.97 and 0.99) for both DFLR and DMAT were estimated under different herbicide treatments and across different locations-seasons (Table 1). Df. degree of freedom, DFLR days to flowering, DMAT days to maturity, PLHT plant height, GY grain yield, * p < 0.05, *** p < 0.001, h 2 narrow sense heritability.
The analysis of variance conducted for each environment showed significant differences (p < 0.001) among accessions in all environments for DFLR except in Marchouch 2016/2017-imazethapyr (environment E), where the season was dry and warm. The flower- ing time was earlier in the environments C, F, I and L where no herbicide treatment was applied at both locations and in different cropping seasons. Means and ranges of DFLR are presented in Table 2. This varied among accessions from 39 days at Marchouch 2016/2017 without herbicide treatments (environment F) to 53 days after sowing (DAS) at Marchouch 2016/2017-imazethapyr treatment (environment E) and the widest range was observed in Terbol 2015/2016-imazethapyr treatment (environment H) where it varied from 93 to 131 DAS. The maturity time was delayed in trials treated by both metribuzin and imazethapyr herbicides over season-locations than those with no herbicide application ( Table 2). The narrowest range of maturity time was observed in Terbol-2015/2016 with no herbicide treatment (environment I) where it varied between 165 and 171 DAS, and the widest range was observed in Marchouch 2014/2015 with no herbicide treatment (environment C) where it varied from 131 to 147 DAS (Table 2).

Plant Height
Combined analysis of variance showed that plant height varied significantly among genotypes and herbicide treatments, but no significant GE and Genotype × Herbicide Treatment interactions were observed across environments (Table 1, Figure 1). However, significant differences among genotypes for plant height were detected in all environments except in Marchouch during 2016/17-imazethapyr treatment (environment E) where severe terminal drought occurred. Narrow sense heritability of plant height was relatively high (0.60) indicating replicability of the traits among accessions in different herbicide treatments and across different locations-seasons (Table 1). h 2 _PLHT varied between 0.01 in Terbol 2015/2016-no herbicide treatment (environment I) and Terbol 2018/2019-metribuzin treatment (environment J) and 0.95 in Marchouch 2014/2015-metribuzin treatment (environment A) (Table 3). However, significant GE interaction was also observed for plant height indicating that the genotypes responded differently in different seasons and locations (Table 1).     All tested accessions had lower plant height under metribuzin and imazethapyr than under no herbicide application ( Table 2). Average plant height varied from 18 cm in Terbol-2015/2016-metribuzin treatment (environment G) and 112 cm in Terbol-2018/2019-no herbicide treatment (environment L). Figure 1 showed that non-treated plants tended to have the highest height, followed by the plants treated with metribuzin and then by those treated with imazethapyr. The plant height of accession IG104039 classified previously as tolerant to both herbicides did not differ significantly under metribuzin or imazethapyr treatments across all environments, and the plant height of accession VF513 classified previously as tolerant to metribuzin did not differ significantly under metribuzin treatment across environments ( Figure 1).
As significant GE interactions were detected for plant height, the stability parameters were assessed to determine the specific response of the tested accessions using the four stability parameters, namely cultivar superiority, static stability, Wricke's eco-valence [24] and Finlay and Wilkinson [25] stability parameter.
The rankings of the accessions based on the plant height stability are presented in Table 4. Considering the cultivar superiority and the ability of the genotypes to have a mean plant height above average across environments, IG13945 had the best plant height performance across environments. As for the ability of the genotypes to maintain a stable performance across environments, FB2601 was considered the most stable based on the static stability parameter, VF845 based on Wricke's eco-valence, as it had the lowest ecovalence value, and IG12659 based on Finaly and Wilkinson stability parameter, as it received the lowest values for these parameters.  Among the genotypes that had plant height above average and performed well in all environments, only four FB2601, IG104374, IG104421 and Flip86-98FB had small fluctuation across environments and were identified as stable by two of the four stability parameters. Among these genotypes classified earlier as moderately tolerant/tolerant to both herbicides, IG104374, IG104421 and Flip86-98FB were identified as stable by Wricke's eco-valence and FB2601 and Flip 86-98FB by the static stability parameter. Wricke's eco-valence followed by the static stability parameter was effective in simultaneously selecting stable genotypes with high plant height unlike the Finlay and Wilkinson's stability parameters which, in this study, identified mostly genotypes with low plant height as being the most stable.
The correlations among the stability parameters are shown in Table 5 The correlation coefficient varied between −0.4 and 0.9 indicating an inconsistency in the classification between the parameters. Significant negative correlation (−0.4) was observed between Cultivar superiority and Finley and Wilkinson's parameter and highly significant and positive correlation was observed between Finley and Wilkinson's parameter and Static stability (0.9) and between Wricke's eco-valence and Static stability (0.7). However, even with strong correlation between methods, genotype ranking can be different. FB2601, Flip86-98 and IG70622 had the most stable plant height as they ranked among the 10 most stable accessions by different stability parameters. Table 5. Correlation coefficients between the stability parameters of the 37 faba bean accessions tested across 12 environments.

Trait
Method Cultivar Superiority

Grain Yield
Combined analysis of variance revealed significant variation among genotypes and herbicide treatments (Table 1). Significant GE and Genotype × Herbicide treatment interactions were also observed for grain yield indicating that the genotypes responded differently to different environments characterized by different herbicide treatments, seasons, and years ( Table 1) (Table 3). This indicates diverse response of each accession to the different herbicide treatments and across the different locations-seasons combinations.
In addition, Genotype × Herbicide Treatment and GE interactions were highly significant indicating that the accessions performed differently under different herbicide treatments and in different locations and seasons (Table 1). Furthermore, Figure 2 indicated that some accessions yielded more under metribuzin and imazethapyr than under no herbicide application. This is also shown by the drown trends of each herbicide applications. The average grain yield was lower in environments treated with herbicides than in environments with no herbicide treatment at both locations in all seasons (Table 2). Figure 2 showed that the non-treated plants tend to have the highest grain yield, followed by the plants treated with metribuzin and then by those treated with imazethapyr. However, the tolerant and moderately tolerant accessions (FB2601 and IG13530) to both herbicides yielded more under metribuzin than under no herbicide treatment across environments, and the tolerant accessions (IG70622 and FB2528) yielded similarly under imazethapyr treatment and under no herbicide treatment across environments.
As significant GE interactions were detected for grain yield, the grain yield stability was assessed based on the same parameters used to evaluate the plant height stability as presented in Table 6. Considering the cultivar superiority and the ability of a genotype to have an above average mean performance across environments, VF845 had the best grain yield performance in all environments. As for the ability of the genotypes to maintain a stable performance across environments FB2528 was considered the most stable based on the static stability parameter, VF810 based on Wricke's eco-valence, as it had the lowest eco-valence value, and IG103043 based on Finlay and Wilkinson stability parameter, as they received the lowest values of these parameters. plants treated with metribuzin and then by those treated with imazethapyr. However, the tolerant and moderately tolerant accessions (FB2601 and IG13530) to both herbicides yielded more under metribuzin than under no herbicide treatment across environments, and the tolerant accessions (IG70622 and FB2528) yielded similarly under imazethapyr treatment and under no herbicide treatment across environments. As significant GE interactions were detected for grain yield, the grain yield stability was assessed based on the same parameters used to evaluate the plant height stability as presented in Table 6. Considering the cultivar superiority and the ability of a genotype to have an above average mean performance across environments, VF845 had the best grain yield performance in all environments. As for the ability of the genotypes to maintain a stable performance across environments FB2528 was considered the most stable based on the static stability parameter, VF810 based on Wricke's eco-valence, as it had the lowest eco-valence value, and IG103043 based on Finlay and Wilkinson stability parameter, as they received the lowest values of these parameters.    Among the genotypes that had grain yield above average in all environments, IG12983, which was tolerant to metribuzin and imazethapyr, was the only one that was able to maintain a stable performance across environments; it was identified as stable by Wricke's eco-valence. Wricke's eco-valence was the only effective parameter in identifying a stable and high yielding genotype, unlike the static stability and Finlay and Wilkinson parameters that identified mostly genotypes with low grain yield as being the most stable.
The correlation between the stability parameters is shown in Table 5. The correlation coefficient varied between −0.6 and 0.6. This result indicates inconsistency in the classification between the parameters. Highly significant and negative correlation (−0.6) was observed between static stability and Finley and Wilkinson's parameter and highly significant and positive correlation was observed between static stability and Wricke's ecovalence (0.6). However, even with strong correlation between methods, genotype ranking can be different between the methods. IG11527 and VF810 had the most stable grain yield as they ranked among the 10 most stable accessions by different stability parameters at the same time.

GGE Analysis
The GGE-biplots presented in Figure 3 accounted for 68.23% and 76.37% of the total variability for the plant height and grain yield, respectively. The environments with low narrow sense heritability values were excluded from the GGE biplot analysis as most of the variations observed are not genetic and might be related to the environmental conditions. The values of narrow sense heritability obtained in the present study shows

Evaluation of Test Environments
GGE-biplot provides a summary of the relationship between test environments. Two environments are positively correlated if the angle between their vectors is less than 90° [26]. Based on this, the plant height biplot (Figure 3a The plant height GGE biplot (Figure 3a) was divided into 6 sections where the 12 environments fell into two of these, and accordingly two mega-environments were

Identification of Mega-Environments and Specific Adapted Accessions
The plant height GGE biplot (Figure 3a) was divided into 6 sections where the 12 environments fell into two of these, and accordingly two mega-environments were identified. The first mega-environment contained three environments of  Figure 4a,b show the ranking of genotypes based on plant height and grain yield performance and stability in 12 environments. The mean yield performance and stability of genotypes were evaluated using an average tester axis (ATA) that passes through the origin [27][28][29]. Based on Figure 4a, 17 accessions were shorter than the average plant height across the 7 environments as they are located on the left side of ATA. The other 20 accessions were taller than the average across the 7 environments as they are located on the right side of ATA; 7 of them were tolerant to both herbicides' treatments across environments. Figure 4a shows that among the accessions tolerant to both herbicides with plant height above the average, FB2601 (16), IG104421 (18) and IG12983 (35) were the most stable as they had the shortest projection to the ATA. On the other hand, despite being moderately tolerant/tolerant to metribuzin and imazethapyr with good plant height performance across seven environments, IG13945 (28) was the least stable among all the accessions given that it had the longest projection to the ATA.

Performance of Tested Accessions
Regarding grain yield performance, Figure 4b shows that 18 accessions yielded less than the average grain yield across the seven environments as they are located on the left side of the ATA. The remaining 19 accessions yielded more than the average grain yield as they were located on the right side of ATA; 6 of them were tolerant to both herbicides' treatments across environments. Among the accessions that yielded more than the average grain yield and are tolerant to both herbicides, FB2648 (14), IG12983 (35) and VF522 (9) were the most stable as they had the shortest projection to the ATA. Hence, based on Figure 4a,b, the metribuzin and imazethapyr tolerant accession IG12983 (35) was considered a superior genotype in terms of plant height and grain yield as it showed a good and stable performance for both traits across environments. Furthermore, the following accessions VF963 (27), FB2648 (14), and IG70622 (33) that are tolerant to both herbicides performed well and were moderately stable across the 7 environments. of genotypes were evaluated using an average tester axis (ATA) that passes through the origin [27][28][29]. Based on Figure 4a, 17 accessions were shorter than the average plant height across the 7 environments as they are located on the left side of ATA. The other 20 accessions were taller than the average across the 7 environments as they are located on the right side of ATA; 7 of them were tolerant to both herbicides' treatments across environments.  Figure 4a shows that among the accessions tolerant to both herbicides with plant height above the average, FB2601 (16), IG104421 (18) and IG12983 (35) were the most stable as they had the shortest projection to the ATA. On the other hand, despite being moderately tolerant/tolerant to metribuzin and imazethapyr with good plant height performance across seven environments, IG13945 (28) was the least stable among all the accessions given that it had the longest projection to the ATA.

Discussion
Multi-environment trials are an integral part of the breeding pipeline to better understand the performance of tested accessions under a wide range of environmental conditions and therefore breeders are able to characterize the mega environments and identify cultivars adapted to specific environments or with broad adaptability [20,30]. In this study, faba bean accessions were evaluated in a wide range of environmental conditions created by different site-season-herbicide treatment combinations.
Faba bean accessions flowered and matured earlier in the environments with no herbicide treatment than in the environments with metribuzin or imazethapyr treatments in all sites and seasons. Past studies also reported a delay in flowering and maturity time of different legume crops treated with metribuzin or imazethapyr [31][32][33][34][35]. This delay might be the result of the growth inhibition of the crops amid their treatment with herbicide [36][37][38]. Faba bean accessions flowered and matured earlier in Marchouch than in Terbol. This might be attributed to cooler and more rainy weather in Terbol as compared to Marchouch. Past studies also reported decline in crop duration under heat and drought conditions in faba bean [39,40], lentil [41,42], chickpea [43], and common bean [44]. The earlier onset of flowering and maturity was observed in the non-treated environment of Marchouch 2016-2017 season; this was expected as it resulted from a combination of an exceptional warm and dry season and no herbicide treatment.
Plant height and grain yield were higher in the environments with no herbicide treatment than in the environments with metribuzin or imazethapyr treatments in all sites and seasons. Several studies also reported reduction in plant height and grain yield of accessions treated by herbicides in chickpea [34,45] and lentil [31,46]. This reduction might be due to the growth and photosynthesis inhibition caused by both metribuzin and imazethapyr [36,38,47]. Plant height and grain yield were higher in Terbol clustered under high rainfall than in Marchouch with low rainfall conditions. The highest average plant height and grain yield were recorded in the experiments with no herbicide treatments at Terbol 2018-2019; this was expected as it resulted from a combination of lowest temperatures, highest rainfall conditions and no herbicide application.
The heritability estimates are effective when combining data from diverse environments as the phenotypic value used to estimate the heritability is the mean value obtained across experiments and replicates [48]. In our study, heritability estimated for grain yield was highly affected by the stress environments followed by those estimated for plant height, and days to flowering and maturity. Similar observations were reported for heritability in faba bean by Toker [15]. Mohamed [49], Abdelmula et al. [50], Ceccarelli [51], and Atlin and Frey [52] concluded that lower heritability was expected in low-yielding environments. Therefore, the selection of faba bean genotypes is best done under optimum environments that are less likely to encounter stress-periods. Furthermore, moderate to high values of narrow sense heritability reported in the present study are important because the response to selection depends on the additive genetic variance captured by the narrow-sense heritability [53] and therefore they make a good basis for further genetic analysis and allow for true replication of a genotype in and across multiple environments [48,54].
Breeding programs focus on the evaluation of the performance and stability of accessions that have traits of economic importance under diverse environments. The stability analysis conducted in the present study allowed the identification of stable and high yielding genotypes across different environmental conditions. A stable genotype should have an above average and stable performance across environments [22,23]. The various stability parameters used in this study ranked plant height and grain yield of genotypes differently at different test environments. The inconsistency in ranking of cultivar superiority, Finlay and Wilkinson and Wricke's eco-valence indices were also reported in faba bean [21], pearl millet [55] and maize [56]. Our results agree with the results reported by Mustapha and Bakari [57] who observed no similarity between static and cultivar superiority but are contrary to the ones reported by Dehghani et al. [58] who observed similarity between Finlay and Wilkinson and cultivar superiority when ranking lentil genotypes.
Our results suggest that some genotypes had stable plant height and grain yield performance based on more than one parameter, but their rankings differed with each parameter. This implies that the comparisons may greatly depend on the parameter used as also observed by Milioli et al. [59] and Westcott [60] and thus more than one parameter should be used to characterize and explore performance of genotypes across environments and enable more reliable selection and recommendation of genotypes [61].
Our results also suggest that the selection for genotypic performance stability based on the static stability, Wricke's eco-valence and Finlay-Wilkinson parameters favor genotypes having plant height and grain yield lower than the population averages. Similarly, Temesgen et al. [21] and Fikere et al. [62] also reported that low-yielding faba bean and lentil genotypes were more stable than high-yielding ones.
Static stability was highly correlated with two other stability parameters for both plant height and grain yield. Seife and Tena [63] found that Wricke's eco-valence was positively correlated with all stability parameters. However, selecting genotypes based on this method exclusively may not be suitable to identify faba bean accessions that are high-yielding and stable. The use of the Finlay and Wilkinson parameter and Static stability as a selection tool would favor superior and stable genotypes for plant height and grain yield, respectively. Temesgen et al. [21] identified high yielded genotypes that show static stability despite the finding that both high yield and static stability rarely occur in multi-location trials. The classification of low yielding genotypes as stable and high yielding genotypes as unstable by the different stability parameters might be due to the type of accessions evaluated. In the present study, the evaluated accessions are pure lines that have narrow genetic base, narrow adaptability and generally are low yielding and unstable due to homozygosity [64,65].
Plant breeders routinely conduct GGE biplot analysis of multi-environment trials to identify ideal test locations, to reduce the cost of breeding and testing strategies, and to identify genotypes that are widely or specifically adapted. The partitioning of the total sum of squares through GGE biplots obtained in our study shows that there were differential plant height and grain yield performances of faba bean genotypes across environments and consequently a high GE interaction. This interaction could reduce the accuracy of genotype evaluation and selection process [66]. Many GGE studies have been carried out in faba bean and other crops but none of them covers the effect of herbicide treatments. The present study employed a GGE biplot to analyze data from multi-location trials carried out across different locations and under different herbicide treatments over three years. Herbicide application is greatly influenced by weather conditions [67][68][69] and therefore evaluation of the environments and genotypes with herbicide treatment is pertinent to identify genotypes with stable herbicide tolerance.
The GGE biplot was used to evaluate the test environments. An environment is considered ideal for genotype testing when it discriminates the genotypes and represents the environments [16]. The presence of correlation between two environments means that similar information about the genotype performance is derived from them [23] and therefore could be an option to reduce the number of test environments and, as a result, to establish a cost-effective breeding program. The correlations observed in our study between two environments are reliable as both plant height and grain yield biplots accounted for more than 60% of the total variation [29]. Yang et al. [70] claimed that a GGE biplot is considered useful if the two PCs account for more than 60% of the (G+GE) variability. As the GGE biplot that included all the environments accounted for only a small percentage (less than 60%) of the total variability, the patterns obtained were considered less useful and a more reliable and informative GGE biplot was obtained after excluding the low heritable environments.
According to Yan and Tinker [26], the test environments that are less discriminating provide little information on the genotype differences and should not be used as test environments. Hence, in this study among the seven test environments, environment G (metribuzin treatment of Terbol 2015/2016) is the ideal environment for plant height evaluation of genotypes as it is highly correlated with other environments (K, I, H, and B) and is highly discriminating. On the other hand, the ideal environment for evaluating grain yield of faba bean genotype is also environment G (metribuzin treatment of Terbol 2015/2016) as it is highly correlated with many other environments (H, A, B and C) and is highly discriminating. Our results suggest that the discriminating ability of environments was highly influenced by the weather conditions as the three environments of Marchouch 2016-2017 where the weather was exceptionally warm and dry were the least discriminating for grain yield evaluation and therefore, when choosing the testing sites for herbicide tolerance it is better to choose a site that is less likely to have a warm and dry spell period during the growing season like Terbol station.
GGE biplot is an effective visual tool for identifying the mega-environment issues and showing the specific adaptation of the genotypes and which cultivar won in which environments [29,71]. A mega-environment is defined as a group of locations that consistently share the same best cultivar(s) [27]. "Which-won-where" plots constructed in the present study grouped the test environments that represent a combination of seasonlocation-herbicide treatment into different mega-environments. However, the grouping did not correspond with the geographic location or herbicide treatment applied; the grouping seemed to be influenced by the weather conditions and the non-repeatability of the winning genotype in the same geographic location or under same herbicide treatments might be the result of the weather fluctuations observed in the same location from one season to another and to the effect that the weather conditions have on the efficiency of herbicide treatments.
To delineate a mega-environment, the consistency of the genotype's performance and the repeatability of the winning genotype in the same locations are necessary [72]. The reason for not meeting this condition in our study might be because plant height and grain yield biplots couldn't capture all the GE variation. Mega-environments are homogeneous groups of locations that reduce research costs by enabling fewer representative environments to be selected for genotype evaluation [73]. However, the identification of mega-environments is not easy. Many studies identified mega-environments in faba bean for grain yield and chocolate spot (Botrytis fabae) disease resistance [74], autumn or spring sowing adaptation [20,30] and resistance to Ascochyta blight (Ascochyta fabae) [75] and many other studies attempted to define mega-environments in spring wheat [76], sugarcane [77] and rice [78]. Our study is the first attempt to identify mega environments that englobe diverse herbicide treatments.
An ideal genotype is defined as one of the highest yielding across the test environments and is stable in performance [16]. The metribuzin and imazethapyr tolerant accessions which showed a good and stable plant height and grain yield performances in the current study are promising and of great importance for faba bean growers. As plant height is highly correlated with the biological yield [39,79], the accessions showing high and stable plant height performance are very important for the regions where faba bean is grown mainly for animal feeding such as Ethiopia and the Mediterranean countries. The environments evaluated in the present study represented different herbicide treatments applied in different sites and seasons. Hence, the genotypes mentioned are less sensitive to the environmental changes, have high yield, and are tolerant to metribuzin and imazethapyr treatments. Some of the high yielding faba bean genotypes reported in the present study were sensitive to environmental changes. Similar results were also obtained in other crops such as soybean [80,81].

Materials
Thirty-seven faba bean accessions with different level of tolerance to one or both herbicides, metribuzin and imazethapyr, were evaluated in this study. Thirty-six accessions were pure lines derived from single plant selections in three consecutive seasons which represent landraces from 21 countries and one cultivar Elisar (Flip 86-98FB) released in Lebanon was used as check (Table S1). Among them, 14 were tolerant or moderately tolerant to both the herbicides @ 250 g ai/ha metribuzin and 75 g ai/ha imazethapyr, six were tolerant to 250 g ai/ha metribuzin but sensitive to imazethapyr, and 16 were tolerant to 75 g ai/ha imazethapyr but sensitive to metribuzin ( [9], ICARDA, unpublished data).  Table 7. Each combination of herbicide treatment, location and season was considered as an independent environment. The experiments were conducted in an incomplete block design with two replicates. The plot size was 2 rows with one-meter length with rows spaced at 45 cm and plants spaced at 20 cm within the row. Terbol station is characterized by a cool high rainfall winter and a moderate wet spring. The rainfall recorded during 2015/2016 and 2017/2018 winter cropping seasons were 343 mm and 810.2, respectively. A supplemental irrigation of 30 mm was provided during 2015/2016 to compensate the dry spell periods. Spring 2017/2018 was warmer than normal spring seasons ( Table 7). The soil in Terbol station is deep and rich clay loam. Marchouch station is characterized by semi-arid environment, low rainfall and moderate temperature during winter and spring seasons. The annual rainfall recorded during 2014/2015 and 2016/2017 cropping seasons were lower than the mean annual rainfall (396 mm) and the spring season of 2016/2017 was relatively warmer than normal spring seasons. The soil at Marchouch is Vertisol mostly silty clay ( Table 7). The experiments were supplied with 250 kg/ha of granulated NPK (15:15:15) during land preparation. The experiments were sown in late November at Terbol and mid-December at Marchouch and harvested by the end of May at both locations. Necessary phytosanitary and agronomic management practices were applied to ensure a good crop stand: lambda-cyhalothrin @ 40 g ai/ha was sprayed to control sitona weevil (Sitona lineatus L.), imidacloprid @ 160 g ai/ha was sprayed to control aphids and azoxystrobin @ 72.8 g ai/ha and difenoconazole @ 45.6 g ai/ha were spayed alternatively to control foliar diseases. To control weeds in all environments, we applied 1200 g ai/ha of pendimethalin as preemergence treatment in addition to post-emergence herbicides. For post-emergence herbicide treatments, two herbicides, namely Metribuzin @250 g ai/ha and Imazethapyr @75 g ai/ha were sprayed uniformly at the stage of inflorescence emergence-BBCH stage 50 [82,83] in addition to a control untreated treatment where hand weeding was applied to ensure an unbiased evaluation of the performance and stability of the faba bean accessions as they were evaluated under the same conditions except for the herbicide treatment.

Recorded Traits
The following traits were measured based on the faba bean ontology described by Maalouf [84]: days to flowering (DFLR) and days to maturity (DMAT) were recorded on plot basis while grain yield per plant (GY) and plant height (PLHT) were recorded on three random plants from each plot and averaged.

Statistical Analysis
The spatial statistical model was applied for variance analysis using the Automatic Spatial Analysis of incomplete block design modules of GenStat 19 edition [85]. Variation among accessions and treatments was assessed in terms of p-values using the Wald statistic, and the best unbiased phenotypic estimates of accessions (A) were estimated with standard error using best linear unbiased prediction values (BLUP) using GenStat software. Narrowsense heritability for each environment (h 2 ) was estimated for the PLHT and GY using the method of residual maximum likelihood (REML) and combined narrow sense heritability (h 2 ) was estimated for the traits based on the combined analysis using REML model and Best unbiased estimated values of Genstat 2019.
In our study the environment is defined as the combination of year-locations and herbicide treatments either metribuzin or imazethapyr. The following four stability parameters were assessed using GenStat software by comparing different treatments and environments: cultivar superiority index [61], which refers to the ability of the accession to perform above the mean in different environments, static stability coefficient; [86] which refers to the consistency of accession's performance across different environments; Wricke's ecovalence [24], which refers to the contribution of the accession to GE interaction; and the index of Finlay and Wilkinson [25], which refers to the response of an accession to different environments by fitting a regression of the environment means for each accession on the average environmental means.
The GGE biplot is an ideal tool for the analysis of data from multi-environment trials (MET); it considers both G and GE interaction effects and graphically displays GE interaction in a two-way table [87]. GGE biplot allows visual examination of the relationships among the test environments, the performance and stability of the genotypes and the mega-environment analysis to recommend specific genotypes for specific mega-environments [16,20,26].
GGE biplot analyses of tested accessions were conducted using the BLUPs obtained under diverse herbicide treatments at two locations: Terbol which is characterized by high rainfall, and Marchouch by low rainfall. The environments with low narrow sense heritability for both GY and PLHT were excluded from the GGE biplot analysis as most of the variations are related to the environmental conditions.
The relationship between the environments was visualized by drawing a vector that connected each environment to the biplot origin:

•
The correlation between two environments was approximated based on the angle between two vectors [26,88]; the smallest the angle between two vectors, the highest is the correlation between the two environments.

•
The discriminating ability of the test environments was evaluated based on the length of the vector of each environment; the longer the environment vector, the more the discriminating ability of the environment.
The best genotypes for each environment and the possibility of existence of mega environments were identified using the "who-win-where" visualization [29,87]; the polygon view of a biplot (convex hull) is the best way to visualize the interaction between the genotypes and the environments [16]; each polygon was formed by connecting the genotypes that are farthest from the biplot origin so that all other genotypes are inside the polygon [20]. The perpendicular lines to the sides of the polygon divide the biplot into sectors. Each sector has a vertex genotype. The vertex genotype is the one having the longest vector and is considered the winning genotype. The mega-environment was identified as the group of environments that share the same winning genotypes following Yan and Rajcan [27].
A genotype is considered superior when it has both high mean performance and high stability across the test environments. The mean yield performance and stability of genotypes were evaluated by an average environment coordination (AEC) method [27][28][29]:

•
The mean performance of the genotype was graphically evaluated based on the line perpendicular to the average tester axis (ATA) that passes through the origin and separates entries with below-average means from those with above-average means; the genotypes located on the right side of this line are taller or have more yield than the ones located on the left side.

•
The stability of the accessions was graphically represented by the projection from the genotype to the ATA; the longer the projection the greater is the GE interaction and therefore the lower the stability of the genotype across environments.

Conclusions
The present study shows significant Genotype × Environment interaction for grain yield and plant height which highlights the need to select genotypes well adapted to specific environment as well as broadly adapted genotypes. The performance and stability of faba bean genotypes were analyzed using four different stability parameters. These stability parameters showed inconsistency in the ranking of genotypes and showed that different stability parameters tend to favor low yielding genotypes. These parameters may not be appropriate, as both breeders and farmers prefer to adopt genotypes that are high-yielding and at the same time perform consistently across environments. GGE identified IG12983 as a superior genotype in terms of plant height and grain yield as it showed a good and stable performance for both traits across locations, seasons, and herbicide treatments which make it attractive to the farmers as it can provide an effective weed management tool. Some accessions had specific adaptation to at least one of the defined mega environments but not to others; IG13945 and IG13906 were specifically adapted to one mega environment in terms of plant height while IG106453, FB2648 and FB1216 were specifically adapted to one mega environment in terms of grain yield. Moreover, some accessions were high yielding but unstable across environments while others were low yielding and stable. To develop superior herbicide tolerant cultivars that are broadly adapted to different mega environments, there is a need to cross the tolerant germplasm identified in the present study with other cultivars adapted to specific environments. These lines could also be crossed with other cultivars to accumulate traits with economical interest in new faba bean varieties. Furthermore, this study suggests conducting herbicide screenings under environments that are less likely to experience drought to avoid the confounding effect of herbicides and drought.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/plants11030251/s1, Table S1: Faba bean accessions with different degree of tolerance to Metribuzin and Imazethapyr used in the present study.