Genetic Progress of Seed Yield and Nitrogen Use Efficiency of Brazilian carioca Common Bean Cultivars Using Bayesian Approaches

Common bean (Phaseolus vulgaris L.) is one of the most important crops worldwide and is considered an essential source of proteins, fibers, and minerals in the daily diet of several countries. Nitrogen (N) is considered the most important nutrient for common bean crop. On the other hand, the reduction of chemical fertilizers is a global challenge, and the development of cultivars with more N use efficiency (NUsE) is considered one of the main strategies to reduce the amount of N fertilizers. Genetic progress of NUsE has been reported in several crops; however, there was still no quantity in common bean. In this study, our goal was to analyze the genetic progress of seed yield (SY) and NUsE-related traits of 40 carioca common bean cultivars release from 1970 to 2017 in eight environments under low (zero) or high N (40 kg ha−1) in top-dressing. Genetic progress, principal component analysis, correlations among traits, and cultivar stability were analyzed using Bayesian approaches. The lowest values of the deviance information criterion (DIC) for the full model tested indicated the presence of the genotype × N × environment interaction for all evaluated traits. Nitrogen utilization efficiency (NUtE) and nitrogen uptake efficiency (NUpE) were the traits that most contributed to discriminate cultivars. The genetic progress of SY under high N (0.53% year−1, 95% HPD = 0.39; 0.65% year−1) was similar to that obtained in low N conditions (0.48% year−1, 95% HPD = 0.31; 0.64% year−1). These results indicate that modern cultivars do not demand more N fertilizers to be more productive. In addition, we observed a high genetic variability for NUsE-related traits, but there was no genetic progress for these variables. SY showed negative correlation with seed protein content (Prot) in both N conditions, and there was no reduction in Prot in modern cultivars. Both modern and old cultivars showed adaptability and stability under contrasting N conditions. Our study contributed to improve our knowledge about the genetic progress of common bean breeding program in Brazil in the last 47 years, and our data will help researchers to face the challenge of increase NUsE and Prot in the next few years.


INTRODUCTION
Common bean (Phaseolus vulgaris L.) is considered the main legume species for human consumption and represents an essential source of proteins, carbohydrates, fibers, and trace minerals in several countries worldwide (Myers and Kmiecik, 2017). Brazil is one of the main world producers and consumers of common beans. In 2019, it produced about 2,360 tons of seeds in an area of 1,619 ha (CONAB, 2020a). Common bean is a basic food used daily in the Brazilian diet with a consumption of 17 kg person −1 year −1 (Ribeiro et al., 2020). To meet this demand, common beans are grown throughout the all year under different cultivation systems and with different technological levels (Heinemann et al., 2016). In general, Brazilian consumers prefer the Mesoamerican carioca (cream seed coat with brown stripes) and black beans, which represent around 70 and 15% of total common bean production, respectively (Pereira et al., 2019).
Although Brazilian common bean production has grown in recent decades, the cultivated area has decreased considerably since 1990 (Supplementary Figure 1). This increase in Brazilian production is mainly due to the use of chemical fertilizers, the adoption of new production technologies, and the release of more productive cultivars (Tsutsumi et al., 2015). According to the Ministry of Agriculture, Livestock, and Food Supply (MAPA, 2020), more than 350 bean cultivars are registered in Brazil, developed mainly by public research institutes, such as the Brazilian Agricultural Research Corporation (Embrapa) (Lopes et al., 2012), Instituto de Desenvolvimento Rural do Parana-IAPAR-EMATER (IDR-Parana) (Moda-Cirino et al., 2012), and Instituto Agronômico de Campinas (IAC) (Carbonell et al., 2012).
Common bean plants can fix atmospheric nitrogen (N 2 ) interacting symbiotically with Rhizobium genus bacteria (Appelbaum, 2018). This symbiotic process (biological nitrogen fixation-BNF) is responsible for around 40% and chemical fertilizer for 60% of the total nitrogen (N) demand for this crop (Argaw et al., 2015;Dwivedi et al., 2015). N is the most important macronutrient influencing common bean yield (Barros et al., 2018). Therefore, the development of cultivars that have greater N use efficiency (NUsE) is considered one of the main strategies to reduce the amount of chemical fertilizers in sustainable farming systems (Mueller et al., 2019).
NUsE can be divided into two components (Moll et al., 1982): (i) N utilization efficiency (NUtE) which is the capacity of the plant to convert N into seed biomass; and (ii) N uptake efficiency (NUpE) which means the capacity of the plant to absorb N from the soil. In legume crops, NUpE is based on the N acquisition from the soil and efficiency of the BNF process (Diacono et al., 2019). Several studies have already described the importance of BNF for common beans; however, studies related to NUsE components are still scarce (Farid et al., 2017;Kamfwa et al., 2019).
Information on the quantification of genetic progress for any crop is important for plant breeders, allowing them to design future breeding program strategies (Stahl et al., 2017;de Faria et al., 2018). In Brazil, genetic progress studies related to common bean seed yield have been extensively reported (Chiorato et al., 2010;de Faria et al., 2013;Barili et al., 2016), but there is no study combining information on the genetic progress of seed yield and NUsE under N contrasting conditions. In addition, a negative correlation between seed yield and protein content (Prot) may have caused its decrease over the years as reported in wheat (Laidig et al., 2017). In the face of this challenge, the question arises whether or not breeding is directing NUsE and Prot improvement.
Bayesian methods have been increasingly used in plant breeding studies (da Silva et al., 2020;Nascimento et al., 2020). The fundamental idea of the Bayesian approach is to describe all uncertainty and variation using probability distributions. Bayesian analysis offers the ability to utilize data containing unbalanced structures, select and study complex models, obtain more precise credibility intervals, and incorporate a priori information (Aczel et al., 2020). In this way, our main goal was to quantify the genetic progress of seed yield and NUsE-related traits in Brazilian carioca common bean cultivars using Bayesian approaches. Our hypothesis was that modern cultivars are more responsive and have improved NUsE when compared to cultivars released in the past since breeding programs are evaluating and selecting superior genotypes under N top-dressing.

Plant Material
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-Parana), Londrina, Brazil, and subsequently multiplied in order to standardize seed germination. Information about the genealogy of the cultivars is available in Supplementary Table 1.

Experimental Design
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 −1 (urea form) at V 3 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 −1 of NPK (04-30-10).

Agronomic Trait Analysis
Five uniform and representative plants were collected at physiological maturation (R 9 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 −1 with moisture of 13%) was obtained after manual removal and mechanical threshing of plants from the two  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.

Fitting Models
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, s²) 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) 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
Correlations and principal component analysis (PCA) were performed using 'BayesianFirstAid' (Bååth, 2014) and 'bPCA' (Smycǩa 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
Genetic progress was calculated for each trait evaluated in both N conditions using the following simple linear regression model: Where: b 0 is the intercept or linear coefficient, b 1 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 s², that is ϵ i~N (0, s²). The parameter s² represents a variability in which the results differ from their predictions based on the model. This model can also be described as y i~N ormal (b 0 + b 1 X i + s²). The difference between genetic progresses (Db 1 ) 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).

Adaptability and Stability Analyses
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 1 is the matrix of genotypes of order n × g, t is the effect vector g × 1 for genotypes, X 2 and d are the matrices for environments of the order n × a and the effect vector a × 1 for environment, respectively. l k is the singular value for the k th principal component, t is the number of multiplicative terms [t ≤ min (g, a) − 1), a k and g 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 s 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 µ, t, d, l, a, g, and s 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 m and s 2 denote mean and variance of the a priori distribution of whatever parameter is shown as subscript: s k~s pherical uniform distribution on the corrected subspace; g k~s pherical uniform distribution on the corrected subspace; and Where: N denotes the normal distribution, N + is the positive normal distribution, and Inv − Scaled − c 2 is the inverse c 2 distribution. In this study, the a priori distributions were noninformative. 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, µ t = 1 g × 0, µ d = 1a × 0 and µ lk = 0, and for the variances s 2 m ,s 2 t ,s 2 d e s 2 l k = 1 × 10 15 . 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).

Fitted Models and Influence of Nitrogen Fertilization
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 and Principal Component Analysis
Correlation coefficients (r) and their respective HPD intervals among traits evaluated are shown in Table 5 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 b 1 parameter and their respective HPD intervals did not overlap the value of zero. Under high N, genetic progress estimate was 13.1 kg ha −1 (95% HPD = 9.3; 17.0 kg ha −1 ) or 0.53% year −1 (95% HPD = 0.39; 0.65%). Under low N condition there was a progress of 10.2 kg ha −1 (95% HPD = 6.1; 14.5 kg ha −1 ) or 0.48% year −1 (95% HPD = 0.31; 0.64%). There was no difference between the a posteriori slope estimates (Db 1 ), indicating that the genetic progress of SY was the same in both N conditions.

Adaptability and Stability
The a posteriori means of genotypic (t) and environmental (d) 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 −1 (95% HPD = 2141.99; 2169.54 kg ha −1 ). 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á ( . 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 Sabia, IPR Quero-quero, IPR Bemte-vi, IAC Sintonia, and IPR Campos Gerais. These cultivars presented high SY mean values and high behavior predictability for different environmental conditions.  Table 1. 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).

DISCUSSION
Our study was the first to quantify the genetic progress of carioca common beans in Brazil over the past 47 years under N  Angular coefficient (b 1 ) of simple linear regression model were considered significant (*) when HPD intervals did not overlap the value of zero or not significant (ns) when there was overlap.
2 Genetic progress differences (Db 1 ) between high and low N conditions were considered significant (*) when there was no overlap of their HPD intervals or not significant (ns) when there was overlap. 7 | A posteriori mean of Bayesian additive main effects and multiplicative interaction (BAMMI) for seed yield (kg ha −1 ) of 40 carioca common bean cultivars in eight environments for the overall mean (µ) and genotypic effects (t i ) with their respective 95% highest posterior density (HPD).

Parameter
Mean 95% HPD Adaptability contrasting conditions. In Brazil, these studies have mainly focused on the impact of breeding programs on plant architecture and yield components. For example, Barili et al. (2016) reported a change in modern cultivar architecture since there was a substitution of prostrate cultivars by erect growth habit. This change allowed mechanical harvesting with reduced losses, lower incidence of disease, and improved grain quality. In addition, the increase in common bean production has also been associated with genetic progress in the number of pods per plant (5.6% year −1 ), seeds per plant (4.5% year −1 ), mass of seeds (2.08% year −1 ), and tolerance to lodging (2.0% year −1 ) (de Faria et al., 2013;Barili et al., 2016). In common bean, there are still no studies quantifying the genetic progress for NUsE. On the other hand, an increase in NUsE has been reported in several other crops (Laidig et al., 2017;Stahl et al., 2017;Mueller et al., 2019).

Interaction Between Factors and Influence of Nitrogen Fertilization
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 macroelement 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 Notavel, 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 Juńior 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.

Agronomic Traits Correlations
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).

Genetic Progress
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 −1 (17.3 kg ha −1 ). 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 2 by symbiosis. However, NUsE genetic progress and their components have been reported in other crops that do not have this N 2 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.

Adaptability and Stability
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 Sabia, 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).

CONCLUSION
Our study is the first to quantify the genetic progress for SY and NUsE-related traits in Brazilian carioca common bean cultivars launched from 1970 to 2017 under N contrasting conditions. Among the traits evaluated, there was genetic progress only for SY under high and low N in top-dressing. The similar genetic progress in both N conditions rejected our hypothesis that modern cultivars are more N-dependent to reach their productive potential. Our results indicate that the challenges for common bean breeders for the coming years will be to improve the NUsE of new cultivars in order to reduce the dependence on N fertilizers as well as to increase the levels and/or quality of Prot in the seeds.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUITIONS
DZ, IM, VM-C, and LG conceived and designed the experiments. DZ, IM, and JD conducted the experiments and collected the data. DZ and GF performed the statistical analyses. DZ wrote the original draft. S-IS, VM-C, JN, GF, CS, and LG read and edited the manuscript. All authors contributed to the article and approved the submitted version.