Grain yield and genotype x environment interaction in bean cultivars with different growth habits

ABSTRACT Breeding of common beans (Phaseolus vulgaris L.) shows restrictions in the genetic advance because of the effect of the environment. Therefore, the behavior of the yield components of genotypes varies according to the crop’s environment. The genotype x environment interaction can cause genotypes with high yields in one location not to behave in the same way in other localities, which limits the recommendation of cultivars for different environments. The objective of this research was to evaluate agronomic traits in new improved bean cultivars in high tropic environments, as well as to determine which cultivars show phenotypic stability for yield. Multi-environment tests were carried out during 2016 and 2017 in two regions of the department of Cundinamarca, in Colombia. Significant differences were found for the genotype x environment interaction and highly significant differences for the evaluation environments and genotypes. The greatest variation was attributed to genetic effects, followed by environmental effects and the genotype x environment interaction. The first two principal components for grain yield showed 88.86% of the variation of the genotype x environment interaction. Cultivars Serrania and Sutagao, of climbing growth habit, were identified as stable and with high yield potential, so they can be considered as a commercial alternative for bean growers. Graphical abstract


Introduction
The common bean (Phaseolus vulgaris L.) is considered one of the most important legumes for human consumption in the world because of its contents of proteins, carbohydrates, vitamins, and minerals (Broughton et al., 2003;FAO, 2018). Although Colombia has all the conditions for its cultivation, the country imports up to 35,000 t to meet the national demand (Fenalce -Federación Nacional de Cultivadores de Cereales y Leguminosas, 2017). Around 120,000 families are linked to the bean production process in the country (Ligarreto et al., 2017) with a planted area of 110,000 ha per year and national average yields of 1.39 t ha −1 (Fenalce -Federación Nacional de Cultivadores de Cereales y Leguminosas, 2017).
In Colombia, bean yield went from 0.945 t ha −1 in 1998 to 1.076 t ha −1 in 2000, and an increase of 0.32 t ha −1 was reported for 2014 (FAO, 2018). According to Fenalce -Federación Nacional de Cultivadores de Cereales y Leguminosas (2017), national production has been maintained despite the reduction in the bean planting area due to the increase in yields per hectare in some high Andean zones of the country. Additionally, Ligarreto et al. (2017) report that yields of 3.2 t ha −1 and 1.7 t ha −1 are achieved for beans with climbing and shrubby habits with the release of new high-yielding cultivars, which contributes to meeting the food needs of the population.
Yield components in crops are closely influenced by the genotype x environment (GxE) interaction (Elias et al., 2016). Yield is a complex variable determined by direct and indirect factors. Among the direct factors are the number of pods per plant and grain weight, which are low heritability traits and their expression is highly affected by the environment (Ligarreto et al., 2017). For cultivar development, it is customary to work in optimal environments for crop growth, which are usually specific. On the other hand, cultivars face different conditions such as planting systems, fertilization levels, pest and disease management, and climatic and edaphic conditions when delivered to producers. These factors impact yield and cause variations between locations or environments (Rodríguez-González. et al., 2011).
GxE interaction can cause genotypes with high yields in one location not to perform well in other localities. This limits the recommendation of cultivars for different environments and promotes the selection of genotypes per environment (Pérez et al., 2005). That is why the evaluation and selection of high-performance and highstability genotypes are important in breeding programs to identify cultivars with general and specific adaptability (Rodríguez-González. et al., 2011) and to provide growers with superior and competitive cultivars at the production level (Nascimento Filho et al., 2010). Furthermore, tests in various environments facilitate the identification of mega-environments (Crossa & Cornelius, 1997;Crossa et al., 2002) to determine the general or specific adaptation of genotypes.
There are various statistical methodologies for the analysis of the GxE interaction, ranging from univariate statistics to multivariate methods. One of the most used multivariate methods is the additive main effects and multiplicative interaction (AMMI) model that allows estimating the stability of cultivars and classifying environments (Carbonell et al., 2004;Shahzad et al., 2019). The model combines the analysis of variance as an additive model with principal components (PC) as a multiplicative model that allows analyzing interactions (Rodríguez-González. et al., 2011). The other method is the site regression (SREG) model that allows studying the incidence of multiple factors on crop behavior (Valencia & Ligarreto, 2010) and is more useful when evaluating environments to know their power of representativeness or genotype discrimination (Rodríguez et al., 2012). Both methods use a graph or biplot that facilitates interpretation (Burgueño et al., 2002;Yan et al., 2000).
The objective of this research was to evaluate agronomic traits of new improved bean genotypes of climbing and shrubby habits in five environments of dry and humid tropics in the high Andean zone of Colombia, as well as to determine which genotypes show high phenotypic stability for yield.

Materials and methods
Multi-environment trials of bean cultivars (genotypes) were carried out during 2016 and 2017 in the regions of Ubate and Guavio in the department of Cundinamarca, in Colombia. The region of Ubate includes approximately 15% of the 6,500 ha that are sown per year in the Department of Cundinamarca. Therefore, it is considered one of the main production sites from areas of less than 3 has sown with climbing beans, whereas in the Guavio region, the crop is cultivated in areas smaller than 1 ha. The environments were identified according to the municipality or locality in which the trials were established. For the region of Ubate, trials were carried out in the municipality of Simijaca at an altitude of 2554 masl (5°30ʹ36.3" N and 73°51ʹ28.0" W) and in the municipality of Guacheta at an altitude of 2563 masl (5°20ʹ45.0" N and 73°43ʹ00.1" W), during February and August 2016. This season coincides with the period between the onset of rainfall in February and that of the dry period in August in which farmers plant bean crops. In the region of Guavio, trials were established in the municipalities of Gachala at an altitude of 1800 masl (4°39ʹ45.37" N and 73°30ʹ58.28" W) and in the municipality of Gama at an altitude of 2350 masl (4°44ʹ22.60" N and 73°35ʹ36.19" W), from August 2016 to February 2017. Rainfall starts to reduce during August and continues to decrease until February. During this period farmers plant bean crops. A fifth environment was established under greenhouse conditions in the Bogota plateau (BP), in the Faculty of Agricultural Sciences of the Universidad Nacional de Colombia, at an altitude of 2620 masl (4° 38ʹ08" N and 74°04ʹ58" W), between January and June 2016 (Table 1).
Imetos 300 meteorological stations were installed in the municipalities of Simijaca, Guacheta, and Gachala. Variables were monitored and recorded in real time using the web-based software Fielclimate® (www.fielcli mate.com, Pessl Intruments©, Weiz, Austria). For the Gama and BP localities, the information was recorded using a temperature and humidity data logger (DT172, CEM, Shenzhen, China) and wireless rain gauges (RGR126, Oregon Scientific, Oregon, USA).
Seven common bean cultivars were used, two shrubby-type cultivars (habit I) called Bacata (with red grains) and Bianca (with white grains) and five climbing or indeterminate habit cultivars (habit IV) named Hunza, Chie, Iraca, Sutagao (red seeded type) and Serrania (cranberry type) cultivated in Colombia. These cultivars were developed by the Edible Legumes program of the Faculty of Agricultural Sciences at the Universidad Nacional de Colombia and constitute the genetic supply for bean cultivation in the farms of local growers. The seeds were sown in plots with four rows of 7 linear meters. The plant spacing for shrubby varieties was 16 cm between plants and 1 m between rows, for a density of 62,500 plants ha −1 according to the distances used by producers in the region of Ubate. In this region, growers sow relay strip intercropping systems using other species such as vegetables, so that when a crop is close to harvest the new crop species is sown interspersed. In the region of Guavio, sowing is carried out using high furrows to facilitate the drainage of excess rainfall. The plant spacing for the cultivars of climbing habit was 30 cm between plants and 1.1 m between rows, for a density of 30,300 plants ha −1 according to the usual practices carried out in climbing bean monoculture in staking systems in Colombia.
Thirty days before sowing, pH was corrected to 5.3 with a dose of 5.70 kg/plot (2 t ha −1 ) of dolomite lime in the Gama, Gachala, and BP environments. Thirty days after emergence (DAE), when weeding and hilling were carried out, the plants were fertilized with a mixture of the following soil fertilizers: i) Diammonium phosphate (DAP, Ecofértil S.A., Colombia) at a dose of 1.10 g/plant (70 kg ha −1 ), ii) urea (Ecofértil S.A., Colombia) at a dose of  1.20 g/plant (80 kg ha −1 ), and iii) potassium chloride (KCl, Ecofértil S.A., Colombia) at a dose of 0.80 g/plant (50 kg ha −1 ). The composition of each fertilizer was as follows: DAP: 18% total nitrogen (18.0% ammonium N) and 46% assimilable phosphorus (P 2 O 5 ); urea: total nitrogen (46.0% urea) and 1.5% maximum Biuret; and KCl: watersoluble potassium (60% K 2 O) and chlorine (45% Cl). The evaluation of the number of pods per m 2 (which corresponded to three plants for climbing-habit beans and six plants for shrubby-type beans per plot), number of grains per pod, and weight of 100 grains was carried out in 10 competing plants within each plot of the bean genotypes assessed. The dry grain yield (based on the two central rows) was evaluated at the end of the crop cycle at 120 DAE for the shrubby cultivars in the Gama, Gachala, and BP environments, and at 135 DAE in Simijaca and Guacheta. Climbing cultivars were harvested at 155 DAE in the Gama, Gachala, and BP environments and at 180 DAE in the other two environments. A randomized complete block design was used with four replicates per treatment. Analyses of variance were performed to determine the effect of the GxE interaction on the studied variables. GxE interaction and the analyses of phenotypic stability and identification of the more discriminatory and representative environments were performed using biplots of the AMMI and SREG models. The R statistical software (RStudio, 2014) was used through an application called GEA-R (genotype x environment analysis with R for Windows) Version 4.0. (Pacheco et al., 2015) developed by the International Maize and Wheat Improvement Center (CIMMYT).
The AMMI stability value (ASV) was calculated for each genotype according to the relative contributions of the principal component axis scores (IPCA1 and IPCA2) to the interaction sum of squares.
The AMMI stability value (ASV) as described by Purchase et al. (2000) was calculated as follows: ASV ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi where IPCA1sq IPCA2sq is the weight derived from dividing the sum of IPCA1 squares by the sum of IPCA2 squares (from the AMMI analysis of variance table). The larger the IPCA score is, either negative or positive, the more adapted a genotype is to a certain environment. Lower ASV values indicate a more stable genotype across environments (Purchase et al., 2000).
The yield stability index (YSI) was also calculated using the sum of the ranking based on yield and the ranking based on the AMMI stability value.
where RASV is the rank of the genotypes based on the AMMI stability value, and RY is the rank of the genotypes based on the average yield across environments. The genotypes with the lowest YSI values are the most desirable, with higher yield and stability (Bose et al., 2014).

Yield components and bean grain yield by environments and genotypes
Highly significant statistical differences (P < 0.01) were observed between the five evaluated environments of the seven improved bean cultivars for the variables number of pods per m 2 , number of grains per pod, weight of 100 grains (g), and grain yield (t ha −1 ). Highly significant differences were also found between cultivars and significant differences were observed for the GxE interaction (Table 2). Table 3 shows the results per environment. The comparison of means shows that there were significant differences. Regarding the number of pods per m 2 , Guacheta and Simijaca were the most favorable environments with mean values of 107.96 and 113.13, respectively, while BP was the least favorable environment with 50.64 pods per m 2 . For the variable grain yield, the best environments were Simijaca and Gama with 3.34 and 3.12 t ha −1 , respectively, which showed significant differences compared to the environments of Guacheta, Gachala, and BP. The weight of 100 grains was higher in Gama (84.30 g) compared to the other environments, followed by Gachala (79.70 g) and Simijaca (74.30 g), but without significant differences with BP (73.50 g).
Regarding the number of grains per pod, the best environment was Gachala with a mean value of 6.00, followed by Simijaca, Gama, and BP with mean values of 5.70, 5.60, and 5.50, respectively, but without significant differences between them. The least favorable environment was Guacheta with 5.30 grains per pod.
Climbing genotypes showed the higher number of pods per m 2 , with values exceeding 100 pods per m 2 . Shrubby bean genotypes Bacata and Bianca reached mean values of 56.27 and 62.58 pods per m 2, respectively, with significant differences compared to the climbing genotypes (Table 3). Bacata and Bianca obtained a lower number of grains per pod close to 4.03 whereas climbing genotypes Hunza and Iraca obtained mean values of 5.33 and 8.03, respectively. Regarding grain size, Bacata stood out among the shrubby bean genotypes with 73.72 g for 100 grains, and Hunza stood out among the climbing genotypes with 89.00 g for 100 grains, showing significant differences with the other genotypes (Table 3). The general mean number of pods per m 2 was 92.96 and the number of grains per pod was 5.60. The weight of 100 grains was 76.70 g and the average grain yield was 2.80 t ha −1 for all genotypes (Table 3).
The genotypes with the highest yield were Iraca and Serrania with 3.42 and 3.30 t ha −1 , respectively, followed by Hunza, Sutagao, and Chie, whereas the genotype with the lowest mean yield was Bacata with 1.77 t ha −1 ( Table 3). The mean grain yield for each genotype in the localities of Guacheta, Simijaca, Gama, Gachala, and BP is shown in Table 4. Significant differences were observed for the genotypes in the five locations.

Analysis of GxE interaction of the AMMI model for grain yield
The main effect of the genotype (G) accounted for 53.14% of the sum of squares of the treatments and the main effect of the environment (E) only accounted for 15.87% of the variation, whereas the remaining 30.99% was due to the GxE interaction (Table 5). The  first two principal components (IPCAs) used for the GxE biplot were statistically significant (P < 0.01) (Figure 1). These IPCAs explained 88.86% of the variation of the GxE interaction. IPCA1 explained 62.03% and IPCA2 26.83%. Cultivars Sutagao and Chie behaved better regarding the grain yield trait in the Simijaca and Guacheta environments. For the Gama environment, the cultivar with the best behavior was Hunza, whereas Bianca showed a correlation close to 1 with the BP and Gachala environments. On the other hand, Iraca and Serrania were not found in the same quadrant with any other environment (Figure 1).
Additionally, genotypes 'Hunza', 'Chie', 'Iraca', 'Sutagao', and 'Bacata' were located on the vertices of the polygon, which corresponded to good behavior in the environments near the respective vertex. The genotypes with the lower contribution to the GxE interaction were 'Serrania', and 'Bianca', and Guacheta and BP were the environments that contributed the lower to this interaction. There is a possible lack of adaptation of 'Bacata' to the Simijaca and Guacheta environments, as well as a contrast between the environment of Gama and Gachala and BP.  Regarding the environments, BP and Gachala were located in the same quadrant with vectors of similar length, as observed for Simijaca and Guacheta. Gama was located alone in one quadrant and its vector showed a greater length. Additionally, it had an angle close to 180° regarding genotypes Iraca and Serrania. The same was observed for the BP and Gachala environments and genotypes Chie and Sutagao (Figure 1). The trend and relationships described above are consistent with the results in the AMMI biplot of the variables number of pods per m 2 , number of grains per pod, and weight of 100 grains. The mentioned biplots are not shown in this paper.
The additive main effect and multiplicative interaction (AMMI) stability value (ASV) ranked the genotypes based on the lowest score (Table 6). Lower scores represent more stable genotypes. Based on the ASV, the most stable genotype for yield was Sutagao since it obtained the highest ASV ranking. Hunza was ranked as the least stable because it showed the highest ASV. The sum of the yield and stability rankings (YSI) ranked Sutagao and Serrania as the genotypes that combined high yield with stability. Although Hunza and Chie were high-yielding genotypes, they were unstable because of their low rank according to the YSI. RY = ranking of grain yields in all environments, IPCA1 score = interaction principal component axis one score, IPCA2 score = interaction principal component axis two score, ASV = AMMI stability value, RASV = ranking of the ASV, and YSI = yield stability index.

Analysis of the GxE interaction of the SREG model for grain yield
Similar results to those of the AMMI model were obtained using the SREG model. The main effect of the environment accounted for 15.85% of the sum of squares of the treatments, and the main effect of the genotype explained 53.14% of the variation, whereas the GxE interaction accounted for 30.99% (Table 5).
The genotypes are grouped in the same way as in the AMMI model. Cultivar Hunza continued to be associated with Gama, just as cultivars Chie and Sutagao show positive correlations when they are close to Simijaca and Guacheta. However, in this case, it was possible to observe a clearer positive correlation compared to that in AMMI with Gama. On the other hand, Bianca and Bacata modified their position when being in a quadrant without any environments while showing negative correlations within the five environments ( Figure 2).
The first two principal components were highly significant and accounted for 91.59% of the total variation (69.85% of IPCA1 and 21.74% of IPCA2). The environment with the highest interaction was Gama and the genotypes with the highest interaction were 'Hunza', 'Iraca', 'Bacata', and 'Bianca', which correspond to the vertices of the polygon (blue dotted line) ( Figure 2). Additionally, cultivars Hunza, Chie, Sutagao, Iraca, and Serrania showed positive interactions with the different environments. The red seeded type cultivars Sutagao, Chie, and Hunza were grouped in the Gama environment and have angles of less than 90° compared to Guacheta and Simijaca. 'Iraca' (red seeded type) and 'Serrania' (cranberry type) were grouped in the BP and Gachala environments. 'Bacata' and 'Bianca' were located far from the five environments ( Figure 2). The proximity of the vectors belonging to the Guachetá and Simijaca environments, as well as Bogota Plateau (BP) with Gachala, were observed.

Discussion
In the mean comparison tests and GxE stability analyses, the environments of Simijaca and Guacheta are grouped regarding the grain yield trait; the same was observed for BP and Gachala. This is an expected result since the first two municipalities are part of the region of Ubate with a predominance of dry climate. On the other hand, the last two municipalities are related because their climatic conditions are similar; Gachala is part of the Guavio region under the influence of the Guavio dam, which causes an average relative humidity of 88%. These conditions are similar to those of the greenhouse in the Bogota plateau, which maintains a high temperature and relative humidity despite being located at 2550 masl (Ligarreto et al., 2017). Furthermore, the high rainfall in Gachala is only comparable to the high-water availability in BP, just as the edaphic conditions are similar. On the other hand, Gama is located at 1800 masl with intermediate precipitation compared to the two groups mentioned above.
The differences in grain yield and its components for the cultivars in each environment are due to both, the aforementioned climate conditions and the edaphic conditions. The high concentration of phosphorus in the soil in Simijaca and BP is highlighted; however, the soil pH in BP could affect the nutrient uptake despite the high concentration of this element, as mentioned by Ligarreto et al. (2017). The suitable soil pH for common bean cultivation varies between 5.5 and 6.5, with BP, Gama, and Gachala being the most acidic soils. The high bean yields in Gama may be due to the high organic matter content, good drainage and soil moisture retention. These factors contribute to proper root anchorage and expansion that allow better nutrient assimilation and, therefore, good plant development. Additionally, the high yields can be attributed to the right amount and distribution of rainfall during the production cycle, especially during pod formation and filling. However, the opposite was observed in Gachala where excess rainfall was recorded.
The variation of the genotypes' behavior in the different locations under study and the significant GxE interaction imply that the identification of a superior genotype cannot be carried out considering the average response of the cultivars. Additionally, it is necessary to consider the effect of the environment on the genotypes' behavior (Rodríguez-González et al., 2011). Regarding bean yield, the effect of genotypes was greater than that of the environments and the interaction, which indicates that genotypes respond similarly to environmental changes (García & Hernández, 2004).
Based on the AMMI model, distribution of environments was observed in three of its four quadrants. In the SREG model, environments were distributed into two quadrants, with different genotypes forming the polygon. The existence of the GxE interaction for grain yield was confirmed. The models do not differ in the approximation of the environments because the existing variation depends mainly on the genotype and the GxE interaction. The coincidence between the Guacheta and Simijaca environments, as well as that between BP and Gachala, suggests that some environments can be omitted for further tests, since the results show a high correlation and, therefore, do not show differences between them. Instead, tests in more contrasting environments could be established (Frutos et al., 2014). Similar observations are reported by Khanal et al. (2014), who identify highly correlated environments with little contribution to the selection of bean materials with canning quality for industrial processing and that can be omitted due to costs and time.
The results obtained from the models show similarities due to their nature since the SREG model expresses the response of the genotypes and the GxE interaction. Additionally, the SREG model is close to AMMI due to the lower environmental effect compared to that caused by the genotypes and GxE interaction (Crespo-Herrera et al., 2017;Frutos et al., 2014). From the AMMI model, genotypes Iraca and Bacata are phenotypically unstable and do not show specific adaptability to any of the studied environments. This differs from genotype Iraca's disposition in the SREG model since it exhibits specific adaptability in BP and Gachala. In this case, the SREG model showed greater discriminating power (Valencia & Ligarreto, 2010).
On the other hand, 'Serranía and Sutagao' are the most stable materials or those with the least contribution to the GxE interaction due to their proximity to the origin of the axes (Murphy et al., 2009;Vargas et al., 1999). They also show a high specific adaptation to the environments of Simijaca and Guacheta, which are neighboring and geographically similar municipalities. Regarding grain yield, 'Hunza' performed better in the environment of the municipality of Gama. The high discrimination of 'Bacata' and 'Bianca' in the SREG model may be due to the fact that no specific adaptation to the studied environments is detected because of their shrubby habit with a lower yield than the other five bean cultivars of climbing habit. This can be observed even when the yield potential of these two cultivars is good as they exceed the national average yield of 1.39 t ha −1 by at least 380 kg.
According to Purchase et al. (2000), AMMI does not provide a quantitative measure of stability, so it is necessary to sort genotypes according to their position in terms of yield stability using the AMMI stability value (ASV) and yield stability index (YSI). The AMMI stability value (ASV), derived from the AMMI model, was a simple measure of yield stability that allowed classifying superior genotypes in the following order: Sutagao, Serrania, Iraca, and Bacata. This, in turn, contributes to the understanding of the interaction between genotypes and the environment to make more reliable recommendations to producers.

Conclusions
The fact that five of the seven genotypes showed a yield higher than the mean registered by all genotypes (2.82 t ha −1 ) and that this data was higher than the national average yield, indicates that the improved bean cultivars evaluated are a competitive offer for the productive sector. Serrania and Sutagao can be particularly highlighted due to their high yield potential, specific adaptation, and stability to the environments described in this research, which can contribute to the economy of small growers and the sustainability of family farming. The AMMI and SREG models and the parameters ASV and YSI were highly effective in establishing the cultivars' specific stability and adaptability to environments since the more stable genotypes such as Sutagao do not always obtain the highest yield. However, SREG was more discriminating in ordering cultivars based on the expression of their yield potential.