Genome-Wide Association Study Identifies Genomic Regions for Important Morpho-Agronomic Traits in Mesoamerican Common Bean

The population growth trend in recent decades has resulted in continuing efforts to guarantee food security in which leguminous plants, such as the common bean (Phaseolus vulgaris L.), play a particularly important role as they are relatively cheap and have high nutritional value. To meet this demand for food, the main target for genetic improvement programs is to increase productivity, which is a complex quantitative trait influenced by many component traits. This research aims to identify Quantitative Trait Nucleotides (QTNs) associated with productivity and its components using multi-locus genome-wide association studies. Ten morpho-agronomic traits [plant height (PH), first pod insertion height (FPIH), number of nodules (NN), pod length (PL), total number of pods per plant (NPP), number of locules per pod (LP), number of seeds per pod (SP), total seed weight per plant (TSW), 100-seed weight (W100), and grain yield (YLD)] were evaluated in four environments for 178 Mesoamerican common bean domesticated accessions belonging to the Brazilian Diversity Panel. In order to identify stable QTNs, only those identified by multiple methods (mrMLM, FASTmrMLM, pLARmEB, and ISIS EM-BLASSO) or in multiple environments were selected. Among the identified QTNs, 64 were detected at least thrice by different methods or in different environments, and 39 showed significant phenotypic differences between their corresponding alleles. The alleles that positively increased the corresponding traits, except PH (for which lower values are desired), were considered favorable alleles. The most influenced trait by the accumulation of favorable alleles was PH, showing a 51.7% reduction, while NN, TSW, YLD, FPIH, and NPP increased between 18 and 34%. Identifying QTNs in several environments (four environments and overall adjusted mean) and by multiple methods reinforces the reliability of the associations obtained and the importance of conducting these studies in multiple environments. Using these QTNs through molecular techniques for genetic improvement, such as marker-assisted selection or genomic selection, can be a strategy to increase common bean production.

The population growth trend in recent decades has resulted in continuing efforts to guarantee food security in which leguminous plants, such as the common bean (Phaseolus vulgaris L.), play a particularly important role as they are relatively cheap and have high nutritional value. To meet this demand for food, the main target for genetic improvement programs is to increase productivity, which is a complex quantitative trait influenced by many component traits. This research aims to identify Quantitative Trait Nucleotides (QTNs) associated with productivity and its components using multi-locus genome-wide association studies. Ten morpho-agronomic traits [plant height (PH), first pod insertion height (FPIH), number of nodules (NN), pod length (PL), total number of pods per plant (NPP), number of locules per pod (LP), number of seeds per pod (SP), total seed weight per plant (TSW), 100-seed weight (W100), and grain yield (YLD)] were evaluated in four environments for 178 Mesoamerican common bean domesticated accessions belonging to the Brazilian Diversity Panel. In order to identify stable QTNs, only those identified by multiple methods (mrMLM, FASTmrMLM, pLARmEB, and ISIS EM-BLASSO) or in multiple environments were selected. Among the identified QTNs, 64 were detected at least thrice by different methods or in different environments, and 39 showed significant phenotypic differences between their corresponding alleles. The alleles that positively increased the corresponding traits, except PH (for which lower values are desired), were considered favorable alleles. The most influenced trait by the accumulation of favorable alleles was PH, showing a 51.7% reduction, while NN, TSW, YLD, FPIH, and NPP increased between 18 and 34%. Identifying QTNs in several environments (four environments and overall adjusted mean) and by multiple methods reinforces the reliability of the associations obtained and the importance of conducting these studies in multiple environments. Using these QTNs through molecular techniques for genetic improvement, such as marker-assisted selection or genomic selection, can be a strategy to increase common bean production.

INTRODUCTION
More than 24 million tons of common beans (Phaseolus vulgaris L.) are produced per year worldwide, and the main producing countries are located in Asia and the Americas (Rawal and Navarro, 2019). This crop is mainly grown by small producers, often in low fertility areas with low-level technology, resulting in low mean productivity (Broughton et al., 2003).
Increasing productivity is one of the main objectives of breeding programs. In this context, understanding the genetic constitution related to the productivity and of the production components are the basis for improvement (Kamfwa et al., 2015). In the cultivation of common beans, productivity is related to several morphological, agronomic, and physiological characteristics. The number of pods per plant (NPP), number of seeds per pod (SP), and seed weight are the primary components related to productivity, but other characteristics are also influential, such as growth rate, the capacity of the seeds to absorb photosynthates and plant architecture (Beebe et al., 2013;Resende et al., 2018;Assefa et al., 2019). Physiological and morphological characteristics such as days to flowering and maturity and resistance to pod shattering also significantly impact adaptability, biomass, and productivity (Zhang et al., 2015;Parker et al., 2020a).
Productivity and its components are quantitative traits and are highly influenced by the environment. Thus, understanding the relationship between these traits is very important for directing strategies and efforts in genetic improvement programs (Singh et al., 1991;Assefa et al., 2019). Traditional selection methods in plant breeding require intensive phenotyping fieldwork, with evaluations in several environments and years, resulting in high cost and a time-consuming process (Ikram et al., 2020). The use of molecular markers can increase efficiency and reduce the costs of phenotyping in plant breeding programs. Using molecular tools makes it possible to identify genomic regions related to the productivity and its components, which can be used in marker-assisted selection (Kamfwa et al., 2015). Genome-Wide Association Studies (GWAS) represent a powerful option for the genetic characterization of quantitative traits and have been widely used to analyze agronomic characteristics in plants (Wang et al., 2012;Contreras-Soto et al., 2017;Zhang K. et al., 2017;Cui et al., 2018;Ward et al., 2019;Ikram et al., 2020).
GWAS is a powerful tool enabling the study of different regions of the genome simultaneously using high-resolution mapping. Multiple polymorphisms occurring naturally within a species are identified in germplasm collections with limited genetic structure; genotypes that have traits of interest for breeding programs can be used preferentially to accelerate the application of GWAS studies (Korte and Ashley, 2013). Important genetic factors are identified based on the presence of linkage disequilibrium (LD), while taking into account the potential confounding effect of the population's structure (Kwak and Gepts, 2009) and kinship relations. Large and highly diverse association panels have unique recombination histories, allowing the detection of small and large genetic effects associated with a particular trait (Oladzad et al., 2019).
Although the statistical power in the detection of Quantitative Trait Nucleotides (QTNs) improves after controlling for the polygenic background of the experimental population under study, most of the small effects associated with complex traits are still not captured by the GWAS single-locus methods (Cui et al., 2018). Single-locus methods perform a one-dimensional scan of the genome; that is, they test one marker at a time, using several rigorous significance test corrections for multiple tests, such as the Bonferroni and False Discovery Rate (FDR) tests. However, these methods can be very conservative in eliminating true QTNs (He et al., 2019b). Multi-locus models are being developed to solve this problem. These models involve a multi-dimensional scanning of the genome, in which the effects of all markers are simultaneously estimated (Cui et al., 2018). The advantage of these models is that it is unnecessary to perform multiple test corrections; therefore, more markers associated with traits of interest are identified (Li et al., 2018).
Common beans of Mesoamerican origin are the most consumed in Brazil, with a preference for the Carioca and Black commercial classes (Gepts et al., 1988;Burle et al., 2010Burle et al., , 2011Delfini et al., 2017). Few GWAS have been directed toward diversity panels of plants of Mesoamerican origin and plants adapted to the Brazilian climatic conditions. Studies of genetic variation in accessions adapted to the target habitats are a powerful and effective approach to investigate the genetic architecture of complex traits, and later these natural allelic variations can be directly employed in breeding programs (Nakano and Kobayashi, 2020). Moreover, multi-locus methods were little explored in the cultivation of common beans. In this context, the present study's objective was to identify genomic regions related to morpho-agronomic traits in Mesoamerican common beans belonging to the Brazilian Diversity Panel (BDP) using the GWAS multi-locus methods.

Genetic Material, Field Experiments, and Phenotyping
In all, 178 Mesoamerican common bean accessions belonging to the BDP were evaluated (Delfini et al., 2021a). This panel consists of accessions that represent a large part of the variability present in Brazil and are adapted to tropical growing conditions. The phenotyping was conducted at the research stations of the Instituto de Desenvolvimento Rural do Paraná (IDR-Paraná), located in the state of Paraná, Brazil. The experiments were conducted in two seasons: the 2018 rainy season in the cities of Londrina (LDA_A18), Ponta Grossa (PG_A18), and Guarapuava (GUA_18); and the 2018/2019 dry season in Ponta Grossa (PG_S19), totaling four environments. The experiment used an incomplete block design with replications in sets. Five sets with two repetitions were used, and each set covered 50 entries, i.e., 46 accessions and four checks. Each plot consisted of four 2 m long rows, spaced 0.50 m between rows and with a density of 12 plants per linear meter. Management and treatments were conducted according to the technical recommendations for the plant's cultivation.
Seven plants from the two lateral lines of the plot were used to evaluate traits such as plant height (PH, in cm), first pod insertion height (FPIH, in cm), number of nodules (NN, count variable), pod length (PL, in cm), total number of pods per plant (NPP, count variable), number of locules per pod (LP, count variable), number of seeds per pod (SP, count variable), total seed weight per plant (TSW, in g), and 100-seed weight (W100, in g). Grain yield (YLD, in kg ha −1 and 13% moisture) was estimated by harvesting the two central lines.

Statistical Analysis of Phenotypic Data
An analysis of variance (ANOVA) was conducted using the PROC GLM function in the SAS software (SAS Institute, 2000). The following mathematical model was used: where is the general mean, Ai is the fixed effect of the i-th environment; S j is the effect of the j-th set; AS ij is the effect of the interaction between environments and sets; R/AS kij is the effect of the k-th repetition within the interaction between the i-th environment and the j-th set; G/S lj is the random effect of the l-th genotype within the j-th set; AG/S mlj is the effect of the interaction of environments and accessions within the j-th set, and e ijklm is the experimental error (Hallauer and Miranda Filho, 1988).
The means adjusted for each accession in each of the environments as well as the overall mean of all environments were obtained through the LSmeans option of the GLM procedure. The heritability (h 2 ) was estimated by the equation: h 2 = σ 2 G / σ 2 P , where genotypic (σ 2 G ) and phenotypic (σ 2 P ) variances were estimated by the following equations: σ 2 G = (QM G − QM E ) and σ 2 P = (QM G / ra) where QMG is the mean square of genotype within sets; QM E is the mean square of error, r is the number of replications, and a is the number of environments. The descriptive analysis was determined by the means adjusted for the two replications of each traits in each environment using the PROC UNIVARIATE function in the SAS software. Pearson's simple linear correlations were calculated and graphically presented using the R software 1 using the "corrplot" package (Wei and Simko, 2017).

Genotyping and Genome Wide Association Study
The genotyping-by-sequencing (GBS) technique was used to obtain the SNPs. The methodology used, as well as the results 1 https://www.r-project.org/ of the population structure and linkage disequilibrium (LD) analyses, are detailed in a previous work (Delfini et al., 2021a). In summary, GBS was conducted using the restriction enzyme CviAII (Ariani et al., 2016(Ariani et al., , 2017 and the data were imputed using Beagle software version 5 (Browning and Browning, 2016). After quality control using VCFtools version 0.1.15 (Danecek et al., 2011), 25,011 SNPs (MAF > 0.05) were used to perform GWAS analyses.
For conducting GWAS, mixed multi-locus models were used with the mrMLM.GUI software version 4.0 0 (Ya- Wen et al., 2019). Four different methods were used: mrMLM (Wang et al., 2016), FASTmrMLM (Tamba and Zhang, 2018), pLARmEB (Zhang J. et al., 2017), and ISIS EM-BLASSO . The critical values for significant associations were LOD ≥ 3 for all methods. Population structure and the kinship matrix were included in these models to minimize the identification of false positive associations and increase the statistical power of the analyses. The result of K = 2 was obtained by the Structure v2.3.4 software (Pritchard et al., 2000) [100,000 burn-in, 100,000 MCMC, and ten repetitions for hypothetical numbers of subpopulations (K) between 1 and 10], while the kinship matrix was obtained using the mrMLM.GUI software version 4.0.
The phenotypic values used were the adjusted means for each of the four environments and the overall adjusted mean (LDA_18, PG_18, GUA_18, PG_19, and LSmeans). In order to obtain more accurate results, only QTNs that presented repeatability, that is, detected at least three times by different methods or environments, were considered truly significant and used in the search for favorable alleles and candidate genes.

Favorable Alleles and Search for Candidate Genes
For each QTN, all accessions were divided into two groups based on the QTN genotype, that is, according to presence or absence of the favorable alleles. A t-test was then conducted to test if there was a significant difference in phenotypic mean between the two groups. Only the statistically stable QTNs between the environments, i.e., those that showed significant difference (P ≤ 0.05) in the phenotypes in at least three of the five environments (LDA_18, PG_18, GUA_18, PG_19, and LSmeans), were used as favorable alleles. The favorable genotype of each QTN was then selected, i.e., the genotype that causes the desired effect according to each trait, and these effects can be positive or negative in the case of PH. Then, the total number of favorable alleles for each trait was accounted for each accession, and, using a boxplot, visualized if the accumulation of these favorable alleles resulted in phenotypes with more desirable traits.
The identification of candidate genes was conducted at a physical distance of 296 kbp above and below the SNP associated with each of the assessed trait. This distance is the point at which the half decay of the LD, calculated with correction by population structure and relatedness (r 2 vs ), occurred (Delfini et al., 2021b). The genes present in the association region with known putative functions according to the GeneOntology (GO) 2 were identified based on the reference genome annotation of Phaseolus vulgaris v.2 published on the Phytozome v10.3 website. 3

Analysis of Variance, Heritability, and Environmental Effect
The analysis of variance showed a significant effect (P ≤ 0.01) of accessions and environments for all traits evaluated ( Table 1). Significant effects were also observed (P ≤ 0.01) for the Genotype x Environment (GE) interactions involving the traits PH, LP, SP, W100, and YLD. The coefficients of variation (CV) varied 3 http://phytozome.net between 5% (PL) and 26% (TSW). As for heritability estimates (h 2 ), the TSW, NN, NPP, and FPIH traits presented moderate values, between 0.54 and 0.68, while high values were detected for the other traits, YLD, SP, and LP had h 2 -values between 0.71 and 0.77, and PH, PL, and W100 had the highest values, 0.88, 0.94, and 0.94, respectively.
SP (4.49-7.55), and W100 (14.69-33.48) did not show significant variations in the minimum, maximum, and mean values for each environment.

Correlation Between Traits
Significant and positive correlations (P ≤ 0.05) were observed between PH, FPIH, and NN traits. The PL, LP, and SP traits also correlated positively with each other (Figure 2). Positive correlations were also observed between YLD and the primary components TSW and W100 (r = 0.31 and 0.32, respectively). In addition, YLD also correlated positively with LP (r = 0.24) and SP (r = 0.27), while NPP correlated positively with TSW (r = 0.65).
The highest number of QTNs was identified in the LSmeans dataset, followed by PG_18, GUA, 18, LDA_18, and PG_19. Among the GWAS multi-locus methods, the ISIS-EM-BLASSO method detected the highest number of SNPs, followed by pLARmEB, mrMLM, and FASTmrMLM. Considering only the 64 reliable QTNs, the environment ranking remained the same, being LSmeans the environment that detected the highest number of QTNs. For the methods, the number of stable QTNs detected was similar among the different methodologies, varying between 53 and 61. Looking at the efficiency of these methods, the number of QTNs considered reliable in relation to the initial number, the FASTmrMLM method stood out from the others (56%), followed by mrMLM (47%), pLARmEB (40%), and ISIS-EM-BLASO (34%).

Identification of Favorable Allelic Variations and Candidate Genes
Among the 64 QTNs considered to be reliable, 39 presented significant results for the t-test in at least three environments and were considered stable ( Table 2). These stable QTNs were used to identify alleles that were considered to be favorable for the traits PH (n = 10), NN (n = 6), PL (n = 6), NPP (n = 2), LP (n = 4), SP (n = 2), TSW (n = 3), W100 (n = 2), and YLD (n = 4). For FPIH, stable QTNs following the established criteria were not observed, thus, only for this trait, three QTNs that presented significant results through the t-test for the overall mean (LSmeans) were used, resulting in a total of 42 stable QTNs (Figure 3).
All QTNs that had a positive effect (increased values) were considered as favorable alleles. The only exception was PH, for which common bean breeding programs search for plants with shorter size, for the purpose of mechanizing the harvest. The accumulation of favorable alleles in the same accession resulted in a gradual increase, or decrease in the case of PH, in all traits (Figure 3). Comparing the mean values between genotypes with zero and those with the maximum number of favorable alleles, PH was the most influenced trait, revealing a difference of 51.7% between the two genotype groups, reducing the values from 92.8 to 48 cm. For those traits where higher values are favorable, the greatest increase was observed for NN (34%, from 9.47 to 12.7), followed by TSW (30%, from 17.4 to 22.7 g), YLD (27%, from 2,688 to 3,419 kg.ha −1 ), W100 (24%, from 19.7 to 24.5 g), FPIH (19%, from 16 to 19 cm), NPP (18%, from 18.4 to 21.7), PL (14%, from 9.3 to 10.6 cm), LP (10%, from 6.5 to 7.14), and SP (10%, from 6.1 to 6.7).
The distance determined using the average LD decay (296 kb) was used to select potential candidate genes at a specific QTN distance. Since the search was conducted in a large genomic region around the QTNs, many genes were identified for the 10 traits evaluated in this study, resulting in 1,528 genes with known putative functions, and of these, 74% were identified more than once in regions of overlap between QTNs. According to GO annotation, the genes were grouped in three functional categories: 55% had a molecular function, 32% had functions related to biological processes, and 13% were cellular components. In the molecular function category, the main functions detected were protein binding, ATP binding, and protein kinase activity; for biological processes, the functions were related to protein phosphorylation, oxidationreduction processes, and transcription regulation, and for cellular components, the functions were related to the membrane and integral components of the membrane.  , not significant at 5% probability level. PH, plant height (cm); FPIH, first pod insertion height (cm); NN, number of nodules; PL, pod length (cm); NPP, total number of pods per plant; LP, number of locules per pod; SP, number of seeds per pod; TSW, total seed weight per plant (gm); W100, 100-seed weight (gm); YLD, grain yield (kg.ha −1 ).

DISCUSSION
Although several studies already identified QTNs associated with morpho-agronomic traits in common beans using GWAS (Nemli et al., 2014;Kamfwa et al., 2015;Moghaddam et al., 2016;Soltani et al., 2016;Nascimento et al., 2018;Resende et al., 2018;Lei et al., 2020;Wu et al., 2020), panels composed exclusively of common beans of Mesoamerican origin and adapted to environmental conditions in Brazil had not been sufficiently explored (Valdisser et al., 2020). Moreover, few studies using the GWAS multilocus approach have been conducted in common bean. The use of the GWAS multi-locus methods has grown in recent years, becoming one of the main tools to identify molecular markers associated with traits of interest, especially for traits considered complex, i.e., controlled by multiple genes of small effect and highly influenced by the environment (Wen et al., 2018;Zhang et al., 2019;Yang et al., 2020).
A GxE (Genotype by Environment) interaction was observed for most of the morpho-agronomic traits evaluated in the present study, indicating that the accessions' differential behavior depends on the evaluation environments. The presence of GxE interaction is frequently observed in GWAS studies, where it interferes with the occurrence of the QTN x Environment interaction (Pan et al., 2018;Fattahi and Fakheri, 2019). The estimates of h 2 obtained in the present study were similar to those observed in literature (Rana et al., 2015;Asfaw et al., 2017). The h 2 is the central parameter of any breeding program, used to estimate the response to selection and explain the  proportion of phenotypic variation due to genetic variations (Falconer and Mackay, 1996). Correlations were observed among the traits related to the pods (PL, LP, and SP), the plant architecture traits (PH, FPIH, and NN), and the production components (TSW, W100, and YLD). Similar results were reported in several studies on common beans (Rana et al., 2015;Asfaw et al., 2017;Nadeem et al., 2020). The positive correlations observed between YLD and LP, SP, TSW, and W100 corroborate the possibility of indirect selection of YLD through these traits. However, Asfaw et al. (2017), studying the relationship between morpho-agronomic traits in 202 accessions of Andean and Mesoamerican origin, reported that the correlation between YLD and W100 occurs only in Andean common beans. Furthermore, the same authors recommended the indirect selection of YLD through NPP, independent of the gene pool. Although no positive correlation between YLD and NPP was found in this study, moderate correlations with TSW and SP were observed.
Quantitative genetics assumes that the genetic correlations between traits can be attributed to gene linkage and/or pleiotropy (Saltz et al., 2017). If pleiotropy is the main reason for genetic correlations between two traits, the same QTN can be identified in both traits. However, if gene linkage is the main reason, an overlap of the location between QTNs is expected. Thus, the pleiotropic QTNs identified between the PL-LP and NPP-TSW traits may be considered one of the causes of the correlations observed in these traits. On the other hand, the high number of QTNs identified in overlapping genomic regions indicates that genetic linkage may be the leading cause of the observed genetic correlations in the other assessed traits. Pleiotropy or linkage effects have been reported for many traits such as productivity, biomass, and plant height (Soltani et al., 2016).
Among the ML-GWAS methods used, the ISIS-EM-BLASSO detected the highest number of SNPs. However, it was the least efficient based on the number of verified QTNs. On the other hand, the FASTmrMLM method, while it detected the lowest number of QTNs detected, was considered the most efficient of the methods evaluated. Several studies comparing the ISIS-EM-BLASSO, pLARmEB, mrMLM, and FASTmrMLM methods have already been performed in many crops and the results are similar to those observed in the present study (Ma et al., 2018;Misra et al., 2018;Zhang et al., 2018;Fang et al., 2020). Although the ML-GWAS methods have similar approaches, the differential identification of QTNs is related to different screening and estimation models of each method . From the present study results, FASTmrMLM can be considered the most reliable method, as it presented a low rate of false-positive associations. The FASTmrMLM method results from an improvement of the mrMLM method, which is a faster, more reliable, with higher statistical power, higher estimation accuracy, and low falsepositive rate (Tamba and Zhang, 2018).
Most of the QTNs identified in this study were observed in only one environment, indicating the presence of frequent QTN x Environment interactions. Several studies have reported previously the presence of these morpho-agronomic traits interaction in common beans, suggesting that the gene expression of these QTNs is influenced by the evaluation environment (MacQueen et al., 2020;Wu et al., 2020). The presence of the QTN x Environment interaction is considered one of the main challenges in selecting QTNs in breeding programs, as these QTNs are more prone to environmental effects. On the other hand, the stable QTNs provided a remarkable demonstration of gradual improvement of all the traits evaluated in this study through the accumulation of favorable alleles.
FIGURE 3 | Accumulation of loci with favorable alleles in relation to adjusted means (LSmeans) for morpho-agronomic traits of common beans detected in accessions belonging to the Brazilian Diversity Panel (BDP). PH, plant height (cm); FPIH, first pod insertion height (cm); NN, number of nodules; PL, pod length (cm); NPP, total number of pods per plant; LP, number of locules per pod; SP, number of seeds per pod; TSW, total seed weight per plant (gm); W100, 100-seed weight (gm); YLD, grain yield (kg.ha −1 ).
Significant increases in productivity (YLD) and its primary components (NPP, SP, TSW, W100) were observed, and alleles that caused a significant PH reduction were identified. The PH is an essential factor in the formation of production components, and at the same time, it promotes or inhibits other components, affecting mainly the resistance to lodging, NN and NPP (Fang et al., 2020). Small-sized plants are associated with a determinate growth habit, less susceptibility to lodging, and shorter cycles (Chang et al., 2018). Over the years, with the technification of agriculture, common bean breeding programs have sought to develop plants with these traits, as they facilitate management and mechanized harvesting, reducing harvest losses and susceptibility to some diseases, allowing an increase in the number of crops per year due to a reduction in the cycle (Teixeira et al., 1999).
Most QTNs identified in this study had a small effect, confirming the complex and quantitative nature of the main morpho-agronomic traits in common beans (Gupta et al., 2020;MacQueen et al., 2020). As most of the traits studied are controlled by polygenes, the effect of each locus individually is relatively small. Nevertheless, it is vital to identify small-effect loci that cumulatively can explain the variation in a trait (Nakano and Kobayashi, 2020). The selection of higher-effect QTNs is preferable for the selection assisted by molecular markers (SAM) (Nadeem et al., 2018;Oladosu et al., 2019). However, the use of small-effect QTNs associated with the traits of interest is considered an important strategy in approaches to genomic selection (GS) since only these QTNs can replace the need for high-density genotyping by random SNPs and thus reduce genotyping costs. Moreover, models of GS using only SNPs known to be associated with the traits of interest showed greater accuracy of prediction since they showed lower background noise in constructing these models (He et al., 2019a;Ali et al., 2020).
Among the candidate gene models identified in this study, nine (Phvul.001G189200, Phvul.001G192200, Phvul.003G039900, Phvul.006G098300, Phvul.007G246700, Phvul.008G013300, Phvul.008G268700, Phvul.008G277352, and Phvul.011G020500) were previously identified in other studies of GWAS for morpho-agronomic traits in common beans Moghaddam et al., 2016;Soltani et al., 2017;Tock et al., 2017;MacQueen et al., 2020). The candidate gene model Phvul.003G039900, associated with W100 in the present study, was also identified for seed weight by MacQueen et al. (2020). The same authors observed an association between the gene model Phvul.006G098300 and PH, whereas this gene was associated with FPIH in the present study. The gene model Phvul.003G039900 has a putative methyltransferase activity function and the gene model Phvul.006G098300 is related to transferase activity and transferring acyl groups other than amino-acyl groups. The candidate gene model Phvul.008G013300, related to PL in this study, was also associated with the weight of seeds as reported by Moghaddam et al. (2016) and has a serine-type endopeptidase activity and proteolysis functions. The candidate gene model Phvul.011G020500 was associated with PH, NPP, and TSW, while this same gene model was associated with the aerial part's biomass trait in the observations of Soltani et al. (2017). Several functions were reported for this gene, such as DNA-binding transcription factor activity, transcription regulator complex, regulation of transcription and cell cycle.
Considering that the genomic regions around the significant QTNs for the different traits assessed in this study overlapped one another, 74% of the identified genes were also detected for more than one trait. Due to the strong LD observed in common beans, it is difficult to assign a gene precisely to a trait, especially when it comes to polygenic traits (Ikram et al., 2020). The genes located in the genomic regions around the identified QTNs may serve as promising targets for studying molecular mechanisms responsible for morpho-agronomic traits in common beans.
In this study, 64 QTNs were identified for ten morphoagronomic traits in common beans. Thirty-nine of them were identified as favorable alleles that can significantly increase trait expression and potentially yield and its components in the cultivation of common beans through allele pyramiding. The results reinforce the importance of conducting phenotyping of individuals in multiple environments, using multiple detection methods to increase the reliability of QTNs obtained in GWAS studies. The QTNs identified proved adequate for implementation in common bean breeding programs, mainly for improving Mesoamerican common beans from the Black and Carioca commercial classes, which are the primary targets in Brazil.

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

AUTHOR CONTRIBUITIONS
JD, VM-C, and LG conceived and designed the study. JD collected plant material, extracted DNA, and performed the genotyping. JD, JS, AN, LR, and DZ performed the phenotyping. JD and LG performed bioinformatics and statistical analyses. JD and DZ drafted the manuscript. JD, VM-C, JS, DZ, PR, PG, and LG edited and revised the final manuscript. All authors read and approved the final manuscript.

ACKNOWLEDGMENTS
We would like to thank the Instituto de Desenvolvimento Rural do Paranaì (IDR-Paranaì) and the University of California, Davis (through the Gepts' Lab) for support this research and the Coordenação de Aperfeiçoamento de Pessoal de Niìvel Superior-Brasil (CAPES) for the scholarship to JD in Brazil and abroad (Finance Code 001).