Forty Brazilian carioca common bean cultivars released from 1970 to 2017 and registered in the National Register of Cultivars of the Brazilian Ministry of Agriculture, Livestock, and Food Supply (RNC-MAPA) were evaluated ( Table 1 ). These cultivars represent a large part of the genetic variability that exists between carioca common bean cultivars grown during the approximately 50 years of history of carioca common bean breeding in Brazil. Seed samples were obtained from the Gene Bank of the Instituto de Desenvolvimento Rural do Paraná (IDR-Paraná), Londrina, Brazil, and subsequently multiplied in order to standardize seed germination. Information about the genealogy of the cultivars is available in Supplementary Table 1 .

Cultivars were evaluated at the Experimental Stations of IDR-Paraná in Londrina and Ponta Grossa during the rainy season of 2017 and in Santa Tereza do Oeste and Ponta Grossa in the dry season of 2018. Soil physical–chemical analyses and other characteristics related to the assessment sites are shown in Table 2 . The experiments were arranged in a randomized complete block design with four replications. Sowing was carried out mechanically, and the experimental plots consisted of four 4-m rows, spaced 0.5 m apart, with 15 seeds per linear meter. Two levels of N top-dressing were used in each environment, totaling eight independent experiments. Experiments under high N were fertilized with 40 kg N ha⁻¹ (urea form) at V₃ stage development, while the experiments under low N did not receive N top-dressing. In all experiments, the crops were fertilized at sowing with 300 kg ha⁻¹ of NPK (04–30–10).

Five uniform and representative plants were collected at physiological maturation (R₉ development stage) from each experimental plot. Seeds and aerial part (stems, leaves, and pod shells) were dried in an oven at 70°C for 72 h to determine later the seeds and shoot dry biomass. These samples were also crushed separately using a Willey MA340 type knife mill (Marconi Laboratory Equipment, Piracicaba, Brazil) to determine seed and shoot N content by the Kjeldahl method (Horwitz, 1975) using a Tecnal TE-0371 digester (Tecnal Scientific Equipment, Piracicaba, Brazil). Seed protein content (Prot, %) was determined by multiplying the percentage of N in the seeds by the factor 6.25 (Horwitz, 1975), while harvest index (HI) was defined by the ratio of shoot dry biomass and total dry biomass (seeds and shoot dry biomass).

Seed yield (SY, kg ha⁻¹ with moisture of 13%) was obtained after manual removal and mechanical threshing of plants from the two central rows of each plot. NUsE components were determined according to Moll et al. (1982). NUpE (mg per g of N absorbed in each plot) was calculated by the ratio between plant total N content and total N fertilizer applied at sowing and top-dressing. NUtE (g of seeds produced per mg N absorbed) was calculated using the ratio between seed dry biomass and plant total N content, while NUsE (in g of seeds produced per g of N applied in each plot) was determined by the product between NUpE and NUtE. The traits Prot, HI, NUsE, NUtE, and NUpE were evaluated only in Londrina and Ponta Grossa during the rainy season of 2017 under high and low N conditions, while the SY was determined in all experiments.

The traits were evaluated by comparing the following models: (i) Full model: considering triple interaction among genotype (G) × nitrogen (N) × environment (E); (ii) Reduced model 1: considering double interaction between G × N; (iii) Reduced model 2: considering double interaction between G × E; and (iv) Null model: considering interaction absence among factors G, N, and E. The Full model follows the mathematical model below:

Where: µ is the overall mean, g_i is the random effect of genotype i, b_j/k/m is the random effect of block j within environment k and within N fertilization m, e_k is the fixed effect of environment k, n_m is the fixed effect of N fertilization m, ge_ik is the random effect of interaction between genotype × environment, gn_im is the random effect of interaction between genotype × N fertilization, en_km is the fixed effect of interaction between environment × N fertilization, gen_ikm is the random effect of interaction among genotype × environment × N fertilization, and ϵ_ijkn ~ N(0, σ²) is the random effect of error associated with each experimental plot.

The marginal a posteriori distributions were performed considering non-informative a priori distributions for all model parameters using R software (https://www.r-project.org/) through the ‘MCMCglmm’ (Hadfield, 2010) package. A total of 1,000,000 values were generated by MCMC (Monte Carlo Markov Chain) process, assuming a burn-in period and thinning interval of 500,000 and five iterations, respectively. The MCMC convergence was verified using the Heidelberger and Welch (1983) criterion through the ‘coda’ package (Plummer et al., 2006).

The models tested were compared by deviance information criterion (DIC) as proposed by Spiegelhalter et al. (2002): D I C = D ( θ ¯ ) + 2 p D , where D ( θ ¯ ) is a point estimate of the deviance obtained by replacing the parameters with their a posterior means estimates in the likelihood function, and p_D is the effective number of parameters in the model. Models with smaller DIC should be preferred to models with higher DIC. However, differences (D) between DIC values of models a and b are given by D = | DIC_a – DIC_b |, and thus, if D < 5, there is no significant difference among the compared models; if 5 ≤ D ≤ 10, the difference is significant; and if D > 10, the difference is highly significant.

Correlations and principal component analysis (PCA) were performed using ‘BayesianFirstAid’ (Bååth, 2014) and ‘bPCA’ (Smyčka and Keil, 2020) packages, respectively. In both analyses, the median scores were reported with their respective 95% highest posterior density (HPD) intervals. Correlation estimates were considered significant when the HPD intervals did not overlap zero. The marginal a posteriori distributions were performed considering non-informative a priori distributions for all model parameters. A total of 100,000 values were generated by MCMC (Monte Carlo Markov Chain) process, assuming a burn-in period and thinning interval of 10,000 and 10 iterations, respectively. The MCMC convergence was verified using the Heidelberger and Welch (1983) criterion through the ‘coda’ package (Plummer et al., 2006).

Genetic progress was calculated for each trait evaluated in both N conditions using the following simple linear regression model:

Where: β ₀ is the intercept or linear coefficient, β ₁ is the slope or angular coefficient to adjust the data of variable y _i as a function of X _i (year of cultivar release), and ϵ_i are errors normally distributed with mean zero and variance σ², that is ϵ_i ~ N(0, σ²). The parameter σ² represents a variability in which the results differ from their predictions based on the model. This model can also be described as y _i ~ Normal (β ₀ + β ₁ X _i + σ²). The difference between genetic progresses (Δβ ₁) at high and low N fertilization was considered significant when there was no overlap between their HPD intervals.

The marginal a posteriori distributions were performed considering non-informative a priori distributions for all model parameters using the ‘MCMCglmm’ (Hadfield, 2010) package. A total of 100,000 values were generated by MCMC (Monte Carlo Markov Chain) process, assuming a burn-in period and thinning interval of 100,000 and 10 iterations, respectively. The MCMC convergence was verified using the Heidelberger and Welch (1983) criterion through the ‘coda’ package (Plummer et al., 2006).

Bayesian additive main effects and multiplicative interaction (BAMMI) method was used to study cultivar stability following the model below:

Where: 1 _n is the vector of the order n × 1, µ is the overall mean, X ₁ is the matrix of genotypes of order n × g, τ is the effect vector g × 1 for genotypes, X ₂ and δ are the matrices for environments of the order n × a and the effect vector a × 1 for environment, respectively. λ_k is the singular value for the k^th principal component, t is the number of multiplicative terms [t ≤ min (g, a) − 1), α_k and γ_k are the singular vectors of k for genotypes and environment, respectively; and ϵ is the vector n of error effect. Vector ϵ has a multivariate normal distribution with zero mean and variance–covariance matrix σ ϵ 2 I n Thus, the vector y also has a multivariate normal distribution.

For the BAMMI model, the estimation of the parameters of the above equation model assumes that the conditional distribution of y, given that µ, τ, δ, λ, α, γ, and σ ϵ 2 is a multivariate normal distribution:

Where: I_n is the identity matrix of order n. The a priori distributions used for the parameters were the same as those presented by Crossa et al. (2011). Subscripted symbols μ and σ ² denote mean and variance of the a priori distribution of whatever parameter is shown as subscript:

λ k ∣ μ λ k , σ λ k 2 ∼ N + ( μ λ k , I e σ λ k 2 ) with the restrictions λ_k > 0 and λ_{k– 1} ≥ λ_k ;

Where: N denotes the normal distribution, N⁺ is the positive normal distribution, and Inv − Scaled − χ ² is the inverse χ² distribution. In this study, the a priori distributions were non-informative. The value zero was used as prior distribution for the mean in all genotypic and environmental effects and high values for the variances, resulting in: µ_µ = 0, µ_τ = 1 _g × 0, µ_δ = 1a × 0 and µ_λk = 0, and for the variances σ μ 2 , σ τ 2 , σ δ 2 e σ λ k 2 = 1 × 10¹⁵.

The a posteriori distribution was estimated by:

A total of 1,000,000 values were generated by MCMC (Monte Carlo Markov Chain) process, assuming a burn-in period and thinning interval of 100,000 and five iterations, respectively. The MCMC convergence was verified using the Heidelberger and Welch (1983) criterion through the ‘coda’ package (Plummer et al., 2006). The BAMMI method was carried out using R script developed by Crossa et al. (2011).

DIC results showed a positive evidence of a triple interaction among genotype, environment, and N levels (G × E × N) for all evaluated traits since the full model presented lower DIC values ( Table 3 ). However, for trait HI, the DIC values difference among full, reduced 1, and reduced 2 models was less than five, indicating that there is no difference among all models tested.

The overall mean for all evaluated traits and their respective HPD intervals are shown in Table 4 . In general, there was a reduction of SY under low N conditions. However, there was an increase in NUsE and NUpE under low N when compared to high N conditions. There were no differences evident between high and low N conditions for the NUtE, Prot, and HI traits since there was an overlap of the HPD intervals in both N conditions. The performance of cultivars is presented in Supplementary Tables 2–4 .

Correlation coefficients (r) and their respective HPD intervals among traits evaluated are shown in Table 5 . In general, there was a moderate agreement between the correlations of high and low N conditions. SY, HI, and NUtE showed a positive correlation in both N levels. On the other hand, SY, HI, and NUtE showed a negative correlation with Prot, which correlated positively with NUsE (r = 0.47, 95% HPD = 0.36; 0.58), and NUpE (r = 0.54, 95% HPD = 0.46; 0.62) under low N. NUsE showed a positive correlation with NUpE under low (r = 0.98, 95% HPD = 0.97; 0.99) and high N (r = 0.87, 95% HPD = 0.84; 0.90) conditions. In addition, we observed a positive and negative correlations between NUsE and NUtE under high (r = 0.23, 95% HPD = 0.12; 0.34) and low N (r = −0.33, 95% HPD = −0.44; −0.22) levels, respectively.

PCA graphics are shown in Figure 1 . The first two principal components (Comp.1 and Comp.2) explained 76.1% (95% HPD = 60.1; 83.2%) and 14.3% (95% HPD = 6.1; 20.2%) of the total variation observed under high N condition. Under low N, Comp.1 and Comp.2 explained 74.2% (95% HPD = 57.4; 82.1%) and 15.2% (95% HPD = 8.0; 23.4%) of the total variation detected, respectively ( Supplementary Figure 2 ). NUsE, NUpE, and SY were the traits that most contributed to discriminate cultivars under high and low N since they showed the highest absolute values for their eigenvectors in Comp.1 and Comp.2 ( Supplementary Figure 3 ). The projection of NUtE and SY vectors showed the same direction in both N conditions. However, their vectors were projected in the opposite direction in relation to Prot, indicating a negative correlation of NUtE and SY to Prot. The IAC Alvorada cultivar was projected more distantly among all cultivars evaluated for both N conditions, probably because this cultivar showed the highest NUpE and NUsE means under high and low N conditions, respectively.

Genetic progress estimates for the traits evaluated under high and low N are shown in Table 6 . The significance of the linear regression model was observed only for SY in both N conditions since the estimates of the β ₁ parameter and their respective HPD intervals did not overlap the value of zero. Under high N, genetic progress estimate was 13.1 kg ha⁻¹ (95% HPD = 9.3; 17.0 kg ha⁻¹) or 0.53% year⁻¹ (95% HPD = 0.39; 0.65%). Under low N condition there was a progress of 10.2 kg ha⁻¹ (95% HPD = 6.1; 14.5 kg ha⁻¹) or 0.48% year⁻¹ (95% HPD = 0.31; 0.64%). There was no difference between the a posteriori slope estimates (Δβ ₁), indicating that the genetic progress of SY was the same in both N conditions.

The a posteriori means of genotypic (τ) and environmental (δ) effects and their respective HPD intervals are shown in Tables 7 and 8 , respectively. Overall mean (µ) of SY from 40 carioca common bean cultivars evaluated in eight contrasting N conditions was 2,155.71 kg ha⁻¹ (95% HPD = 2141.99; 2169.54 kg ha⁻¹). Twenty-one cultivars showed a posteriori positive means with HPD intervals that did not overlap the zero value ( Table 7 ). Among these cultivars, the ones that stood out the most were: IPR Sabiá (376.61 kg ha⁻¹; 95% HPD = 363.95; 389.91 kg ha⁻¹), IPR Quero-quero (319.24 kg ha⁻¹, 95% HPD = 306.23; 332.79 kg ha⁻¹), IPR Bem-te-vi (288.76 kg ha⁻¹; 95% HPD = 275.63; 302.05 kg ha⁻¹), IAC Sintonia (274.92 kg ha⁻¹; 95% HPD = 261.39; 287.95 kg ha⁻¹), IPR Campos Gerais (235.44 kg ha⁻¹; 95% HPD = 222.72; 24.06 kg ha⁻¹), and BRS Notável (214.96 kg ha⁻¹; 95% HPD = 201.86; 228.57 kg ha⁻¹). It is worthwhile to mention that not all environments under high N conditions were considered as favorable environments since they did not show a posteriori positive means ( Table 8 ).

Genotypic and environmental scores with their respective HPD intervals for the first two principal components (Comp.1 and Comp.2) are shown in Figure 2 . Comp.1 and Comp.2 explained 43.33 and 21.47% of G × E interaction, respectively. The HPD intervals overlapping in the central point indicate the presence of genotypic or environmental stability. In addition, HPD intervals overlapping among genotypes or environments indicate similar responses among them. The genotypic scores showed that 25 cultivars have high stability since their HPD intervals overlapped with zero values on both axes ( Figure 2A ). The environmental scores allowed us to observe that PG17–LN was the environment which contributed more to G × E interaction since HPD interval did not overlap zero on the axis of Comp.1 ( Figure 2B ).

Comp.1 and Comp.2 coordinates allowed us to infer about the genotype specific adaptation to certain environments ( Figure 2 ). Considering only genotypes and environments with significant contribution to G × E interaction, seven cultivars (Carioca, IAPAR 14, IAPAR 72, IAPAR 81, BRS Horizonte, IPR Eldorado, and TAA Bola Cheia) showed specific adaptation to environment PG17–LN. The cultivars that showed to be more adapted and stable were: IPR Sabiá, IPR Quero-quero, IPR Bem-te-vi, IAC Sintonia, and IPR Campos Gerais. These cultivars presented high SY mean values and high behavior predictability for different environmental conditions.

The environmental correlations of SY between high and low N conditions are shown in Figure 3 . In general, correlation estimates were positive and significant in all evaluated environments since their HPD intervals did not overlap the zero. Moderate correlation estimates were observed in PG18 (r = 0.50, 95% HPD = 0.37; 0.62), while weak correlations were observed in STO18 (r = 0.39, 95% HPD = 0.23; 0.53), LD2017 (r = 0.33, 95% HPD = 0.18; 0.46), and PG2017 (r = 0.29, 95% HPD = 0.15; 0.43).

The best adjustment of the full model for all agronomic traits evaluated indicated that cultivars have a differential behavior based on the combination of N fertilization levels and environments in which the cultivars were submitted. The interaction between genotypes and environments is frequently observed in common bean breeding programs in Brazil (Barili et al., 2016; Pereira et al., 2018) and worldwide (Bruno et al., 2017; Caproni et al., 2018).

SY increased under high N top-dressing (17.51%), corroborating a previous study that reported N as the macro-element that most influences SY in common bean crop (Barros et al., 2018). Similar results were also observed for 16 Brazilian carioca common bean cultivars which observed an increase of approximately 21.71% in plants under high N condition (Leal et al., 2019).

The majority of the NUsE and NUpE evaluated in carioca common bean cultivars in this study were higher under low N compared to the high N condition. Studies have already been reported that plants improve their NUsE under low N availability, especially because under high N conditions this element can be lost by leaching, denitrification, and volatilization processes (Shibata et al., 2017).

Several studies also reported the influence of environmental factors on Prot in common beans (Pereira et al., 2017; Fidelis et al., 2019). A study using 40 common bean genotypes (carioca and black types) observed 18.33 and 19.36% of Prot in common beans cultivated in low and high N condition, respectively, where 17 genotypes did not show an increase in Prot in the presence of N top-dressing (Fidelis et al., 2019). In the present study, although the overall mean of Prot was not significantly altered by N top-dressing, seven cultivars (BRS Requinte, BRS Horizonte, IAC Alvorada, FTS 65, IPR Campos Gerais, BRS Notável, and IPR Quero-quero) showed an increase in Prot under the high N condition.

We did not observe significant differences for HI values between both N conditions. However, some studies observed an increase of HI through N top-dressing fertilization (Fageria and Baligar, 2005; Fageria et al., 2014). HI is a measure of biomass partitioning, which indicates the fraction of total above-ground biomass allocated in the seeds and is considered to be one of the most determinant traits of common bean tolerance to several abiotic stresses (Assefa et al., 2015; Amanuel et al., 2018). In general, the genotypes that are more tolerant to abiotic stresses have a higher HI value, that is, a greater ability to partition photo-assimilates from vegetative organs to seeds (Assefa et al., 2015). In common bean, HI values have varied between 0.3 and 0.6, being considered a complex trait and strongly influenced by environment (Pinto Júnior et al., 2018). According to Beebe et al. (2008), the increase in HI is a key strategy in the improvement of beans, mainly under abiotic stress conditions.

NUtE was the only component of NUsE to present a positive correlation with SY in both N conditions, indicating that this trait is the main component of NUsE in common bean crop, especially when the breeding program aims to increase SY. Several metabolic and physiological mechanisms are associated with NUtE in plants, including greater photosynthetic efficiency per unit of N, improved partition of carbohydrates, storage, and N remobilization through senescent tissues (Han et al., 2015). NUtE increase is often associated to the stay-green character since genotypes with this characteristic maintain leaves photosynthetically active for a longer period (Zhang et al., 2020). HI has also been associated with NUtE since its increase contributes to a more efficient production of seed biomass by total biomass accumulated in plants (Han et al., 2015). In this study, we also observed a positive correlation between HI and NUtE in both N conditions.

We observed a negative correlation between NUtE and NUpE under low N condition, similar to what has been previously described in maize (Gallais and Hirel, 2004). These authors reported that NUtE is related to protein degradation in senescent leaves, and they assumed that the use of N occurs mainly when absorption is reduced or interrupted during stresses and/or natural senescence. In addition, the genetic variability of NUsE is a function of the existing variability of NUpE in high N condition, while in low availability it is a function of NUtE. Our results partially corroborated this study since we observed a genetic variability in NUpE with a greater contribution of NUsE in both N conditions.

The observed negative correlations between SY and Prot may be related to altered patterns of carbon and N metabolism, and this result is corroborated by different studies carried out in common bean (Razvi et al., 2018) and other crops, such as, oilseed rape (Stahl et al., 2017) and wheat (Thorwarth et al., 2019). Several studies investigated the genetic basis of this negative correlation and showed that pleiotropic effects, environmental conditions, and management techniques can influence the relationship between these variables (Rozbicki et al., 2015; Thorwarth et al., 2019). Another study reported genetic correlations between Prot and SY (r_g = 0.51; p < 0.01) similar to those obtained in the present study using 140 common bean recombinant inbred lines (RIL) under contrasting N condition (Farid et al., 2017).

Our genetic progress related to SY was similar under low (0.48%, 95% HPD = 0.31; 0.64%) and high N (0.53%, 95% HPD = 0.39; 0.65%), indicating that modern cultivars do not demand more N fertilization to be more productive than cultivars released earlier. Around ten years are needed to develop a new common bean cultivar (Moreira et al., 2010). As the number of selected lines is reduced, the range of environments in which they are tested is wider. Among all these trials, moderate N stresses surely occur unintentionally. Thus, the selection process may already mix high and low N environments explaining in part the similar genetic progress observed in the present study under both N conditions.

Genetic progress estimates for SY using Brazilian carioca and black common beans cultivars have already been reported (Chiorato et al., 2010; de Faria et al., 2013; Barili et al., 2016), but not under low N condition. For example, de Faria et al. (2013) assessing lines and cultivars representative of 22 years of the Embrapa breeding program reported progress of 0.72% year⁻¹ (17.3 kg ha⁻¹). In the IAC bean breeding program, genetic progress was 1.07% year^–1 (13.17 kg ha^–1) between 1989 and 1996. Differences between germplasm, environments, methodologies, and statistical methods are the main reasons for the differences observed among studies (de Faria et al., 2018).

We did not detect a significant difference in Prot content among modern and old cultivars. Our hypothesis is that selection on Prot may only result in the elimination of low Prot lines and not in increasing Prot, likewise reported in wheat (Cormier et al., 2013) and oilseed rape (Stahl et al., 2017) crops. Common bean breeding program objectives were clearly to increase SY, and the concern and interest in biofortification is recent. Since there is a negative correlation between SY and Prot, an alternative could be to improve protein composition to increase the levels of essential amino acids (mainly methionine, cysteine, and tryptophan) and/or decrease protein digestibility inhibitors (Rezende et al., 2017).

Although there is wide genetic variability observed for NUsE, NUtE, and NUpE in this study, no genetic progress was detected for these traits. Possibly the common bean breeding programs are not focusing in the improvement of NUsE since common beans are able to fix N₂ by symbiosis. However, NUsE genetic progress and their components have been reported in other crops that do not have this N₂ fixation ability, such as, wheat (Laidig et al., 2017), maize (Mueller et al., 2019), and oilseed rape (Stahl et al., 2017). In addition, the evaluation and selection for NUsE is considered an extremely expensive process and may not be implemented in breeding programs soon. High-throughput phenotyping methods are currently being developed (Neilson et al., 2015; Kefauver et al., 2017) and can become an important tool in the improvement of NUtE and/or NUpE in addition to molecular selection on genes or quantitative trait loci (QTL) (Cormier et al., 2013).

In our study, the genetic progress of HI was absent and indicated that superior genotypes selection was not performed focusing on this agronomic trait. However, the genetic progress of HI has been historically reported in several other crops (Cormier et al., 2013; Mueller et al., 2019; Todeschini et al., 2019). In a meta-analysis including eleven different crops, the authors reported a positive relationship between SY and HI in cereal (maize, wheat, barley, and oat), oilseed (soybean, canola, and flax), and pulse crops (pea, chickpea, and lentil), with the exception of potato crops only (Fan et al., 2017). The same authors affirmed that this positive linear relationship indicates that SY improvement was achieved in part through HI, and that plant breeders should focus greater attention on this agronomic trait to develop new cultivars.

The development of genotypes with environment adaptive plasticity, good stability, and high SY mean is one way to mitigate G × E interaction effects (Resende et al., 2019). Statistical methods to study adaptability and stability have been developed and widely used by plant breeding programs (Nascimento et al., 2020). The BAMMI method stands out for its power to explain G × E interaction compared to other methods (Teodoro et al., 2019). Another advantage of the BAMMI model is the presence of HPD intervals for genotypic and environment scores (Oliveira et al., 2015). These credibility intervals lead to greater precision to infer genotypic and environment stability by eliminating subjective mean scores in relation to the proximity to the central point of the biplot (coordinates 0 and 0).

In the present study, the large number of cultivars considered stable is due to the selection for stability in several environment conditions by their respective breeding programs. In addition, there was no relationship between stability and release year since both modern and old cultivars showed high stability. The modern cultivars IPR Sabiá, IPR Quero-quero, IPR Bem-te-vi, IAC Sintonia, and IPR Campos Gerais were the most adapted and stable in the present study. These results confirm that common bean breeding programs are developing cultivars that combine SY and stability for different environmental conditions.

We observed that not all environments under high N condition were considered favorable environments by BAMMI. In addition, the positive correlations between high and low N for all environments explain the similar genetic progress among them for SY. These results indicated that even though common bean breeding programs have made selection under N fertilization, modern cultivars do not require high N levels to achieve their genetic potential. Similar results were already reported in maize (Haegele et al., 2013) and wheat (Cormier et al., 2013).