Identiﬁcation of the Optimum Environments for the High Yield and Quality Traits of Lentil Genotypes Evaluated in Multi-Location Trials

: Lentil is a versatile and proﬁtable pulse crop with high nutritional food and feed values. The objectives of the study were to determine suitable locations for high yield and quality in terms of production and/or breeding, and to identify promising genotypes. For this reason, ﬁve lentil genotypes were evaluated in a multi-location network consisting of ten diverse sites for two consecutive growing seasons, for seed yield (SY), other agronomic traits, crude protein (CP), cooking time (CT) and crude protein yield (CPY). A signiﬁcant diversiﬁcation and specialization of the locations was identiﬁed with regards to SY, CP, CT and CPY. Different locations showed optimal values for each trait. Locations E4 and E3, followed by E10, were “ideal” for SY; locations E1, E3 and E7 were ideal for high CP; and the “ideal” locations for CT were E3 and E5, followed by E2. Therefore, the scope of the cultivation determined the optimum locations for lentil cultivation. The GGE-biplot analysis revealed different discriminating abilities and representativeness among the locations for the identiﬁcation of the most productive and stable genotypes. Location E3 (Orestiada, Region of Thrace) was recognized as being optimal for lentil breeding, as it was the “ideal” or close to “ideal” for the selection of superior genotypes for SY, CP, CT and CPY. Adaptable genotypes (cv. Dimitra, Samos) showed a high SY along with excellent values for CP, CT and CPY, and are suggested either for cultivation in many regions or to be exploited in breeding programs.


Introduction
Lentil (Lens culinaris Medik.; Fabaceae) is among the major cool season annual grain legume crops, and is best adapted to semi-arid regions of the Mediterranean Basin and the Indian subcontinent. It is among the oldest crops, together with cereals like barley and wheat [1], and is probably the oldest domesticated grain legume because carbonized remains of lentil were found in the Franchthi cave in Greece, dated 11,000 BC, and from Tell Mureybit in Syria, dated 8500-7500 BC [2,3]. Lentil is the fourth major cool season grain legume in the world after dry pea (Pisum sativum), chickpea (Cicer arietinum) and broad bean (Vicia faba) produced for human consumption, and it constitutes an excellent source of complex carbohydrates, protein, minerals, vitamins and dietary fibers for humans, as well as being highly valuable as feed and fodder for livestock [4][5][6][7][8]. As well as its prominent role in the food and feed industry, lentil also grown for the diversification and intensification of cropping systems [9]. It exhibits ecological advantages as a rotational crop in cereal-based cropping systems by acting as break crop for the better control of pests, weeds and diseases, and the better management of herbicide residue, and at the same time plays an important role in maintaining soil fertility due to its capacity to fix atmospheric nitrogen (N 2 ) [10]. For instance, the range in quantity of N 2 fixed by Lens esculents Medik. through bacteria (Rhizobium) reaches 40-68 kg ha −1 year −1 [11,12].
In 2019, lentil was grown on more than 4.8 M ha across 52 countries, accounting for an annual production of 5.7 Mt [13] whereas, between 2010 and 2015, it accounted for 6% of the total dry pulse production, with an average yield of 926 kg ha −1 [14]. However, despite its agronomic and nutritional value, until recently, lentil seed production remained at low levels and attracted much less attention from plant breeders than cereal grains [15,16].
For lentil genotypes to become commercially successful, high production levels and yield stability are required under variable farming conditions. However, lentil production is controlled by several biotic and abiotic factors such as drought, salinity, high temperatures and mineral deficiency, and is trialed in many regions of the world under constant pressure from biotic and abiotic stresses, in the context of climate change and the food security crisis [17]. In order to achieve this goal, multi-location trials are conducted annually with the purpose of identifying superior genotypes for the target locations. As lentil can be grown in a wide range of environments, the yield of the tested genotypes across different locations and over the years can be different because of its high genotype (G) × environment (E) interactions. In plant breeding programs, the structure of the G × E interactions is very important in order to achieve a high rate of adoption of new improved cultivars [18]. It is often difficult to determine the pattern of genotypic responses across environments. Moreover, lentil, which has a relatively a narrow genetic base, may exhibit a low buffering capacity against environmental changes, which challenges this response [19].
The G plus G × E (GGE) biplot model provides a powerful solution to this problem as it graphically displays the genotype (G) plus the interaction G × E (Environment) of multi-environment data in a way that facilitates visual evaluation and mega-environment identification; that is, a group of locations that consistently share the same best genotype or genotypes. In addition, the GGE biplot technique is useful in selecting superior genotypes and test environments for a given mega-environment [20,21].
The development of new lentil cultivars with desirable seed quality attributes is possible because primary traits such as total starch, protein, 1000-seed weight, seed color, and other quality characteristics are significantly controlled by genetic factors, as well as by local pedoclimatic conditions [22]. It is also widely accepted that G affects seed lentil chemical properties considerably [23][24][25]. The seed protein of lentil is an important nutritional trait that has particular meaning for the food and feed industry. Protein quantity and protein yield are complex quantitative traits that depend on a number of factors and are highly influenced by genetic and environmental fluctuations [26]. The average crude protein in lentil seeds is about 26% on a dry weight basis [27], although several studies have reported seed crude protein in cultivated lentil ranging from 10.5-36.4% dry weight depending on the number of tested lentil genotypes, environments and methods of analysis [28,29].
The seed cooking quality is another important factor for lentil's utilization as a food crop because the seeds are largely consumed in their cooked form. This key culinary attribute of lentils is associated with the ease and the low cost of food preparation [5]. The cooking time, which is defined as the time needed from the start of boiling until the seeds become soft, constitutes one of the most important quality characteristics, and is believed to be a heritable trait, although the environmental and storage conditions are also known to be influencing factors [30,31]. In addition, high air temperatures and soil moisture are important in increasing the cooking time of lentils [32]. Furthermore, different cooking times were found when the lentil was cultivated in different locations, and this was mainly attributed to different levels of soil fertility [33][34][35]. The duration of the cooking time is of central importance in breeding programs, because a longer time to cook of a new variety may be penalized in the marketplace owing to the lower demand for such a product. In contrast, a quick-cooking variety may achieve premiums because of cooking convenience and reduced energy cost requirements [36].
Lentil is an important pulse crop in Greece, traditionally consumed as a soup with whole seeds, and it shows an increasing pace of production; however, its production levels do not meet the domestic consumption requirements, which are largely covered by imports. The dietary value of lentil is currently gaining recognition, while a remarkable shift of consumers to this crop is attributed to its high nutritional and culinary value. Lentil, as with many other pulses, can be regarded as key component of lowering the environment footprint of food production and consumption, and as an important ingredient of a sustainable diet. The former is mainly attributed to lentil's ability to fix atmospheric N, both for its own use and for the use of subsequent crops, thus replacing the expensive synthetic N fertilizers and thereby reducing greenhouse gas emissions. The latter emanates of the fact that lentil seeds are broadly consumed as sources of dietary proteins by over one billion people and constitute an excellent source of plant-based proteins for human food and animal feed. Nowadays, high lentil seed yields, stability and seed quality traits are attainable by the choice of the most adapted cultivar in a given pedoclimatic environment. In Greece, the majority of the cultivated varieties are commercial cultivars; however, there is a small number of lentil landraces with interesting agronomic, physiological and quality traits. The edaphic and climatic conditions play a key role for the achievement of elevated production and enhanced seed quality, although their interaction with the genotype has not yet been thoroughly explored. The objectives of this research-work were to identify: (i) the optimum locations in terms of high production and quality by investigating G × E interactions through a multi-location experimental network, (ii) suitable locations for lentil breeding for yield and quality traits, and (iii) the most suitable genotypes among those registered in the Greek National Varietal Catalogue to be cultivated by lentil growers, and for possible exploitation by future breeding programs targeting high yield and/or seed quality traits.

Genetic Material
Five lentil genotypes (four commercial cultivars and one improved population) were evaluated for two consecutive growing seasons (2018/2019 and 2019/2020, hereafter mentioned as 2019 and 2020, respectively) at ten (10) locations from various pedo-climatic zones in Greece where lentil is widely cultivated.
The commercial cultivars were 'Samos' (G1), 'Dimitra' (G2), 'Thessalia' (G3) and 'Elpida' (G4), and the local population was the '03-24L' (G5). All of the cultivars were developed by the Institute of Industrial and Forage Crops (IIFC) in Larissa, Greece, while the population (G5) originated from ICARDA and was improved by IIFC. G5 was improved with the method described by Tokaltidis and Vlachostergios [37]. The genotypes were selected based on their popularity among growers, their high yield potential and their differentiation in flowering and maturity earliness. Their basic characteristics with respect to cotyledon color, seed size, and earliness can be found in Table 1.

Test Locations, Pedoclimatic and Growing Conditions
The 10 test locations (listed with codes from E1-E10; Figure 1) included 5 Regions across Greece, as follows: south, the region of Sterea (St); central, the region of Central Greece (Ce); north, the regions of Central Macedonia (CM) and West Macedonia (WM); and northeast, the Region of Thrace (Th).

Test Locations, Pedoclimatic and Growing Conditions
The 10 test locations (listed with codes from E1-E10; Figure 1) included 5 Regions across Greece, as follows: south, the region of Sterea (St); central, the region of Central Greece (Ce); north, the regions of Central Macedonia (CM) and West Macedonia (WM); and northeast, the Region of Thrace (Th). Their geographical co-ordinates are given in Table 2. All of the locations had altitudes lower than 150 m above sea level, except of locations E1, E5 and E10, with elevations of 500 m, 624 m and 476 m respectively.
The climate in Greece is typically Mediterranean, characterized as semiarid by mild and rainy winters, relatively hot and dry summers and prolonged sunshine periods for most of the year. However, due to its specific geographical position in the Mediterranean and its diverse topographical relief, Greece is characterized by various climatic zones. The study locations included, basically, three types of climates according to the Köppen-Geiger classification [38,39]. BSk (arid, steppe, cold) in E1, E8 and E9; Csa (temperate, dry summer, hot summer) in E2 and E3; a combination of BSk/Csa in E6 and E7; and Cfa (temperate, without dry season, hot summer) in E4, E5 and E10 (Table 2). Their geographical co-ordinates are given in Table 2. All of the locations had altitudes lower than 150 m above sea level, except of locations E1, E5 and E10, with elevations of 500 m, 624 m and 476 m respectively.
The climate in Greece is typically Mediterranean, characterized as semiarid by mild and rainy winters, relatively hot and dry summers and prolonged sunshine periods for most of the year. However, due to its specific geographical position in the Mediterranean and its diverse topographical relief, Greece is characterized by various climatic zones. The study locations included, basically, three types of climates according to the Köppen-Geiger classification [38,39]. BSk (arid, steppe, cold) in E1, E8 and E9; Csa (temperate, dry summer, hot summer) in E2 and E3; a combination of BSk/Csa in E6 and E7; and Cfa (temperate, without dry season, hot summer) in E4, E5 and E10 (Table 2). In the 2019 growing season (November to July), the precipitation ranged from 367.2 mm (E3) to 676.3 mm (E2), whereas the 2020 growing season was wetter, with precipitation ranging from 402.5 mm (E5) to 731.9 mm (E2). Similarly, the precipitation between April and May (the early reproductive till pod filling growth stages) was lower in 2019, ranging from 68.6 mm (E4) to 144.7 mm (E5) in comparison to wetter the 2020 growing season, which ranged from 99.0 mm (E6) to 186.4 mm (E1). Interestingly, among the 10 test locations, the northern locations had the lowest season (November to July) precipitation (i.e., at E3, 367.2 mm in 2019, and at E5, 402.5 mm in 2020, respectively), whereas the season precipitation was higher in Ypato-Thiva (St) (E2) in both years (676.3 mm in 2019 and 731.9 mm in 2020), which was the southernmost location ( Table 2). The season (November to July) mean temperature was relatively higher in 2020 compared to 2019. The lowest temperatures (10.1 • C and 11.1 • C) were recorded in E5, and the highest (14.1 • C and 15.1 • C) in E2 in the growing seasons 2019 and 2020, respectively.
Under the influence of soil genesis factors (parent material, time, climate, organisms, and topography), soils in the broader tested area have been classified in the orders of Fluvisols, Cambisols, Vertisols, Luvisols and Regosols [40,41]. Briefly, Fluvisols are mainly young azonal soils in alluvial deposits (floodplain), lacustrine (lake) and marine deposits. Cambisols are moderately developed soils with weak horizon differentiation, and they generally make good agricultural land and are used intensively. Vertisols are productive soils if properly managed, but become very hard in the dry season and sticky in the wet season. Luvisols, commonly found in colluvial deposits of limestone weathering, are fertile soils suitable for a wide range of agricultural uses. Finally, Regosols are young mineral soils with little profile development.
In all of the trials, the genotypes were arranged in plots following a randomized complete block design (RCBD) with four replications. The plots consisted of seven (7) plant rows of 4m length each, with a distance of 0.25 m between the rows. The seeding rate was 80-120 kg ha −1 and the sowing was continuous in the furrows, with the uniform distribution of the seed. The sowing was performed during the last week of November 2018 and November 2019 in all of the locations, with the exception of location E5, where the sowing date was at the last week of February due to the very low temperatures that prevail during the winter period.
The trials at the 10 locations were accounted as rain-fed environments. In the two growing seasons (2019 and 2020) in each location, the trials were initiated on the same position of the field or in a nearby position of the same field. A balanced basal fertilization with 160 kg ha −1 (0N-46P 2 O 5 -0K 2 O) was incorporated pre-sowing, and the phytosanitary practices such as weed removal were performed by hand; diseases and pest control were followed according to the local recommendations.

Yield and Agronomic Measurements
The agronomic traits that were recorded from the three inner plant rows in each plot incuded the seed yield (SY), plant height after flowering (PH), number of pods per plant Sustainability 2021, 13, 8247 7 of 23 (PP) and 1000-seed weight (1000 SW). The description of the growth stages (phenological development of 50% of the plants in each plot) was conducted using the extended BBCH scale [48]. The PH and PP were determined as the average of ten (10) randomly selected plants per plot. PH was recorded after the flowering was completed (BBCH growth stages 69-71), while the PP was determined at physiological maturity (BBCH stage 89). In order to assess the SY, the plots were harvested by hand and threshed using a laboratory thresher (Wintersteiger LD350) when the plants of each genotype reached their physiological maturity. The SY was calculated on a plot basis (3 m 2 per plot) and converted to t ha −1 after adjusting to 13% seed moisture content.

Protein and Cooking Time
The cooking time of each lentil sample was estimated according to the tactile method [49], with some modifications. Briefly, 5 g lentil seeds were transferred in beakers (100 mL capacity) covered with aluminum foil containing 50 mL of boiling distilled water (ratio 1:10). The beakers were then placed on a hotplate, and after approximately 20 min of cooking, the softness of the cooked seeds was checked every 30 seconds. The cooking was considered complete when the seeds were soft enough to make uniform transparent mass (with no opaque core in the seeds) when placed between two glass slides and pressed against each other with no lateral movement. The optimum cooking time in minutes was determined as an average of triplicate determinations.
For the assessment of the seed crude protein concentration (% CP) and seed crude protein yield (CPY) at physiological maturity, samples consisting of the 3 inner plant rows of the plots were cut, bagged and dried to a constant weight in a dryer room. These samples were threshed whilst the seed was weighed for yield and monitored for moisture content. The seed crude protein concentration (% CP) was determined on uniformly ground samples obtained from all plots for each genotype using the Kjeldahl method, and expressed as the % CP = total N × 6.25. The seed crude protein yield (CPY) was the product of the SY and % CP content. The data presented for the % CP and CPY are expressed on a 0% moisture basis.

Statistical Analysis
The sources of variation and the appropriate F-ratios were estimated through analysis of variance (ANOVA) for balanced combined randomized complete block experimental designs (RCBD) over the years (Y) and locations (L), considering genotypes (G) and L as fixed treatment effects (variables of interest), and Y as random effects [50]. In order to estimate the contribution of each source of variation to the variability of each trait studied, the ANOVA treatment sum of squares (SSTRMT = SS G + SS L + SS Y + SS G×L + SS Y×G + SS L×Y + SS G×L×Y ), including G, L, Y, and their respective two-and three-way interactions, was partitioned into its components (SS G , SS L , SS Y , SS L×Y , SS G×Y , SS G×L , SS Y×L×G ) as the percentage of the sum of squares of the SS TRMT .
The Shapiro-Wilk test for normality showed that the variables were normally distributed, whereas the ANOVA's assumptions for the equality of the error variances and residual normality were met by Levene's test and the normal quantile-quantile (QQ) plot method, respectively [51]. The trial means were compared using Fishers' least significant difference test (LSD) at p ≤ 0.05, whereas Pearson phenotypic correlation coefficients (r) between the agronomic traits, yield and quality were calculated in order to identify the most prominent correlations. The statistical analyses were conducted using the statistical software JMP 11.0.0 (SAS Institute, 2013) [52] and IBM SPSS package v. 23 (IBM Corp., 2015) [53].
A genotype (G) plus genotype × environment (GGE) biplot model was used: (i) To generate graphs showing the polygon view of the pattern (i.e., which genotype had the highest yield in which location, or which genotype had the highest value in which trait). The combinations of G × Y were equated to G × L in order to reveal cross-over interactions across the two years of experimentation. (ii) For the ranking of the genotypes against the "ideal" (the ideal is the virtual highest yielding and stable genotype defined by the center of concentric circles on the biplot graph).(iii) On the same biplots, test environments were ranked against the "ideal", which was the virtual environment with the highest discriminating ability and representativeness for each trait [21]. The GGE biplot analysis was performed by the software package GGE Biplots in R version 1.0-8 [54].

Seed yield, Yield Components and Plant Height
The genotype (G), location (L), and the two-way (Y × L, G × Y, G × L) and three-way interaction (Y × L × G) were highly significant (p < 0.01) for SY, PP, and 1000 SW. Less consistent were the results for PH: although the year (Y) effect was significant only for PH, this trait was further affected by the L and G, and the two way (Y × L and G × L) interactions (Table 3). Table 3. Analysis of variance (ANOVA), sum of squares (SS) and partitioning of the treatment SS (SS TRMT ) of the five lentil genotypes tested across the ten locations for two growing seasons (2019 and 2020). df, degree of freedom; SY, seed yield (t ha −1 ); PP, number of pods per plant; 1000 SW, thousand-seed weight; PH, plant height (cm); CP, crude protein concentration (%); CT, cooking time (min); CPY, crude protein yield (t ha −1 ); ns: non-significant (p > 0.05); * and ** are the significance at p < 0.05 and 0.01, respectively.

Source
For SY, PP and PH, L contributed to most of the variation-54.3%, 65.4% and 60.7%, respectively-whereas Y accounted for only 0.0-2.0% of SS TRMT . The G × L contribution to the variation was also much higher in comparison to G × Y. The SS G×L /SS G×Y ratio was 4.6-fold higher for SY, and for this reason the GGE biplot with the polygon view based on the G × L interaction was used to study the genotypic adaptation across the two years. Among the agronomic traits, the 1000 SW was highly genetically controlled, given that the G effect accounted for the 82.2% of SS TRMT . Over the Y and G, E8 and E3, followed by E7, were the most productive locations, with a mean SY of 2.32 t ha −1 , 2.27 t ha −1 and 1.96 t ha −1 , respectively.
Accordingly, of the yield components, location E8 had the highest PP (77.6 pods plant −1 ), whereas location E3 exhibited the highest 1000 SW (48.5 g). The tallest plants over Y and G were recorded in location E2, with PH = 39.5 cm ( Table 4). The growing season of 2019 was more productive in SY by 9.3% in comparison to 2020, whereas there was small differentiation for the rest of the agronomic traits. Over Y and L, G1 followed by G2 were the best performing genotypes in SY (1.73 t ha −1 and 1.65 t ha −1 , respectively, Table 4). G1 and G2 had the highest PP (59.0 and 57.5, respectively), whereas G4 had the  The two annual GGE biplots, based on the SY data of the 2019 ( Figure A1; Appendix A) and 2020 growing seasons ( Figure A2; Appendix A), were constructed with the polygon view. The polygon view of the biplot presents which genotypes performed the best in one or more locations, in order to reveal the presence or absence of crossover G × L interactions, which are suggestive of the absence or existence of different mega-environments. However, the data show that the location grouping pattern was not repeatable, but seemed to be year dependent. Thus, e.g., G3 was the winning genotype in the different environments in each growing season, and there was not repeatable grouping. The across-the-years GGE biplot indicated three regional groupings, as follows: G1 for WM, West Macedonia; CM, Central Macedonia; Th, Thrace; G2 for St, Sterea; and G5 for Ce, Central Greece; Figure 2.
Across the ten locations, the closest to the "ideal" for SY and stability was G1, followed by G2; however, G1 was the most stable ( Figure 3a). The "ideal" locations in terms of representativeness and discriminating ability were E3 and E4, followed by E10 (Figure 3b).

Protein and Cooking Time
The year (Y), G, L, and the two-way (Y × L, G × Y, G × L) and three-way interactions (Y × L × G) were highly significant (p ≤ 0.01) for CP, CT and CPY ( Table 3). The cooking time was genetically controlled, with G explaining 38.6% of the SS TRMT , whereas CP and CPY were mainly affected by L (60.3% and 51.9% of the SS TRMT, respectively). However, the SS G×L /SS G×Y ratio was 38-fold for CP and sevenfold for CT, and for this reason the GGE biplot with the polygon view based on the G×L interaction could be useful to study the pattern observed. Across the ten locations, the closest to the "ideal" for SY and stability was G1, followed by G2; however, G1 was the most stable ( Figure 3a). The "ideal" locations in terms of representativeness and discriminating ability were E3 and E4, followed by E10 ( Figure  3b).
(a)  Across the ten locations, the closest to the "ideal" for SY and stability was G1, followed by G2; however, G1 was the most stable ( Figure 3a). The "ideal" locations in terms of representativeness and discriminating ability were E3 and E4, followed by E10 ( Figure  3b). . GGE biplot based on the seed yield for (a) the ranking of five lentil genotypes against the "ideal", (b) the ranking of the ten locations (abbreviated, E1-E10) against the "ideal" (the ideal is the most representative and discriminating) for two growing seasons (2019 and 2020).

Protein and Cooking Time
The year (Y), G, L, and the two-way (Y × L, G × Y, G × L) and three-way interactions (Y × L × G) were highly significant (p ≤ 0.01) for CP, CT and CPY ( Table 3). The cooking time was genetically controlled, with G explaining 38.6% of the SSTRMT, whereas CP and CPY were mainly affected by L (60.3% and 51.9% of the SSTRMT, respectively). However, the SSG×L/SSG×Y ratio was 38-fold for CP and sevenfold for CT, and for this reason the GGE biplot with the polygon view based on the G×L interaction could be useful to study the pattern observed.
G2 was the best performing in CT (23.05 min) over the years, and for each growing season ( Table 5). The GGE biplot also indicated that G2 was the single "winner" genotype ( Figure A5; Appendix A). The ideal locations for CT were E3 and E5, followed by E2 (Figure A6; Appendix A). The mean CPY ranged from 0.20 (E10) t ha −1 to 0.51 t ha −1 (E3, E8), with an average of 0.36 t ha −1 across the locations ( Table 5). The high-yielding genotypes or locations exhibited a respectively high CPY. Regarding the GGE biplot analyses, the "ideal" locations, exhibiting high discriminating ability and representativeness for the selection of superior genotypes and for the examined traits, the rankings were as follows: for SY E3, E4 > E10; for CP E1, E3 > E7; and for CT E3, E5 > E2 (Figure 4).   Figure 3. GGE biplot based on the seed yield for (a) the ranking of five lentil genotypes against the "ideal", (b) the ranking of the ten locations (abbreviated, E1-E10) against the "ideal" (the ideal is the most representative and discriminating) for two growing seasons (2019 and 2020).
G5 and G2 were significantly the highest in CP (26.23% and 26.20%, respectively), whereas G1 exhibited the lowest value (24.76%, Table 5). The GGE biplot also indicated that G2 was the best performing genotype in E4, E6 and E10, whereas G5 was better in E1, E3, E7, E8 and E9 ( Figure A3; Appendix A). The "ideal" locations for CP were E1 and E3, followed by E7 ( Figure A4; Appendix A). Table 5. Comparisons of the means for the crude protein concentration (CP), cooking time (CT) and crude protein yield (CPY) of the five lentil genotypes across the ten locations over two years (2019-2020). G2 was the best performing in CT (23.05 min) over the years, and for each growing season ( Table 5). The GGE biplot also indicated that G2 was the single "winner" genotype ( Figure A5; Appendix A). The ideal locations for CT were E3 and E5, followed by E2 ( Figure A6; Appendix A). The mean CPY ranged from 0.20 (E10) t ha −1 to 0.51 t ha −1 (E3, E8), with an average of 0.36 t ha −1 across the locations ( Table 5). The high-yielding genotypes or locations exhibited a respectively high CPY. Regarding the GGE biplot analyses, the "ideal" locations, exhibiting high discriminating ability and representativeness for the selection of superior genotypes and for the examined traits, the rankings were as follows: for SY E3, E4 > E10; for CP E1, E3 > E7; and for CT E3, E5 > E2 (Figure 4).

Discussion
Knowledge of the relative contribution of the genotype, location and growing season, and their interaction on SY, yield components, 1000 SW, CP and CT is important for the implementation of breeding programs and guidance for local farmers and private companies.
Due to high concern about the higher productivity of grain legumes, environmental sustainability, healthy diets of plant-based proteins and food/feed security, there is increasing interest about protein crops like lentils [25].

Seed Yield, Yield Components and PH
The combined ANOVA for the SY, PP and 1000 SW of the five genotypes across the 10 locations for the two growing seasons (2019 and 2020) revealed highly significant single treatment effects (L, G) and their combinations (two-way Y × L, G × Y, G × L and threeway Y × L × G interactions), except for the Y. The treatment effects were less consistent for PH (Table 3).
Dealing with the G × E interaction in breeding programs is crucial, and could be managed in three ways [55]: it can be ignored, reduced or exploited. When the G × E is ignored, genotype recommendations are based on the mean performance across all of the testing environments. In the other two cases, stability analysis is required as a step where a decision has to be made regarding the division into sub-groups. Then, a genotype recommendation is made separately for each subgroup (reduction) or for particular environments (exploitation). Based on the SY data of the 2019 ( Figure A1; Appendix A) and 2020 growing seasons ( Figure A2; Appendix A), the two annual GGE biplots of the five lentil genotypes in the ten locations (E1-E10) were constructed with the polygon view. The polygon view of the biplot presents which genotypes performed the best in one or more locations, thus enabling the identification of mega environments [18,20]. Genotypes G1, G2, G3 and G5 in Figure A1 (Appendix A), and G1, G2, G3 and G4 in Figure A2 (Appendix A) are located at the top of the polygons. These genotypes are the highest or the lowest yielding genotypes in some or all of the locations, as they are located at the maximum distance from the biplot center. However, the annual polygon view reveals the absence of a unique mega-environment in the tested area. The annual biplots indicated a non-repeatable location grouping pattern and inconsistency in the performance of the genotypes, which is not suggestive of the existence of different mega-environments. For example, G3 (cv. Thessalia), which is located on the top of the polygons in both growing seasons (2019, 2020) was the winning genotype in different locations in each growing season, and there was not a repeatable grouping over the years. Then, it would be more suitable to select for superior lentil genotypes based on the combination of high average SY and SY stability.
Across Y and L, the mean SY among the genotypes ranged from 1.37 t ha −1 to 1.73 t ha −1 , and fell within the SY range of 0.78 t ha −1 to 2.89 t ha −1 , in compliance with previous studies in Mediterranean conditions [56]. G1 ranked in the top for SY and PP, indicating a highly productive and stable genotype across the locations and years. On the other hand, G4 was a large-seeded genotype exhibiting significantly higher 1000 SW, which ranked at the bottom, with 21% less SY than G1, corroborating previous studies about the lower yield potential of large-seeded genotypes compared to their small-seeded counterparts [57]. The negative correlation between SY and 1000 SW could be attributed to the lower PP of large-seeded genotypes [7].
The large SS for L (54.3%) indicated that locations were diverse and that differences among location pedoclimatic conditions caused most of the variation in SY which was widely ranged from 0.91 t ha −1 to 2.27 t ha −1 . Data across Y and G, indicated that E8 (Komotini-Mesochori) and E3 (Orestiada) (both in the same Region of Thrace) were the most productive locations, with mean SYs of 2.32 t ha −1 and 2.27 t ha −1 , respectively (Table 4). Seed yield is one of the most variable traits, and even though it is genetically controlled by many genes with small effects, it is also manipulated by many agronomic factors and environmental variables, such as soil characteristics and climatic conditions. In location E8, the higher SY may be attributed to the highest mean PP (77.6 pods plant −1 ), whereas in E3 it may be attributed to the highest mean 1000 SW (48.5 gr) (Table 4), which is one of the most important parameters related to ultimate crop yields. These results are in accordance with Tyagi and Khan [58], who concluded that PP is the most important trait responsible for SY in lentil. Their correlation studies indicated that PP was positively and significantly correlated with SY at both the phenotypic and genotypic levels. Similar results were also recorder by Salehi et al. [59] for twenty lentil genotypes under normal and drought stress conditions, whereas a recent field evaluation of high yielding lentil genotypes revealed a highly significant correlation between agronomic parameters (number of branches plant −1 , seeds pod −1 and 1000 SW) and the SY of lentil [60].
The pedoclimatic conditions recorded in these locations were also favorable for sustaining high lentil SYs. From the flowering to pod filling stages, the amount of precipitation (>105 mm) in locations E3 and E8 (April-May precipitation) may have secured adequate soil moisture for the proper growth of the lentil genotypes. A deficiency of water during any growth stage in lentil, as in any legume species, often results in a loss of SY [61]. More particularly, the seedling and flowering stages are the most sensitive to water availability and drought stress [59]. Moreover, the soil types classified in these locations revealed very fertile soils which belong in the orders of Fluvisols (E3) and Luvisols (E8) ( Table 2) (soils with intrinsic high fertility), typically with a high cation exchange capacity (CEC) and high levels of clay (37-49% in E3 and 34-47% in E8, data not shown). These soils are well-drained (100-150 cm) (E8) or moderately well-drained (50-100 cm) (E3), exhibit an absence of a cultivation pan, show good soil structure and consistency, and show no (E3) or weak (E8) erosion (data not shown, derived from classification map symbols). Conversely, over Y and G the lowest SY (0.91 t ha −1 ) was recorded in E10 (Petrana), a location which is characterized by calcareous soils (the presence of calcium carbonate from 30-39%; Table 2), which may be rather problematic under certain conditions for lentil crop production (i.e., the low availability of soil P and Zn).
With regards to location representativeness and discriminating ability, E4, followed by E3 and E10, showed the maximum proximity to the ideal environment for selection (Figure 3b), and therefore these locations were suggested as suitable locations for lentil breeding when targeting high SY. It is noteworthy that the locations suggested for breeding are not always the most productive (e.g., E10), but those in which phenotypic differences are maximized and thus selection is more efficient.
The location effect was the main source of variation for PP (65.4%), while the PP values ranged from 25.3 pods per plant −1 to 77.6 pods per plant −1 among the environments. A high significant correlation (r = 0.97, p < 0.01, Table 6) of SY and PP showed that PP is a reliable criterion to select lentils for seed yield, confirming previous results [58][59][60]. The location and genotype effect were similar for PP and PH, while their correlation was positive and relatively high (r = 0.67, p < 0.05, Table 6).
Seed size is an important seed quality trait that determines the market categories and the relative market value of lentil cultivars. Low L and G × L effect (5% and 5.7%, respectively) were detected for the 1000 SW, which was mainly affected by G, as it explained 82.2% of the total variation. Therefore, the selection of the appropriate cultivar must be the criterion for seed size, as this trait is independent of the growing environment effect.
The tallest plants across Y and G were recorded in E2 (Ypato: the southernmost location among the tested locations) (PH = 39.5 cm) ( Table 3). This was the location with the highest precipitation and mean temperature among the test locations in both seasons. These conditions may result in prolonged growth seasons and the accumulation of large amounts of vegetative assimilates. However, although various parameters of lentil, including the PH, are influenced by both the genetic and the environmental factors, the changes in the genetic architecture seem to be the major determinant of PH [60].

Crude Protein
Lentils are an excellent source of plant proteins and an alternative to animal and soy protein for the food and feed industries [29]. The CP values ranged between the genotypes from 24.76% to 26.23%, while the red lentil G5 ranked first, followed by G2. In general, a higher seed protein content is observed in lentils with red cotyledons. Subedi et al. [26] tested 34 genotypes in 10 locations and found 3% higher protein content in red lentils compared to other categories. The GGE biplot analysis illustrated that G5 had higher and consistent CP values, and adapted better in five locations (E1, E3, E7, E8, E9), while G2 showed better adaptation in three (E4, E6, E10) ( Figure A3; Appendix A). A clear cut explanation as to why these two genotypes exhibited a higher CP in these specific locations cannot be given. However, even though these locations exhibit diverse pedoclimatic conditions, some specific soil traits may partly explain these results. For example, other locations (i.e., E1, E7) exhibited sufficient soil P content, which may have favored the genotypes' N-fixing ability, and other locations (i.e., E3 and E8) have soils with intrinsic high fertility (Fluvisols, Luvisols), high cation exchange capacity (CEC) and high levels of clay. Moreover, most of these locations contained medium or adequate quantities of CaCO 3 , which accounts for an important fraction of the C present in soils, linking the longterm geological C cycle with the faster biogeochemical cycling of soil organic carbon [62]. These conditions may have favored the mobilization of N to the seed, and thus may have increased the seed CP content.
The location explained 60.3% of the total variation, and was the main source of variability. Similar results were reported by Subedi et al. [26], although other researchers found low variation across locations [63]. The fact that the CP was highly affected by the location could be attributed to the low diversity between the tested genotypes, pedoclimatic differences between the locations, and the location effect on the N-fixing ability [64]. Moreover, CP is a quantitatively inherited trait [65], and thus it is strongly influenced by the environmental fluctuations [26,66]. In grain legumes, environmental variability affects the pod filling period, which is responsible for the increasement of seed CP [67].
The CP values ranged among the locations, from 23.05% to 28.87%. The higher CP was recorded in E6, which was characterized by very low precipitation during the crucial period of pod filling (April-May) in both years. It is possible then that the water deficit negatively impaired the starch synthesis, resulting in an increase of seed protein.
The mean CPY ranged from 0.20 to 0.51 t ha −1 , with an average of 0.36 t ha −1 across the locations. This mean value was similar to the one reported by Lizarazo et al. [63] for Northern Europe regions (0.4 t ha −1 ), but was 38% lower than the one reported by Subedi et al. [26] in the USA. CPY followed the ranking of locations for SY. Therefore, E3 and E8 were the most productive locations for both SY and CPY. Locations E3 and E8, which were the most productive, were located in the same region (Thrace), and had similar elevation, soil properties and precipitation during the crucial period of anthesis and pod filling ( Table 2).
The protein content and SY showed a negative correlation that was in accordance with other studies [31,68]. This indicates an independent selection for the two characteristics that makes a simultaneous selection for both traits very difficult, as the rate of genetic gain in one trait is reduced by the other. In this case, independent culling or truncation selection could be the most suitable breeding method for simultaneous selection for SY and CP [69]. The ideal environment for CP breeding was E1 and E3, followed by E7.

Cooking Time
Cooking time is an important commercial trait for pulse crops, which is characterized by high heritability values [30,31,68,70,71]. In the present study, genotype was the main contributor to CT variance (38.6%), whereas significant differences were detected among the genotypes. The mean CT across the years and locations ranged from 23.05 min to 28.03 min between the genotypes. G2 had the lowest CT and showed considerable stability across the locations and growing seasons. In comparison, G2 required 12% less time to reach the optimum cooking point than the second genotype in the ranking list and 17.7% less time than the genotype which ranked in the bottom. Other studies reported significant genotypic variability for CT [34]. More specifically, Erkine et al. [31] reported an average CT of 33 min in a collection of 24 genotypes, Jood et al. [72] reported values ranging between 38 min and 43 min, and Vandenberg [73] recorded CTs varying from 15 min to 20 min. The highest CT was observed for the red lentil G5. The same genotype was evaluated for CT by Theologidou et al. [74], and was characterized as unstable and highly affected by environmental fluctuations. Generally, red lentils are characterized by poor cooking quality, and this is the reason why they are consumed dehulled [75], as dehulling reduces their CT to as low as 15 min [76]. However, it should be underlined that the CT of G5 was commercially acceptable as it was less than 30 min and was comparable to cultivars like G3, which is categorized as a fast-cooking genotype among Greek varieties [77].
The location effect explained 19.1% of the total variation. The mean CT across the years and genotypes ranged between locations from 24.96 to 29.58 min. The lower CT was recorded in E9, followed by E1. Both locations pertained to the BSk climate category (arid, steppe, cold); however, their soil characteristics and altitude were quite different. Locations like E3 were expected to have high CT because of their very low precipitation [77]. Ultimately, a specific pattern linking climatic and soil conditions to CT was not identified, and this could be attributed to the significant G × L interactions.
However, it is remarkable that the range of mean CTs between the genotypes or locations was near 5 min. This CT range is not judged to be high according to the market criteria, and was attributed to the fact that most of the tested genotypes were commercial varieties that have already been screened for low CT before their release into the market.
A significant negative correlation was observed for CT with SY. Ninou et al. [7] investigated the indirect effect on the seed quality characteristics of an intense lentil breeding project for high yield [78], and reported the selection of higher-yielding plants with reduced CT compared to the parental landraces. The relatively high percentage of G in the total variability suggested that the phenotypic differences observed for CT between genotypes are highly controlled by genetic factors, and that phenotypic selection could be an effective means to develop fast-cooking lentil varieties. GGE-biplot diagrams ( Figures A5 and A6; Appendix A) illustrated that E3 and E5, followed by E2, showed the highest representativeness and discriminating ability for CT, and thus could be suggested as the most suitable locations for breeding for CT. Regarding the correlation of CT with seed size, expressed by 1000 SW, our findings showed a negative correlation. The references in this issue are contradictory, as Erkine et al. [31] found a high positive correlation (r = 0.919) between CT and seed size, and suggested that seed size can be used to predict CT, whereas Theologidou et al. [74] and Ninou et al. [7] concluded the opposite results.

Market and Breeding Implications
A significant diversification and specialization of locations was identified with regards to the traits of major interest (SY, CP, CT). As a consequence, different locations were found to be suitable for each trait (Figure 4). The purpose of the cultivation was the criterion that would determine the selection of the location in which the lentil will be cultivated. As such, when the purpose of the cultivation was a high SY, disregarding the CP content or CT, locations E8 and E3, followed by E7, were the best. This is important for seed companies, or when high lentil production is needed for food or feed products that are not related to CP content or CT. In the case where the purpose of cultivation was the production of lentil varieties with a high CP content, the most suitable locations were E6 followed by E9 and E2. These locations could be selected for cases where high protein food or feed are required. The best quality for CT was found in lentil varieties cultivated in location E9, followed by E1. Those locations were considered suitable for lentil production directed for wide human consumption, or to niche markets that are interested in lentils certified as PGI (Protected Geographical Indication) or PDO (Protected Designation of Origin). It should also be emphasized that E9 was the optimum location for lentil production of high quality, as it combined optimum performance for both CP content and CT.
Breeding must be applied in locations that are characterized by a high discriminating ability and representativeness. Therefore, when the purpose of cultivation is the selection for SY, the "ideal" locations were E4 and E3, followed by E10, whereas when the purpose was high CP, the "ideal" locations were E1 and E3, followed by E7. As far as CT is concerned, the "ideal" locations were E3 and E5, followed by E2. The above clearly indicates that location E3 was serving the three purposes, as it was "ideal" or close to "ideal" for the selection of superior genotypes for SY, CP and CT, and therefore is suggested as the optimum location for lentil breeding.
A comparative investigation of the general performance of the five genotypes for the three basic traits (SY, CP, CT) revealed that G2 had a genotypic profile characterized by the top CP and CT values, and high SY, while G1 was defined as the higher yielding cultivar for SY, with very good CT across the ten diverse test locations. Both cultivars indicated a wide adaptation capacity, and thus are suggested for cultivation in many regions. It is noteworthy that G2 had the highest CP content and CPY, and it could be reclaimed by food or feed industries that are interested in alternative high-protein products. This is even more important because, nowadays, progress in food processing technologies has opened new avenues for the application of lentil protein in novel food formulations [25]. Finally, G2 is considered a promising genetic material to be included in future breeding programs as a donor parent for quality and productive traits, whereas G1 would be suitable for seed production companies.

Conclusions
(i). Significant G × L interactions were detected for SY, CP and CT. Different locations were found to be suitable for each trait. The purpose of the lentil cultivation is the key point that will determine the location in which the crop will be established. (ii). Breeding must be applied in locations that exhibit a high discriminating ability and representativeness, regardless of their productivity potential. Location E3 (Orestiada, Region of Thrace) was identified as "ideal" or close to "ideal" for the selection of superior genotypes for SY, CP, CT and CPY, and is recommended as the optimum location for breeding lentils for major agronomic and quality traits. (iii). Suitable genetic material was identified either for cultivation by lentil growers or for further improvement by breeding programs. Genotype G2 (cv. Dimitra) showed the top CP, CT and CPY values along with high SY, whereas G1 (cv. Samos) was the most productive cultivar for SY, with very good CT and high CPY. Both cultivars indicated a wide adaptation capacity, and thus are suggested for cultivation in many regions.  Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to privacy restrictions.   Figure A2. GGE biplot with the "which-won-where" pattern based on the seed yield of the five lentil genotypes (abbreviated G1-G5) across the ten locations (abbreviated E1-E10), based on the 2020 growing season. Figure A2. GGE biplot with the "which-won-where" pattern based on the seed yield of the five lenti genotypes (abbreviated G1-G5) across the ten locations (abbreviated E1-E10), based on the 202 growing season.