The Mesoamerican common bean diversity panel (BDP) used in this study was developed to gather accessions adapted to tropical conditions, reflecting the diversity of Brazil’s domesticated Mesoamerican gene pool (Delfini et al., 2021a). This panel includes cultivars, landraces, and breeding lines, all conserved in the gene bank of the Instituto de Desenvolvimento Rural do Paraná IAPAR-EMATER (IDR-Paraná). This study analyzed a subset of the Brazilian diversity panel of common beans (BDP), including 178 Mesoamerican common bean accessions.

The experiments were conducted during the 2018 rainy season at the research stations of the IDR-Paraná in the municipalities of Londrina (23°22’8”S; 51°9’48” W) and Guarapuava (25°23’51”S; 51°32’36” W), located in the Paraná state, Brazil. The experimental design adopted was a randomized block design with four replications. Each plot consisted of four rows, 2.00 m in length, with 0.5 m spacing between rows and a density of 12 plants per linear meter. Pest, disease, and weed control was carried out according to crop recommendations. After physiological maturity, each experimental plot was harvested, and a sample of 100 g of grains without physical or insect damage and free of disease was taken. These samples were stored in a cold chamber at a temperature of 5.6°C and 33% humidity until laboratory analyses were performed.

Samples of 50 g of grains from each experimental plot were ground in a mill (Perten 3100) to obtain a homogeneous powder using a 60-mesh sieve. The bean flour samples were packaged and stored at -18°C until analysis. The moisture content of the ground grains was determined by drying at 105 ± 3°C until constant weight, using triplicate repetitions of approximately 2 grams each, with the results expressed as a percentage (Brasil, 2009).

For the determination of total phenolic content (TPC), total flavonoid content (FLA), and antioxidant activity (DPPH), a metabolic extract was obtained following the methodology of González and Vázquez (2008), with minor modifications. The extract was prepared by extracting 1.0 g of fresh samples with 10 mL of 70% ethanol (v/v), homogenized in a blender (500 W) for 20 seconds, followed by centrifugation at 1013 x g (Excelsa 2 Fanem model 205N) for five minutes.

To determine TPC, 1.0 mL of the metabolic extract was mixed with 1.0 mL of methanol, 1.0 mL of 0.2 N Folin-Ciocalteu reagent, and 1.0 mL of 10% sodium carbonate (w/v). The mixture was allowed to rest for 30 minutes in the dark at room temperature (25°C). Subsequently, absorbance was measured at a wavelength of 765 nm using a Micronal spectrophotometer (AJX1600). Gallic acid was used as a standard at 10.0–100.0 mg L-1 concentrations, and results were expressed as mg of gallic acid equivalents per 100 g of sample (mg GAE 100 g-1) (Swain and Hillis, 1959).

FLA was determined by mixing 1.0 mL of the metabolic extract with 1.0 mL of 5.0% aluminum chloride (w/v) and 2.0 mL of methanol, allowing the mixture to rest for 30 minutes in the dark. Absorbance was then measured using a spectrophotometer (Micronal, AJX-1600) at 425 nm. Quercetin was used as a standard at concentrations of 50.0–500.0 mg L-1, and results were expressed as mg of quercetin equivalents (QE) per 100 g of sample (mg QE 100 g-1) (Lee et al., 1995).

For determining antioxidant capacity (DPPH) by scavenging the 2,2-diphenyl-1-picrylhydrazyl (DPPH•) radical, 50.0 μL of the metabolic extract was mixed with 1.0 mL of 100 mM acetate buffer (pH 5.5), 1.0 mL of methanol, and 0.5 mL of ethanolic DPPH solution (250.0 μM). The tubes were kept at room temperature in the dark for 15 minutes. The absorbance of the remaining DPPH• radical was measured using a spectrophotometer (Thermo-Genesys) at a wavelength of 517 nm. The calibration curve was prepared using Trolox (6-hydroxy-2,5,7,8-tetramethylchroman-2-carboxylic acid) as a standard, and results were expressed as μmol of Trolox equivalent antioxidant capacity (TEAC) per 100 g of fresh sample (Brand-Williams et al., 1995).

Phenotypic data were subjected to variance analysis, and heritability was subsequently determined. Histogram, density distribution, and boxplot graphs were generated for the groups of bean accessions evaluated, categorized by grain color (carioca, black, and colored) and genetic material type (cultivars, lines, and landraces). T-test was used to evaluate the difference between the bean groups. All analyses were performed using R software (R Core Team, 2024) (https://www.r-project.org), using the packages ggplot2 (Wickham and Wickham, 2016), FactoMineR (Lê et al., 2008), and cluster (Maechler, 2019) packages.

The genotyping-by-sequencing (GBS) technique was used to obtain the SNPs. The methodology used, as well as the results 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) and the data were imputed using Beagle software version 5 (Browning et al., 2018), and only SNPs anchored to chromosomes in the common bean reference genome (Schmutz et al., 2014) were used.

A total of six multi-locus genome-wide association studies (GWAS) were implemented in the MLM.GUI 4.0 software (Wen et al., 2018) was used to detect significant QTNs for the target traits: mrMLM, FASTmrMLM, FASTmrEMMA, pLARmEB, pKWmEB, and ISIS-EM-BLASSO. Population structure (Q) and kinship matrix (K) were included in the model to reduce false positives and enhance analytical power. STRUCTURE results (K = 2) were used for Q, and a kinship matrix was calculated using the mrMLM.GUI 4.0 software. For all methods, default parameters were used, and a logarithmic odds score (LOD) ≥ 3 was used as the critical value for significant associations. The phenotypic data used for GWAS were the adjusted means from each of the two environments. To obtain more reliable results, only QTNs that showed repeatability (i.e., detected at least twice using different methods or in different environments) were considered significant and were used in searches for favorable alleles and candidate genes.

For each QTN, all accessions in the panel were initially divided into two groups based on the QTN genotype, and alleles associated with a positive effect on the phenotype (i.e., those associated with an increase in the analyzed trait) were identified. Subsequently, the number of favorable alleles for each accession was determined, and the association between the accumulation of these superior alleles in the same accession and trait was verified. Boxplots were generated using R software for better visualization of the results.

The search for potential candidate genes focused on QTNs detected using multiple methods or environments. The search radius (physical distance) was limited according to the decline of half the LD corrected for population structure and relatedness (r2vs). Genes within the association region were identified based on the reference genome annotation of common bean Phaseolus vulgaris v.2, published on the Phytozome v10.3 website (http://phytozome.net). Subsequently, genes with putative functions related to the traits of interest were selected as candidate genes, as indicated by Gene Ontology (GO) annotation (http://www.geneontology.org/).

The analysis of variance showed significant effects for genotypes and genotype × environment interactions (p<0.01) for all traits evaluated, while no significant effect was observed for the environment ( Table 1 ). The experimental coefficient of variation ranged from 3.67 to 6.40%, and heritability ranged from 93.1 to 99.8%, indicating a low influence of environmental variation. Significant correlations (p<0.01) were observed among traits, with values of 0.74 (TPC × FLA), 0.48 (TPC × DPPH), and 0.50 (FLA × DPPH).

No differences were observed among the cultivar groups (cultivars, lines, and landraces) for TPC ( Figure 1 ). Regarding the color groups (Carioca, black, and colored), the Carioca group had the highest average value (165.70 mg GAE 100 g^-¹), followed by the colored and black types, with 154.31 and 132.67 mg GAE 100 g^-¹, respectively. For the carioca group, the highest TPC values were found in the genotypes FEB200, FEB178, TAA Gol, Carioca1070, and LP07 (332.42, 290.17, 277.83, 249.90, and 249.36 mg GAE 100 g^-¹, respectively), while in the black group, the highest values were found in BAT76, DOR445, ARC2, Awauna, and BlackHawk (354.36, 309.73, 272.58, 228.78, and 210.16 mg GAE 100 g^-¹, respectively).

Total flavonoid content (FLA) ranged from 108.55 to 712.73 mg QE 100 g^-¹, with no differences observed among the cultivar groups. The carioca group (422.18 mg QE 100 g^-¹) had higher average FLA values than colored and black beans (371.47 and 315.57 mg QE 100 g^-¹, respectively). In the carioca group, the highest FLA values were observed in the genotypes Carioca1070, LP15, FEB200, LP14, and LP13 (653.34, 626.26, 615.80, 606.14, and 603.07 mg QE 100 g^-¹, respectively), while in the black group, they were DOR445, BAT76, BlackHawk, ARC2, and LP09 (687.66, 671.44, 602.26, 574.89, and 567.31 mg QE 100 g^-¹, respectively).

The average value for antioxidant activity (DPPH) was 803.80 μmol TEAC 100 g^-¹ (range: 444.91–944.34 μmol TEAC 100 g^-¹). The average values for landraces, breeding lines, and cultivars were 788.11, 816.53, and 797.92 μmol TEAC 100 g^-¹, respectively, and for the carioca, black, and colored grain groups, they were 813.75, 807.36, and 776.65 μmol TEAC 100 g^-¹, respectively. The highest antioxidant activity values were observed in the carioca group genotypes BRS Ametista, Aporé, IAPAR 72, LP31, and FEB178 (936.41, 931.53, 893.61, 893.32, and 893.28 μmol TEAC 100 g^-¹, respectively), and in the black group genotypes Diamante Negro, BRS Campeiro, BRS Esplendor, ICAQuetzal, and LP29 (944.34, 923.41, 903.60, 894.75, and 892.72 μmol TEAC 100 g^-¹, respectively).

Linear regression analysis showed no significant effects (p>0.05) for all traits in both the carioca and black groups, indicating a lack of genetic progress ( Figure 2 ). The coefficient of determination (R²) was 0.00 for most traits.

In total, 147 unique QTNs were identified as significant for the three biochemical traits studied (TPC, FLA, and DPPH), with 41 of these QTNs being detected at least twice by one of the seven multi-locus methods used or in different environments analyzed (LDA, GUA, and LSmeans), as shown in Table 2 and Supplementary Figure S1 . Table 2 presents the QTNs that demonstrated repeatability. The analyses revealed variation in the number of QTNs associated with each trait: 11 QTNs for TPC, 17 for FLA, and 13 for DPPH. One pleiotropic QTN was identified, showing a significant association at least twice for more than one trait, specifically shared between TPC and FLA. In addition to this pleiotropic QTN, seven overlapping QTNs were found within a genomic region defined by linkage disequilibrium (LD) of ± 296 kb. Two QTNs associated with FLA and DPPH were located on chromosome Pv2, where an overlap occurred. Other overlapping QTNs, all associated with DPPH, were identified on chromosome Pv06. Additionally, on chromosome Pv11, overlapping QTNs were observed for all three traits: TPC, FLA, and DPPH.

QTNs were identified on all 11 common bean chromosomes, with chromosomes Pv02 and Pv11 having the highest number of QTNs, with five each. The phenotypic variation explained (PVE) by these QTNs varied significantly, ranging from very low values, close to 2%, to 16.37%. Although some QTNs showed negligible effects, the observed effects varied from -135.18 to 79.81.

Different PVE and QTN effect results were obtained for the same QTN using different methods, with more significant results for PVE. For instance, the FASTmrMLM method identified some QTNs with effects and PVE close to zero. ISIS EM-BLASSO identified the highest number of significant SNPs among the six methods used, followed by pLARmEB, FASTmrMLM, mrMLM, pKWmEB, and FASTmrEMMA. However, when focusing only on QTNs that demonstrated repeatability, the mrMLM and FASTmrMLM methods stood out, followed by FASTmrEMMA, ISIS EM-BLASSO, pKWmEB, and pLARmEB. The FASTmrMLM method achieved 99% efficiency, indicating that only 1% of the SNPs detected by this method did not show repeatability. In contrast, the mrMLM, ISIS-EM-BLASSO, and pLARmEB methods achieved efficiencies of 70, 52, and 10%, respectively. Although the ISIS EM-BLASSO method initially identified the highest number of SNPs, many of these SNPs were discarded in subsequent stages of the study.

The Guarapuava environment identified the highest initial number of significant SNPs, with 53, followed by Londrina and LSmeans, with 49 and 45, respectively. However, when considering the more stable SNPs, Guarapuava had the highest number, followed by LSmeans, while the Londrina environment had the lowest number of significant SNPs.

Favorable alleles that demonstrated repeatability were identified in each QTNs, meaning alleles were associated with increased biochemical content. Subsequently, an analysis was conducted to investigate whether the accumulation of these alleles corresponded to an increase in TPC, FLA, and DPPH content ( Figure 3 ). The results indicated a progressive increase in TPC, FLA, and DPPH levels in bean grains, proportional to the number of favorable alleles in the genotype.

Initially, genes in the genomic regions identified around the QTNs were searched. Some QTNs were located within the candidate genes, while others were at varying distances, ranging from 6.6 to 287 kbp (kilobase pairs) from the QTNs. This highlights the intricate and multifaceted nature of genetic regulation, emphasizing the need for further research and exploration to fully understand genetic associations.

The identified gene functions were associated with various processes, including redox reactions, protein phosphorylation, transcriptional regulation, metabolic processes, transmembrane transport, and metal ion transport, reflecting the diversity of biological functions of genes linked to QTNs. All the biochemical compounds studied had at least one candidate gene whose function was related to their composition, as indicated in Table 3 . The highest number of candidate genes was found for DPPH, followed by FLA, with 8 and 2 candidate genes, respectively, while only one gene was associated with TPC.