Skip to main content

GWAS of pod morphological and color characters in common bean



Common bean (Phaseolus vulgaris L.) is an important legume species which can be consumed as immature pods and dry seeds after re-hydration and cooking. Many genes and QTL, and epistatic interactions among them, condition pod morphological traits. However, not all them have been mapped or validated nor candidate genes proposed. We sought to investigate the genomic regions conditioning pod morphological and color characters through GWAS.


Single and multi-locus genome wide association analysis was used to investigate pod traits for a set of 301 bean lines of the Spanish Diversity Panel (SDP). The SDP was genotyped with 32,812 SNPs obtained from Genotyping by Sequencing. The panel was grown in two seasons and phenotypic data were recorded for 17 fresh pods traits grouped in four pod characters: pod length, pod cross-section, pod color, and number of seeds per pod. In all, 23 QTL for pod length, 6 for cross-section, 18 for pod color, 6 for number of seeds per pod and 9 associated to two or more pod characters were detected. Most QTL were located in the telomeric region of chromosomes Pv01, Pv02, Pv04, Pv08, Pv09 and Pv10. Eighteen detected QTL co-localized with 28 previously reported QTL. Twenty-one potential candidate genes involving developmental processes were detected underlying 11 QTL for pod morphological characters, four of them homologous to A. thaliana genes FIS2, SPL10, TTG2 and AML4 affecting silique size. Eight potential candidate genes involved in pigment synthesis, were found underlying five QTL for pod color.


GWAS for pod morphological and color characters in the bean Spanish Diversity Panel revealed 62 QTL, 18 co-localized with previously reported QTL, and 16 QTL were underlain by 25 candidate genes. Overall 44 new QTL identified and 18 existing QTL contribute to a better understanding of the complex inheritance of pod size and color traits in common bean and open the opportunity for future validation works.

Peer Review reports


Common bean (Phaseolus vulgaris L.) is an important legume species domesticated in two different areas of Latin America representing distinct Mesoamerican and Andean gene pools [1]. Cultivated genotypes of common bean exhibit wide diversity for growth habit, flower color, and shape, size and color of pods and seeds. Immature pod phenotypic diversity involves variation in length and curvature (straight vs curved), cross section (diameter, flat, round, sieve size), and color (yellow, green, purple) before the seeds start to develop. Immature pods of some bean genotypes are consumed as fresh green beans (syn. Garden, green, pole, snap, haricot or French beans) when the pods have reached maximum length while the seed is still forming, in contrast with dry beans that are consumed as mature seeds after re-hydration and cooking. In all, 2.29 Mha were destined to snap bean crop in 2019 while 33.8 Mha were used in the dry bean crop ( The snap bean group includes different market classes such as ‘string snap bean’ referring to types where the pod suture strings must be removed before consumption; ‘yellow wax’ and ‘green bean’, referring to yellow and green pod, respectively; ‘Romano type’ with a very large and flat pod; and ‘blue lake type’ with dark green pods that remain stringless and fibreless [2]. Furthermore, snap bean can be classified according to processing adaptation: frozen, canned, or fresh market.

Different studies have reported on the genetic control of pod morphological characters. Classical genetic studies in common bean described major genes controlling the cross section (Ea and Eb genes [3, 4]), pod membrane (Fa, Fb, and Fc genes [3, 4]), parchment pod (Ia, Ib genes [3,4,5];), stringless pod (St gene [6];), twister pod (Tw gene [7]) and straight pod (Da, Db [3, 4]). St was mapped to chromosome Pv02 [8], near the common bean ortholog of PvIND, a gene controlling pod dehiscence [9]. Pod shattering (dehiscent pod) is an important trait associated with seed dispersal which was modified to indehiscent pod during domestication [8]. Recent studies indicate that two major quantitative trait loci (QTL) located on chromosomes Pv03 (PvPdh1 [10]) and Pv05 (qPD5.1-Pv [11]) also influence the pod shattering trait.

A few studies report on the quantitative inheritance of pod length, thickness and width, and identification of QTL controlling these traits mapped across all 11 bean chromosomes [12,13,14,15]. Hagerty et al., used a dry bean x snap bean recombinant inbred population, to map: St (pod suture string) to Pv02; overlapping pod wall fiber, width, and thickness to Pv04; and pod length to Pv09 [14]. Murube et al., using two nested populations, found four genomic regions located on chromosomes Pv01, Pv02, Pv07 and Pv11 with overlapping QTL for pod size characters and number of seeds per pod [15].

The genetic control for color of immature pods is influenced by the Y and Arg genes: Y Arg exhibits green pod, y Arg yellow wax pod, Y arg greenish gray (silvery) pod, and y arg white pod [4]. The y allele conferring yellow pod color was mapped to Pv02 by Koinange et al. [8]. The B gene which regulates the production of precursors of anthocyanins pathway above the level of dihydrokaempferol formation also resides on Pv02 [16]. The genes Pur and Ro influence a range of pod colors from rose to purple pods [17]. The Ace gene produces shiny pod [18]. Myers et al. identified quantitative trait nucleotides (QTNs) associated with CIE L*, a*, b* color space values for pod color on Pv02, Pv03 and Pv05 in a panel of 149 snap bean accessions [19].

In summary, many genes and QTL, and epistatic interactions among them, condition pod morphological traits. However, not all genes have been mapped nor candidate genes proposed. Moreover, QTL need to be validated in different genetic backgrounds and environments before they can be implemented directly in plant breeding or used to search for underlying candidate genes. The reference genome for P. vulgaris [20] provides the framework for fine mapping genes and QTL conditioning pod morphological traits and to identify candidate genes. The reference genome combined with high throughput genotyping, improving statistical programs for detecting marker - trait associations, and access to diversity panels which have greater variation than bi-parental populations, enhances opportunities to identify putative genomic regions controlling specific traits [21].

A Spanish common bean diversity panel (SDP) of 308 lines was established from the local Spanish germplasm collection that included landraces and old and elite cultivars used for pod consumption [22]. The main aim of this work was to investigate genomic regions controlling pod size and color traits through genome wide association analysis (GWAS) of the SDP. Results will contribute to discovery of new genomic regions associated with pod characters, validation of reported QTL, and identification of candidate genes for the investigated traits.


Phenotypic variation, correlations and heritability

A total of 301 SDP lines were successfully characterized for the 17 morphological traits. The results show a wide and continuous variation for the 16 quantitative traits evaluated (see Figure S1 and Table S1). For instance, PL and NSP, two traits related to yield, ranged between 7.1 to 26.4 cm and 2.2 to 8.3 seeds, respectively. The SDP exhibited wide variation for color with green (241 lines), yellow (38), purple (3), green mottled (16) and yellow mottled (2) pods. Pod color measured by the CIE scale exhibited wide variation as well for the L*, a*, and b* vectors. For example, b* varied from − 4.38 to 40.8. The H2 estimations for the 16 quantitative traits were high ranging from 0.31 for PSW to 0.91 for PLW (see Table S1).

Correlation analyses indicated significant relationships between many evaluated traits (Fig. 1). There were significant and positive correlations among the six pod section traits and a significant negative correlation for PSH/PSW, PSC and PSW. Most of the six pod length variables were significantly correlated except PL/PLC with PLP and PLA. Correlation analyses also revealed significative correlation among section and length traits except in five cases; PSC with PLP, PL and PLC and PSW with PLW and PL/PLC. NSP was significantly correlated with four pod length traits (PLA, PLP, PL and PLC). Finally, the three pod color variables (L*, a*, b*) were also significantly correlated.

Fig. 1

Corrplot showing the Pearson correlation among the 16 quantitative pod traits evaluated (see Table 1). Non-significant correlations (α = 0.05) are indicated with X

Characterization and detection of SNPs

Sequencing of the GBS libraries yielded approximately 418 million reads in total for the 301 SDP lines. About 76.3% of the reads were successfully aligned to the common bean reference genome, 21.5% of the reads mapped to more than one locus, and 23.7% were unmapped. The NGSEP genotyping pipeline produced 346,819 biallelic SNPs in the 11 chromosomes and scaffolds of the reference genome. 32,812 SNPs distributed across the eleven bean chromosomes were retained after filtering parameters (Figure S2). Most of these SNPs were present in coding regions (51.1%) and represented 46.1% silent mutations, 32.6% missense, 5.4% non-sense and 15.9% prime UTR regions. While intronic and intergenic regions contained 31.2 and 17.8%, respectively. A genome-wide transition/transversion (Tr/Tv) ratio of 1.17 was observed.


SL-GWAS (MLM) revealed 63 significant QTNs, 57 of them grouped in 9 genomic regions (QTI): 7 for pod length, 1 for cross-section, and 1 for pod color. Six QTN showed a single association (Table S2). QTNs were not detected for PCOL and NSP. Interestingly, twenty-eight QTNs for pod morphological traits were detected at the distal end of chromosome Pv01 (45,582,871 to 48,454,962) and 9 QTNs for the color vector b* in the telomere of chromosomes Pv07 (32026373–32,413,401).

ML-GWAS, using the six multi-locus models in the mrMLM package, revealed 103 QTN (Tables S3, S4, S5). QTNs were not detected for the index PSC. The mrMLM method detected the most associations (37) while FASTmrEMMA detected the fewest associations (21). In all, 14 significant QTNs were found for pod section traits (Table S3) and the QTN number per character ranged from 5 for PSH and only one for PSA and PSW. For pod length traits, 52 QTN were detected (Table S4), with 18 of them identified by at least two different GWAS methods. The number of QTNs ranged from 2 for PL/PLC a and 10 for PLP, PLA, PLC and NSP. These QTNs were mostly located in the telomeric regions of Pv01 and Pv02. Concerning pod color measured by CIE space, a total of 27 QTNs were detected (11 for the vector L*, 9 for a* and 7 for b*) while 10 QTNs were detected for pod color measured visually as a qualitative character (Table S5). These QTNs were mostly located in telomeric regions of chromosomes Pv02 (7) and Pv08 (5). 31 QTNs revealed by ML-GWAS were grouped in eleven QTI. Three QTI were identified by both methods (SL-GWAS and ML-GWAS): Chr01:48090873–48,454,962; Chr02: 47302543–47,669,811; Chr07:32026373–32,413,401.

The 166 QTNs detected 62 QTL, 23 for pod length, 6 for cross-section, 18 for pod color and 6 for number of seeds per pod as well as 9 QTL associated with multiple characters (Table 1). Most QTL were located on chromosomes Pv02 (12), Pv04 (7), Pv08 (7), and Pv10 (11) (see Fig. 2).

Table 1 Quantitative trait loci (QTL) for pod morphological characters and pod color detected in this association study (GWAS). The description of specific associations trait-SNP are reported in Tables S2, S3, S4 and S5. QTL in bold indicates associations detected by different types of analysis (SL-GWAS & ML-GWAS)
Fig. 2

Circle plot showing the comparison of genomic positions for the QTL identified in this work (light blue track) with previously reported QTL for pod morphological traits (green track). Yellow track shows the positions of the meta-QTL detected in this study

Co-location of QTL

Genomic positions for 96 previously reported QTL [10, 12, 14, 15, 19] for pod morphological traits in common bean were examined for overlap with the QTL identified in this work (Fig. 2). There were 15 genomic regions where a reported QTL and QTL detected in this study for pod traits overlapped (Table 2; Fig. 2). These regions were located on seven chromosomes (Pv01, Pv02, Pv03, Pv04, Pv05, Pv06 Pv08 and Pv11). The beginning of chromosome Pv02 (542087–959,169) only co-located QTL for pod color, whereas the other overlapping QTL were associated with pod morphological traits or both morphological and pod color traits.

Table 2 List of detected QTIs and QTNs that showed a co-location with reported QTL

In silico genome exploration

In silico analysis of the annotated genes underlying the 62 QTL for pod traits revealed 25 potential candidate genes in 16 QTL (Table 3). There were 12 candidate genes underlying six QTL for pod size traits. Five of these candidate genes were Cytochrome P450 and five were WRKY or MYB transcription factors, proteins involved in multiple processes like responses to biotic and abiotic stresses, development, differentiation, metabolism, defense, and pigment synthesis [23, 24]. Seven of these genes were involved in controlling pod development: Phvul.001G229900, Phvul.001G221500, Phvul.002G016100, Phvul.004G144900, Phvul.006G076800, Phvul.006G077200 and Phvul.010G010200 [25]. Four genes, homologous with genes controlling silique size in A. thalina [26] were detected underlying QTL associated with both size and color traits: Phvul.001G262600 with SPL10 gene in QTL NSPCol01_51; Phvul.002G141800 with FIS2 gene near QTL PodL02_29.1; Phvul.006G074600 with TTG2 gene in QTLPodLSN06_18.4; and Phvul.008G019500 with AML4 gene in QTL NSPCol08_1.7. A DELLA protein gene Phvul.001G230500, controlling various aspects of plant growth and development, including flowering, and pod setting and development [27], is a candidate for QTL PodLCol01_48.4.

Table 3 Potential candidate genes underlying quantitative trait intervals (QTI) and single quantitative trait nucleotides (QTN) detected by SL-GWAS and ML-GWAS. *, candidate genes located in the external border of investigated regions

Concerning QTL for pod color traits, six candidate genes encoding Cytochrome P450 proteins underlie the QTL PodCol02_2.4, PodCol02_43.6, PodLCol02_47.6 and PodCol10_38.7, and five candidate genes encoding MYB TRANSCRIPTION FACTOR underlie the QTL PodLCol01_48, PodCol07_32 and PodCol08_60.2 (Table 3). Also, Phvul.001G261500 encoding a Flavonoid 3′-hydroxylase is a candidate gene for QTL PodColN01_51, and a cluster of genes encoding Flavone/flavonol 7-O-beta-D-glucoside malonyltransferase bordered QTL PodLCol08_2.7.


Pod morphology and color are important traits in common bean because they influence consumer preference for pods which are eaten as green beans for many genotypes. This study identified genomic regions controlling pod traits in the Spanish Diversity Panel. This panel encompasses wide genetic [22] and phenotypic variation for pod color, pod size, pod cross section, and number of seeds per pod (see Figure S1). For instance, variation in pod length ranged between 26.5 and 7.5 cm, for SDP203 a Romano type with a very large green pod and SDP138 with a very short and flat green pod, respectively. Pod color varied from green to yellow to purple and by quantitative classification (CIE scale). The experimental design used (randomized complete block with a repetition per season) may affect the accuracy of trait estimation, particularly in a large trial. However, the recorded traits have high heritability and for most traits, H2 estimates were high, suggesting a few major genes were involved. Results of correlation analysis support the grouping of the traits in four characters (pod length, pod cross section, and pod color traits, and seeds per pod). Most traits within a pod character were significantly correlated (Fig. 2). Number of seeds per pod (NSP), a major yield component [28], was only significantly correlated with pod length, pod area, pod perimeter and pod length curved.

Most existing methods used in association studies are based on single marker association in genome-wide scans with population structure and consider stringent methods to control false positive rate [29, 30] so that some associations may not be detected by single locus models [29, 30]. ML-GWAS showed a total of 103 associations with pod traits (14 for pod cross section characters, 42 for pod length character, 10 for number of seeds per pod and 37 for pod color) while SL-GWAS revealed 63 associations (3 for pod cross section, 50 for pod length, 10 for number of seeds per pod and 10 for pod color). All these association were grouped in a 62 QTL; 23 QTL involved in pod length characters, 6 in pod cross section characters, 18 in pod color, 6 in NSP and 9 in two more characters (Table 1; Fig. 2).

We observed that 18 QTL were co-located with earlier described QTL for pod size in various populations. QTL located in the same position across different studies and populations supports QTL validation and disposition of robust QTL. Five remarkable chromosome regions for pod traits are detailed below:

  • Pv01 (50–51 Mb) where PodColN01_51 overlaps with QTL PP1.2XB, PL1XC* and PWI1XC* [15] for pod length. Within this region is Phvul.001G262600, a homologue to the Arabidopsis SPL10 gene (AT1G27370), which was proposed as candidate gene for silique length in Brassica napus [26]. The SPL genes are also implicated in the regulation of anthocyanin biosynthesis [31], which may explain why the same region possesses QTL involved in both pod color and pod size control.

  • Pv02 (0.54–0.95 Mb) where PodCol02_0.8 overlap with the reported QTL a*2.1, L*0.9 and L*2.1 [19] for pod color.

  • Pv02 (29.1–30.2 Mb) where PodL02_29.1 overlaps the reported QTL PL2.1XC for pod length [15]. Phvul.002G141800, a homologue of the Arabidopsis gene FIS2 (At2G35670) that represses seed development in the absence of pollination, is a candidate gene for pod length in this region. FIS2 was also related to silique size by Wang et al. [26].

  • Pv02 (48.6–49.6 Mb) where QTL NSP02_48.7 and PodLCol02_49.4 overlap the reported QTL PL2.2XC, NSP2XCand E-PP2XB for pod length, number of seeds per pod and pod perimeter [15].

  • Pv06 (18.3–19.3 Mb) where PodLSNS06_18.4 overlaps the QTL PWI6XB for pod width [15]. Phvul.001G173700, a homologue of the Arabidopsis gene TTG2 (At2G37260), which affects seed size and weight in Arabidopsis and underlies a QTL for silique length in Brassica napus [26], is a candidate gene in this region.

For pod color, more QTL were detected when measured as a quantitative variable (CIElab scale) than as a qualitative trait (5), suggesting that the former evaluation provides additional information. Six QTL were associated with both qualitative and quantitative pod color characters: NSPCol01_51, NSPCol08_1.7, PodCol02_0.8, PodCol04_47.8, PodCol04_7.2, PodCol10_43.4 and PodLCol02_47.6. Using a qualitative assessment, gene Y controlling yellow wax, was mapped to the proximal end of Pv02 [8]. A quantitative assessment detected significant QTL for a* and L* color variables in the same proximal location on Pv02 [19]. Herein PodCol02_0.2, PodCol02_0.8 and PodCol02_2.4 QTL were similarly located, supporting the relevance of this region to pod color. Myers et al. proposed Phvul.002G004400 (a pentatricopeptide repeat) as a candidate gene for the Y gene [19]. However, between PodCol02_0.8 and PodCol02_2.3 reside other genes with functions that could be involved with pigment synthesis such us Phvul.002G014700 and Phvul.002G014800 encoding for a Isoflavone 2′-hydroxylase, and Phvul.002G022800 and Phvul.002G022900 encoding a Cytochrome P450-Related protein. Together with chlorophylls and carotenoid, flavonoids are one of the major pigments in higher plants, and some of them can influence yellow coloring [32].

Cytochrome P450, one of the largest gene families in plants are involved in different cellular processes including the synthesis of pigments [23, 24, 33]. Two other Cytochrome P450 (Phvul.002G263700, Phvul.002G263900) genes associated with QTL PodCol02_43.6 further suggested a possible role for them in influencing pod color. In fact, the gene B, involved in the pigment production in seed coats, and tightly linked to gene I (conferring resistance to BCMV [34]), is mapped in a similar distal telomeric region of Pv02. This resistant locus was characterized and located in bean in the bean genome near Phvul.002G323200, Chr02:48805820–48,810,839 [35]. It has been reported that the actual genes influencing a trait were often up to 2 Mbps away from the peak SNP detected by GWAS [36]. Similarly, close to the QTL PodCol08_2.7 were 10 genes (2.567.331–2.636.603pb) with a flavone/flavonol 7-O-beta-D-glucoside malonyl-transferase function related with the pigment synthesis (Phvul.008G031900, Phvul.008G032000, Phvul.008G032100, Phvul.008G032200, Phvul.008G032400, Phvul.008G032450, Phvul.008G032501, Phvul.008G032551, Phvul.008G032600, and Phvul.008G032700). Bordering QTL PodCol07_32 and PodCol10_40.2 were candidate genes Phvul.007G206200 (MYB61 [26]), Phvul.010G117200 (ARF18 [3637]), Phvul.010G117100 (CYP78A9 [38]), and Phvul.010G118700 (NTS1 [39]) which have homology with genes involved in the control of silique (see Table 4). Finally, undelaying to the QTL PodLCol01_48, PodCol07_32 and PodCol08_60.2 were found genes codifying MYB transcription factors (Phvul.001G229900, Phvul.001G221500, Phvul.007G206200, Phvul.008G262700). MYB proteins are key factors in regulatory networks controlling development, metabolism including the synthesis of anthocyanins [40].

Table 4 List of the 17 pod traits analysed. The code assigned to each character is indicated in parentheses. * 1, measured from digital images; 2, manually measured

In summary, GWAS revealed new and known genomic regions with QTL influencing pod size, pod color and number of seeds per pod. The 44 newly identified regions involved in the genetic control of pod size or color should be verified in future genetic analysis. The eightteen regions overlapping to previously reported QTL provide relevant information for the development of breeding programs and genetic analysis focused on these characters.


Plant material

The Spanish Diversity Panel (SDP) of 308 bean lines was described by Campa et al. [22]. Briefly, the SDP includes: 220 landraces, mostly from the updated Spanish Core Collection; 51 elite cultivars, mostly cultivated in Europe for snap bean consumption; and 37 lines representing traditional old cultivars and well-known breeding lines. The sequenced bean genotypes, G19833 [20] and BAT93 [41] were included as representatives of the Andean and Mesoamerican gene pools. The panel exhibits wide phenotypic variation for pod traits (see Figure S3). The population structure and linkage disequilibrium, described previously by Campa et al. [22], indicates two main groups corresponding to the Andean and Mesoamerican gene pools and a third group with admixture of both gene pools.


The SDP was phenotyped in the greenhouse at Villaviciosa, Spain (43°2901 N, 5°2611 W; elevation 6.5 m) during two seasons (spring 2017 and 2018). Each year represented a single replicate of a single 1-m row plot including 8–10 plants per line. The experiment design was a randomized complete block. Standard agronomic practices for tillage, irrigation, fertilization, and weed and insect control were followed to ensure adequate plant growth and development. Phenotyping was conducted for a set of 17 pod traits grouped in four main pod characters: pod length, pod cross section, pod color and number seeds per pod (see Table 4). Fresh pods were harvested at the beginning of R8 stage when pods had reached maximum length and seeds began to enlarge. Twelve quantitative characters included pod longitudinal (PLP, PLA, PLW, PL, PLC, PL/PLC) and cross section (PSP, PSA, PSW, PSH, PSH/PSW, PSC) dimensions that were obtained from 10 scanned fresh pods per line with the help of Tomato Analyzer software [42]. The external pod color was visually recorded as green, yellow, mottled green, mottled yellow and purple. To record variation within phenotypic classes, the fresh pod color was also quantified with Tomato Analyzer software measuring three vectors in the CIE scale: L* detects the brightness from 0 (black) to 100 (white), a* represents color from green (negative values) to red (positive values), and the b* measures blue (negative values) to yellow (positive values). In parallel, the fresh pod color was visually recorded as green, yellow, mottled green, mottled yellow and purple. Finally, the number of seeds per pod (NSP) was manually recorded as an average from 10 pods.

Mean values were adjusted identifying outliers through the coefficient of variation (CV). CV over 25% were not accepted. The phenotypic variation for individual traits was visualized by frequency distributions generated by ggplot2 [43]. Pearson’s correlation coefficients among the traits were also investigated using the package ‘corrplot’ [44]. The broad-sense heritability (H2) for each trait was estimated using the package ‘heritability’ [45]. The heritability and subsequent statistical analyses of the phenotypic data were conducted in R platform [46].


Genotyping-by-sequencing (GBS), as described by Elshire et al. [47], was conducted at BGI-Tech (Copenhagen, Denmark) using the ApeKI restriction enzyme. A GBS sequencing library was prepared by ligating the digested DNA to unique nucleotide adapters (barcodes) followed by PCR with flow-cell attachment site tagged primers. Sequencing was performed using Illumina HiSeq4000 and 100x Paired-End. The sequencing reads from different genotypes were deconvoluted using the barcodes and aligned to the Phaseolus vulgaris L. v2 reference genome (!info?alias=Org_Pvulgaris).

SNP discovery and genotype calling were conducted using NGSEP-GBS pipeline [48, 49]. Maximum base quality score was set to 30 and minimum quality for reporting a variant was set to 40. All SNP markers detected with less than 50% missing values and a minor allele frequency (MAF) 0.05 were retained to perform imputation with ImputeVCF module into NGSEP, which is a reimplementation of the Hidden Markov Model (HMM) implemented in the package fastPHASE ( Annotation of variants was performed using the command Annotate by NGSEP. The distribution of the SNPs along chromosomes was visualized with the CMplot package ( of the R project [46]. SNPs were named considering physical position in the bean genome: chromosome and genomic position (bp).

Genome-wide association analysis

Association analyses were carried out using both single-locus-GWAS (SL-GWAS) and multi-locus-GWAS (ML-GWAS) models for all traits. SL-GWAS was conducted in Tassel V5.1 [50] using the mixed linear model (MLM) approach with the PCA (3) and Kindship matrix as cofactor. ML-GWAS was performed with the mrMLM v4.0 package ( representing six different statistical models for traits with multi and polygenic effect (mrMLM; FASTmrMLM; ISIS EM-BLASSO; FASTmrEMMA; pLARmEB; pKWmEB). Five PCs generated from GAPIT were included as covariates and an identity-by-state kinship matrix was created using the Efficient Mixed Model Association (EMMA) algorithm implemented in GAPIT R package [51].

Association analysis were carried out in the three data set; two seasons (spring 2017 and 2018) and the mean of two seasons. Critical threshold of significance was -log(p) > 5 for SL-GWAS and LOD > 5 for ML-GWAS. Significant trait-SNP (QTN) associations were considered when detected in the three analysis. Quantitative trait intervals (QTI) were defined when several QTNs were located at distance less than 0.3 Mbp. Significant QTN were classified as QTL according to the pod character (PodL, PodS, PodCol and NSP) and named considering the genomic position (chromosome and position Mbp).

QTL alignment

For QTL alignments, published mapping data from four independent studies that reported QTL for pod morphological traits in common bean [10, 12,13,14,15, 19] were considered. Physical position was used to investigate the correspondence between the genomic regions identified in this work with the previously reported QTL. The physical position of QTL from the literature were based on flanking or underlying markers which were aligned with the bean reference (G19833) genome sequence v2.1 using the BLASTN algorithm ( Marker sequences were obtained from PhaseolusGenes ( or tag sequences containing the SNP supplied by the GBS analysis. ShinyCircos package [52] was used to visualize the position of each QTL in the bean genome from the underlying markers.

Candidate genes mining

Potential candidate genes were investigated in the bean genome v2.1 ( through exploration of the functional annotation of the genes underlying the detected QTL. In the case of single QTN, a window ±75,000 bp from the QTN position was considered. Genes with a known function in developmental processes were considered. In addition, homologous genes to genes involving in the control of silique traits in Arabidopsis thaliana model species were examined [26, 36].

Availability of data and materials

All data generated or analysed during this study are included in this published article [and its supplementary information files] or are available from the corresponding author on reasonable request.


  1. 1.

    Bellucci E, Bitocchi E, Rau R, Rodriguez M, Biagetti E, Giardini A, et al. Genomics of origin, domestication and evolution of Phaseolus vulgaris. In: Tuberosa R, Graner A, Frison E, editors. Genomics of Plant Genetic Resources. Dordrecht: Springer; 2014. p. 483–5007.

    Google Scholar 

  2. 2.

    Wallace L, Arkwazee H, Vining K, Myers JR. Genetic diversity within snap beans and their relation to dry beans. Genes. 2018;9(12):587. Published 2018 Nov 28.

    CAS  Article  PubMed Central  Google Scholar 

  3. 3.

    Lamprecht H. Beitrage zur Genetik von Phaseolus vulgaris. II. Über Vererbung von Hülsenfarbe und Hülsenform. Hereditas. 1932;16:295–340.

    Article  Google Scholar 

  4. 4.

    Lamprecht H. The inheritance of the slender-type of Phaseolus vulgaris and some other results. Agri Hort Genet. 1947;5:72–84.

    Google Scholar 

  5. 5.

    Lamprecht H. Weitere koppelungsstudien an Phaseolus vulgaris mit einer Übersicht über die Koppelungs gruppen. Agri Hort Genet. 1961;19:319–32.

    Google Scholar 

  6. 6.

    Prakken R. Inheritance of colors and pod characters in Phaseolus vulgaris L. Genetica. 1934;16(3-4):177–294.

    Article  Google Scholar 

  7. 7.

    Baggett JR, Kean D. Inheritance of twisted pods in common bean (Phaseolus vulgaris L.). J Amer Soc Hort Sci. 1995;120(6):900–1.

    Article  Google Scholar 

  8. 8.

    Koinange EMK, Singh SP, Gepts P. Genetic control of the domestication syndrome in common bean. Crop Sci. 1996;36(4):1037–45.

    Article  Google Scholar 

  9. 9.

    Gioia T, Logozzo G, Kami J, Spagnoletti Zeuli P, Gepts P. Identification and characterization of a homologue to the Arabidopsis INDEHISCENT gene in common bean. J Hered. 2013;104(2):273–86.

    CAS  Article  PubMed  Google Scholar 

  10. 10.

    Parker TA, Berny MY, Teran JC, Palkovic A, Jernstedt J, Gepts P. Pod indehiscence is a domestication and aridity resilience trait in common bean. New Phytol. 2020;225(1):558–70.

    CAS  Article  PubMed  Google Scholar 

  11. 11.

    Rau D, Murgia ML, Rodriguez M, Bitocchi E, Bellucci E, Fois D, et al. Genomic dissection of pod shattering in common bean: mutations at non-orthologous loci at the basis of convergent phenotypic evolution under domestication of leguminous species. Plant J. 2019;97(4):693–714.

    CAS  Article  PubMed  Google Scholar 

  12. 12.

    González AM, Yuste-Lisbona FJ, Saburido S, Bretones S, DeRon AM, Lozano R, et al. Major contribution of flowering time and vegetative growth to plant production in common bean as deduced from a comparative genetic mapping. Front PlantSci. 2016;7:1940.

    Article  Google Scholar 

  13. 13.

    Yuste-Lisbona FJ, González AM, Capel C, García-Alcázar M, Capel J, de Ron AM, et al. Genetic variation underlaying pod size and color traits of common bean depends on quantitative trait loci with epistatic effects. Mol Breed. 2014;33(4):939–52.

    Article  Google Scholar 

  14. 14.

    Hagerty CH, Cuesta-Marcos A, Cregan P, Song Q, McClean P, Myers JR. Mapping snap bean pod and color traits, in a dry bean × snap bean recombinant inbred population. J Amer Soc Hort Sci. 2016;141(2):131–8.

    CAS  Article  Google Scholar 

  15. 15.

    Murube E, Campa A, Song Q, McClean P, Ferreira JJ. Toward validation of QTL associated with pod and seed size in common bean using two nested recombinant inbred line populations. Mol Breeding. 2020;40(1):7.

    CAS  Article  Google Scholar 

  16. 16.

    Beninger CW, Hosfield GL, BAssett MJ, Owens S. Chemical and morphological expression of the B and asp seed coat genes in Phaseolus vulgaris. J Amer Soc Hort Sci. 2016;125(1):52–8.

    Article  Google Scholar 

  17. 17.

    Lamprecht H. Über die vererbung der roten vülsenfarbebei Phaseolus vulgaris. Agri Hort Genet. 1951;9:84–7.

    Google Scholar 

  18. 18.

    Yen DE. A shiny-podded mutant in pole bean (Phaseolus vulgaris L.). N Z J Sci Tech Sec A. 1957;38:820–4.

    Google Scholar 

  19. 19.

    Myers JR, Wallace LT, Mafi Moghaddam S, Kleintop AE, Echeverria D, Thompson HJ, et al. Improving the health benefits of snap bean: genome-wide association studies of total phenolic content. Nutrients. 2019;11(10):2509.

    CAS  Article  PubMed Central  Google Scholar 

  20. 20.

    Schmutz J, McClean PE, Mamidi S, Wu GA, Cannon SB, Grimwood J, et al. A reference genome for common bean and genome-wide analysis of dual domestications. Nature Genetic. 2014;46(7):707–13.

    CAS  Article  Google Scholar 

  21. 21.

    Burghardt LT, Young ND, Tiffin P. A guide to genome-wide association mapping in plants. Curr Protoc Plant Biol. 2017;2(1):22–38.

    Article  PubMed  Google Scholar 

  22. 22.

    Campa A, Murube E, Ferreira JJ. Genetic diversity, population structure, and linkage disequilibrium in a Spanish common bean diversity panel revealed through genotyping-by-sequencing. Genes. 2018;9(11):518.

    CAS  Article  PubMed Central  Google Scholar 

  23. 23.

    Jun X, Xin-yu W, Wang-zhen G. The cytochrome P450 superfamily: key players in plant development and defense. J Integr Agric. 2015;14(9):1673–86.

    CAS  Article  Google Scholar 

  24. 24.

    Ambawat S, Sharma P, Yadav NR, Yadav RC. MYB transcription factor genes as regulators for plant responses: an overview. Physiol Mol Biol Plants. 2013;9(3):307–21.

    CAS  Article  Google Scholar 

  25. 25.

    Gómez-Martín C, Capel C, González AM, Lebrón R, Yuste-Lisbona FJ, Hackenberg M, et al. Transcriptional Dynamics and Candidate Genes Involved in Pod Maturation of Common Bean (Phaseolus vulgaris L.). Plants. 2020;9(4):545.

    CAS  Article  PubMed Central  Google Scholar 

  26. 26.

    Wang X, Chen L, Wang A, Wang H, Tian J, et al. Quantitative trait loci analysis and genome-wide comparison for silique related traits in Brassica napus. BMC Plant Biol. 2016;16(1):71.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  27. 27.

    Marciniak K, Przedniczek K. Gibberellin Signaling Repressor LlDELLA1 Controls the flower and pod development of yellow lupine (Lupinus luteus L.). Int J Mol Sci. 2020;21(5):1815.

    CAS  Article  PubMed Central  Google Scholar 

  28. 28.

    Adams MW. Basis of yield component compensation in crop plants with special reference to the field bean, Phaseolus vulgaris. Crop Sci. 1967;7(5):505–10.

    Article  Google Scholar 

  29. 29.

    Safdar LB, Andleeb T, Latif S, Umer MJ, Tang M, Li X, et al. Genome-wide association study and QTL meta-analysis identified novel genomic loci controlling potassium use efficiency and agronomic traits in bread wheat. Front Plant Sci. 2020;11:70.

    Article  PubMed  PubMed Central  Google Scholar 

  30. 30.

    Zhang YM, Jia Z, Dunwell JM. Editorial: the applications of new multi-locus GWAS methodologies in the genetic dissection of complex traits. Front Plant Sci. 2019;11(10):100.

    Article  Google Scholar 

  31. 31.

    Gou JY, Felippes FF, Liu CJ, Weigel D, Wang JW. Negative regulation of anthocyanin biosynthesis in Arabidopsis by a miR156-targeted SPL transcription factor. Plant Cell. 2011;23(4):1512–22.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  32. 32.

    Falcone-Ferreyra ML, Rius SP, Casati P. Flavonoids: biosynthesis, biological functions, and biotechnological applications. Front Plant Sci. 2012;28(3):222.

    CAS  Article  Google Scholar 

  33. 33.

    Tanaka Y, Brugliera F. Flower colour and cytochromes P450. Philos Trans R Soc Lond Ser B Biol Sci. 2013;368(1612):20120432.

    CAS  Article  Google Scholar 

  34. 34.

    Kyle MM, Dickson MH. Linkage of hypersensitivity to five viruses with the B locus in Phaseolus vulgaris L. J Hered. 1988;79:308–11.

    Article  Google Scholar 

  35. 35.

    Vallejos CE, Astua-Monge G, Jones V, Plyler TR, Sakiyama NS, Mackenzie SA. Genetic and molecular characterization of the I locus of Phaseolus vulgaris. Genetics. 2006;172(2):1229–42.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  36. 36.

    Brodie A, Azaria JR, Ofran Y. How far from the SNP may the causative genes be?. 2016. Nucleic Acids Res. 2016;44(13):6046–54.

  37. 37.

    Liu J, Hua W, Hu Z, Yang H, Zhang L, Li R, Deng L, Sun X, Wang X, Wang H. Natural variation in ARF18 gene simultaneously affects seed weight and silique length in polyploid rapeseed. Proc Natl Acad Sci U S A. 2015;112(37):E5123–32.

  38. 38.

    Shi L, Song J, Guo C, Wang GZ, Yang P, Chen X, et al. A CACTA-like transposable element in the upstream region of BnaA9.CYP78A9 acts as an enhancer to increase silique length and seed weight in rapeseed. Plant J. 2019;98(3):524–39.

    CAS  Article  PubMed  Google Scholar 

  39. 39.

    Kaur J, Akhatar J, Goyal A, Kaur N, Kaur S, Mittal M, et al. Genome wide association mapping and candidate gene analysis for pod shatter resistance in Brassica juncea and its progenitor species. Mol Biol Rep. 2020;47(4):2963–74.

    CAS  Article  PubMed  Google Scholar 

  40. 40.

    Dubos C, Stracke R, Grotewold E, Weisshaar B, Martin C, Lepiniec L. MYB transcription factors in Arabidopsis. Trends Plant Sci. 2010;15(10):573–81.

    CAS  Article  PubMed  Google Scholar 

  41. 41.

    Vlasova A, Capella-Gutiérrez S, Rendón-Anaya M, Hernández-Oñate M, Minoche AE, Erb I, et al. Genome and transcriptome analysis of the Mesoamerican common bean and the role of gene duplications in establishing tissue and temporal specialization of genes. Genome Biol. 2016;17(1):32.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  42. 42.

    Brewer MT, Lang L, Fujimura K, Dujmovic N, Gray S, Knaap E. Development of a controlled vocabulary and software application to analyse fruit shape variation in tomato and other plant species. Plant Physiology. 2006;141:15–25.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  43. 43.

    Wickham H. ggplot2: Elegant graphics for data analysis. New York: Springer-Verlag; 2016. ISBN 978–3–319-24277-4,

    Book  Google Scholar 

  44. 44.

    Wei T, Simko V. R package "corrplot": Visualization of a correlation matrix (Version 0.84). 2017. Available from

    Google Scholar 

  45. 45.

    Kruijer W, Boer MP, Malosetti M, Flood PJ, Engel B, Kooke R, et al. Marker-based estimation of heritability in immortal populations. Genetics. 2015;199(2):379–98.

    Article  PubMed  Google Scholar 

  46. 46.

    R Core Team. The R project for statistical computing. 2019. Available at: www.R-projectorg/. Accessed: September 2019.

  47. 47.

    Elshire RJ, Glaubitz JC, Sun Q, Poland JA, Kawamoto K, Buckler ES, et al. A robust, simple genotyping-by-sequencing (GBS) approach for high diversity species. PLoS One. 2011;6(5):e19379

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  48. 48.

    Perea C, De La Hoz JF, Cruz DF, Lobaton JD, Izquierdo P, et al. Bioinformatic analysis of genotype by sequencing (GBS) data with NGSEP. BMC Genomics. 2016;17(S5):498.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  49. 49.

    Lobaton JD, Miller T, Gil J, Ariza D, de la Hoz JF, Soler A, et al. Resequencing of common bean identifies regions of inter–gene pool introgression and provides comprehensive resources for molecular breeding. Plant Genome. 2018;11(2):170068.

    CAS  Article  Google Scholar 

  50. 50.

    Bradbury PJ, Zhang Z, Kroon DE, Casstevens TM, Ramdoss Y, Buckler ES. TASSEL: software for association mapping of complex traits in diverse samples. Bioinformatics. 2007;23(19):2633–5.

    CAS  Article  PubMed  Google Scholar 

  51. 51.

    Lipka AE, Tian F, Wang Q, Peiffer J, Li M, Bradbury PJ, et al. GAPIT: genome association and prediction integrated tool. Bioinformatics. 2012;28(18):2397–9.

    CAS  Article  PubMed  Google Scholar 

  52. 52.

    Yu Y, Ouyang Y, Yao W. ShinyCircos: an R/shiny application for interactive creation of Circos plot. Bioinformatics. 2018;34(7):1229–31.

    CAS  Article  PubMed  Google Scholar 

Download references


This work was supported in part by grant AGL2017-87050-R of the AEI-Spanish Government. Authors thank M. Bueno, J.A. Poladura, and F. Díaz for their technical assistance.


This work was supported in part by grant AGL2017–87050-R of the State Research Agency (AEI)- Spanish Government.

Author information




CG, performed the phenotyping and in silico analysis. AC, performed the genotyping. AS, performed the genomic analyses. PM, contributed to writing the manuscript. JJF, conceived and prepared the manuscript, and conducted statistical data analysis. All authors read and approved the last version.

Corresponding author

Correspondence to Juan Jose Ferreira.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no conflict of interest.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

García-Fernández, C., Campa, A., Garzón, A.S. et al. GWAS of pod morphological and color characters in common bean. BMC Plant Biol 21, 184 (2021).

Download citation


  • Phaseolus vulgaris L.
  • Spanish diversity panel
  • Phenotyping
  • Association
  • Meta QTL candidate genes