ORIGINAL RESEARCH
published: 25 April 2022
doi: 10.3389/fpls.2022.830896
Improving Association Studies and
Genomic Predictions for Climbing
Beans With Data From Bush Bean
Populations
Edited by:
Eric Von Wettberg, Beat Keller 1†, Daniel Ariza-Suarez 1,2, Ana Elisabeth Portilla-Benavides 2,
University of Vermont, United States Hector Fabio Buendia 2, Johan Steven Aparicio 2, Winnyfred Amongi 3, Julius Mbiu 4,
Reviewed by: Susan Nchimbi Msolla 5, Phillip Miklas 6, Timothy G. Porch 7, James Burridge 8,
Clare Mukankusi 3Anupam Singh, , Bruno Studer 1* and Bodo Raatz 2†
University of Southern California,
1Molecular Plant Breeding, Institute of Agricultural Sciences, ETH Zurich, Zurich, Switzerland, 2Bean Program, International
United States
Center for Tropical Agriculture (CIAT), Cali, Colombia, 3 Bean Program, International Center for Tropical Agriculture (CIAT),
Lorenzo Raggi,
Kampala, Uganda, 4 Tanzania Agricultural Research Institute (TARI), Dodoma, Tanzania, 5Department of Crop Science and
University of Perugia, Italy
Horticulture, Sokoine University of Agriculture, Morogoro, Tanzania, 6Department of Agriculture, Agriculture Research Service
*Correspondence: (USDA-ARS), Prosser, WA, United States, 7Department of Agriculture, Agriculture Research Service (USDA-ARS), Tropical
Bruno Studer Agriculture Research Station, Mayaguez, PR, United States, 8Department of Plant Science, The Pennsylvania State
bruno.studer@usys.ethz.ch University, University Park, PA, United States
†Present addresses:
Beat Keller, Common bean (Phaseolus vulgaris L.) has two major origins of domestication, Andean
Crop Science Group, Institute of
and Mesoamerican, which contribute to the high diversity of growth type, pod and seed
Agricultural Sciences, ETH Zurich,
Zurich, Switzerland characteristics. The climbing growth habit is associated with increased days to flowering
Bodo Raatz, (DF), seed iron concentration (SdFe), nitrogen fixation, and yield. However, breeding
Limagrain Vegetable Seed,
La Menitŕe, France efforts in climbing beans have been limited and independent from bush type beans.
To advance climbing bean breeding, we carried out genome-wide association studies
Specialty section: and genomic predictions using 1,869 common bean lines belonging to five breeding
This article was submitted to
Plant Breeding, panels representing both gene pools and all growth types. The phenotypic data were
a section of the journal collected from 17 field trials and were complemented with 16 previously published trials.
Frontiers in Plant Science
Overall, 38 significant marker-trait associations were identified for growth habit, 14 for
Received: 07 December 2021 DF, 13 for 100 seed weight, three for SdFe, and one for yield. Except for DF, the results
Accepted: 25 February 2022
Published: 25 April 2022 suggest a common genetic basis for traits across all panels and growth types. Seven
Citation: QTL associated with growth habits were confirmed from earlier studies and four plausible
Keller B, Ariza-Suarez D, candidate genes for SdFe and 100 seed weight were newly identified. Furthermore, the
Portilla-Benavides AE, Buendia HF,
Aparicio JS, Amongi W, Mbiu J, genomic prediction accuracy for SdFe and yield in climbing beans improved up to 8.8%
Msolla SN, Miklas P, Porch TG, when bush-type bean lines were included in the training population. In conclusion, a large
Burridge J, Mukankusi C, Studer B population from different gene pools and growth types across multiple breeding panels
and Raatz B (2022) Improving
Association Studies and Genomic increased the power of genomic analyses and provides a solid and diverse germplasm
Predictions for Climbing Beans With base for genetic improvement of common bean.
Data From Bush Bean Populations.
Front. Plant Sci. 13:830896. Keywords: genome-wide association studies (GWAS), genomic selection, population structure, pleiotropy, growth
doi: 10.3389/fpls.2022.830896 habit, common bean (Phaseolus vulgaris L.), climbing and bush type bean
Frontiers in Plant Science | www.frontiersin.org 1 April 2022 | Volume 13 | Article 830896
Keller et al. Improving Genomic Predictions for Beans
1. INTRODUCTION in an increased yield of lines among and across the different
growth types (White and Izquierdo, 1989; Keller et al., 2020).
Common bean (Phaseolus vulgaris L.) was domesticated about However, within growth type I and type II, DPM was not
8,000 years ago in two geographic regions, resulting in the related to yield in near-isogenic lines (White et al., 1992). This
Andean and the Mesoamerican gene pools (Gepts et al., 1986; suggests that the relationship between yield, DF, and DPM can be
Bitocchi et al., 2013). Within the two gene pools, several partially uncoupled.
groups, including climbing and bush type beans, were identified The higher SdFe in climbing beans is promising to combat
through genetic or phenotypic characterization (Gepts et al., iron deficiency in human nutrition, which causes anemia,
1986; Rodriguez et al., 2016). Seven domestication events for the increases morbidity, and leads to economic losses (Boccio and
common bean were discovered by investigating a genetic locus Iyengar, 2003). About 30% of the global population suffers from
for flowering determinacy (Kwak et al., 2008, 2012). Flowering anemia, especially women and children in developing countries
determinacy defines the first criterion for growth type: the (Black et al., 2008; Stein, 2010). Increasing SdFe in legumes is
determinate growth type forms a reproductive terminal bud, a possible avenue to improve nutritional quality in the human
whereas the indeterminate growth types produce a vegetative diet (Petry et al., 2015; Rehman et al., 2019). In the last years, a
one (Singh, 1981). The second criterion describes the bush vs. few biofortified lines with higher SdFe were successfully released
climbing growth habit. Both criteria were used to distinguish four (HarvestPlus, 2022). Iron biofortified beans showed higher phytic
growth types: type I (determinate bush), type II (indeterminate acid concentrations, which decreased the relative but not absolute
bush), type III (indeterminate semi-climber), and type IV iron absorption (Petry et al., 2014). However, the SdFe of
(indeterminate climber) (Singh, 1981). Since growth type is climbing beans has not been investigated intensively and SdFe
associated with flowering determinacy, it also affects vegetative is negatively correlated with yield (Kelly and Bornowski, 2018).
growth and the length of the crop cycle (González et al., 2016). In Such tradeoffs, including the longer DPM, which is associated
recent decades, major breeding efforts have been directed toward with the climbing habit, need to be taken into account when
the erect growth habit of bush type beans since this habit enables improving climbing beans.
a faster cultivation cycle without staking of the plants and a single, Efficient breeding for multiple traits requires detailed
automated harvest (Teixeira et al., 1999; Ronner et al., 2018). knowledge about the underlying genetic architecture of the target
The joint diversity of common bean growth types may offer new traits and their correlations with other key characteristics. The
insights to improve not only climbing but also bush type beans. genetic architecture of traits can be investigated by genome-wide
Although largely neglected in breeding programs, climbing association studies (GWAS) and genomic prediction models
beans offer three main advantages over bush type beans: first, (Crossa et al., 2017; Cortes et al., 2021). For common bean,
climbing beans reach higher yield per area, with up to 5 t ha-1 GWAS have been carried out successfully and QTL for various
(Rosales-Serna et al., 2004; Checa et al., 2006; Barbosa et al., traits were tagged with molecular markers (Miklas et al., 2006;
2018). Second, they have a higher symbiotic nitrogen fixation Wu et al., 2020). Recently, genomic predictions were evaluated
capacity, with up to 92 kg of N fixed ha-1 (Graham, 1981; Bliss, for agronomic traits in an elite Andean breeding panel (Keller
1993; Checa et al., 2006; Barbosa et al., 2018). Third, they achieve et al., 2020). This study revealed the prediction abilities (PAs)
higher seed iron content (SdFe), with up to 10 mg/100 g (Blair for the genomic estimated breeding values (GEBV) based on
et al., 2010; Blair, 2013; Petry et al., 2015; Mukamuhirwa and genomic data only, as proposed by Meuwissen et al. (2001).
Rurangwa, 2018). Indeed, the production of climbing beans can This allows breeders to efficiently select superior lines before
be more profitable, and they are preferentially adopted in higher they enter field trials (Crossa et al., 2017). In general, PA is
and/or drought prone regions by small holder farmers in Uganda increased for more heritable traits and in lines, which are closely
and Rwanda (Ronner et al., 2018; Katungi et al., 2019). However, related to the training population (TP). In order to deal with
these advantages come with the cost of staking the plants and a different populations or breeding panels, multivariate models
longer vegetative period, partly due to the indeterminate growth were suggested to account for population structure (Lehermeier
type (White et al., 1992). et al., 2015). Following another approach, the TP can be
To shorten the cultivation cycle, Kornegay et al. (1992) optimized for genetic relatedness to improve the accuracy of the
suggested crossing type I and II bean lines with type IV lines GEBV (Akdemir and Isidro-Sánchez, 2019; Sarinelli et al., 2019).
to achieve an increase in yield while selecting against climbing These approaches have great potential to improve PAs, given the
ability. However, the climbing growth habit is tightly linked to increasing availability of genotypic and phenotypic data (Spindel
plant development and productivity. The relation of days to and McCouch, 2016). The efficient use of prediction models
flowering (DF), vegetative growth, and plant production was based on appropriate TPs and available data sets is a key factor
investigated in two Andean (type I) x Mesoamerican (type IV) in the development of new and adapted climbing bean cultivars.
recombinant inbred line populations (González et al., 2016). In this study, comprehensive genetic analyses across all bean
In those mixed climbing and bush type populations, the QTL growth types and the two gene pools were carried out using
containing the PvTFL1y flowering gene explained 32% of the 1,869 lines belonging to five breeding panels. We hypothesized
variation for DF, 66% for vegetative growth (length of the main that (i) climbing and bush type beans show, despite specific
stem), and 19% for the rate of plant production, including growth habit loci, overall a high genetic similarity, that (ii)
traits such as yield and seed weight (González et al., 2016). In combined analysis will increase the power of GWAS as well as
general, more days to physiological maturity (DPM) resulted the genomic predictions across all growth types and that (iii)
Frontiers in Plant Science | www.frontiersin.org 2 April 2022 | Volume 13 | Article 830896
Keller et al. Improving Genomic Predictions for Beans
predictions for climbing beans with an optimized TP including the days from planting until 50% of the plants in the plot had at
bush beans will outperform predictions which were based on the least one open flower. The seed yield per plot was normalized to a
climbing bean growth types only. Hence, we evaluated whether moisture content of 14% and extrapolated to yield per hectare.
breeding programs can supplement their trial data with already The weight of 100 seeds (100SdW) was measured separately.
existing data to improve GWAS and genomics predictions. The The growth type was assessed according to the four categories
overall objective was to provide molecular markers and genomic described by Singh (1981). The growth habit described climbing
prediction models which can be used to speed up the selection of ability and differentiated only between bush (types I and II) and
new bean lines across all growth types and gene pools. climbing types (type III and IV). Further traits such as DPM,
pod harvest index (PHI), SdFe, SdZn, and canning quality were
2. MATERIALS AND METHODS phenotyped only for the VEC. Analogous to DF, DPM represents
the days until 50% of the pods in one plot had lost their green
2.1. Germplasm pigmentation. For the PHI, the seed dry weight of 20 pods at
The germplasm used in this study consisted of five bean breeding harvest was divided by the corresponding pod dry weight. The
panels across all growth types: the newly composed climbing SdFe and SdZn were assessed on dried and ground seeds as
bean panel (VEC), the Andean diversity panel (ADP), the panel described by Stangoulis and Sison (2008). The SdFe and SdZn
of progeny from two-way crosses between five Andean and five content was then quantified by the X-ray fluorescence method
Mesoamerican parents (AxM), the Mesoamerican introgression using the X-Supreme 8000 instrument (Oxford Instruments,
panel (MIP), and the elite Andean breeding panel (VEF), UK) (Guild et al., 2017). The canning quality was assessed
whereas the latter four are bush type bean panels known from by a trained sensory panel at Michigan State University as
previous studies. described by Cichy et al. (2014). Briefly, the beans were soaked
in distilled water, canned at 100◦C, sealed, and stored for 2
2.1.1. Climbing Bean Panel weeks. Upon opening, the canning quality was assessed by a
The VEC comprised climbing bean lines of growth types III trained consumer panel and expressed as an overall score from 1
and IV. The lines were selected for grain quality, commercial (unacceptable appearance) to 5 (excellent appearance). The rated
seed type, disease resistance, SdFe, seed zinc concentration criteria included color, bean splitting, free starch clumps, and
(SdZn), and agronomic performance. They represent the genetic brine clarity after cooking (Cichy et al., 2014).
variation used in climbing bean breeding at the International
Center for Tropical Agriculture (CIAT). The VEC was composed 2.2.1. Field Trials for Climbing Beans
mainly of lines from the Andean gene pool, but it also included The field trials for the VEC were carried out at five
a line of the Mesoamerican gene pool (G2333) and a few lines locations in three countries: Darién (3◦53′31′′N 76◦31′00′′W,
of admixed origin. In total, the VEC was comprised of 290 lines altitude of 1,491 m a.s.l.), Palmira (3◦30′03.0′′N 76◦21′03.5′′W,
including twelve breeding groups, four genebank accessions, and altitude 965 m a.s.l.), and Popayán (2◦25′39′′N 76◦37′17′′W,
six cultivars (Supplementary Table 1). altitude of 1,750 m a.s.l.) in Colombia; Kawanda (0◦24′49′′N
Field trials including all VEC lines were conducted in 32◦31′59′′E, altitude of 1,190 m a.s.l.) in Uganda; and Kagera
Colombia, Uganda, and Tanzania to collect phenotypic data for (1◦24′56.5′′S 31◦46′48.8′′E, altitude of 1,320 m a.s.l.) in Tanzania
this study. (Supplementary Table 2). Each line was arranged in an alpha
lattice design with three replicates.
2.1.2. Bush Type Bean Panels Regarding the field trials in Colombia, the experimental units
Four bush type breeding panels were used in this study: (i) the consisted of one row with a row-to-row distance of 0.95 m
ADP consisting of 352 Andean bush cultivars and breeding lines and with seven seeds sown manually per meter row length.
from public and private breeding programs as described by Cichy Row length per plot differed between locations with 2.5 m,
et al. (2015). The ADP genotyping-by-sequencing (GBS) data was 2.2 m, and 2.0 m used in Darién, Palmira, and Popayán,
available via the ARS Feed-the-Future Grain Legumes Project respectively (Supplementary Table 2). Climbing beans (mainly
(arsftfbean.uprm.edu/bean/?p=472); (ii) the AxM as described by type IV) required a trellis to support the plant. Wooden poles of 3
Mayor Duran et al. (2016) and Mayor Duran (2016); (iii) the m height were distributed in the field in squares of 5 x 5 m. Wires
MIP consisting ofMesoamerican breeding lines with interspecific were spanned between the poles at a height of 2.3 m. Each bean
introgressions from P. acutifolious, P. dumosus, and P. coccineus plant climbed on a string hanging from the wire. Eight strings
in their pedigrees and their parental lines as described by Diaz per plot were deployed. Plant protection was carried out when
et al. (2021); (iv) the VEF as described by Keller et al. (2020). needed using good agricultural practices.
For the VEF, MIP, and AxM, phenotypic and genotypic data of The soil type in Darién and Popayán was an Inceptisol (Typic
605, 217 and 200 lines were available, respectively. For the ADP, Dystrandept) with about 70 g/kg and 140 g/kg of soil organic
field trials were conducted in Mozambique, Tanzania, and in the matter, respectively, and a pH of around 5 (Barbosa et al., 2018).
United States to collect phenotypic data for this study. The soil type in Palmira was Mollisol with a pH between 7.0
and 7.5.
2.2. Phenotyping In Uganda and Tanzania, the plots were planted in single 3.0 m
Agronomic traits were evaluated in the VEC and ADP as rows with about 30 seeds planted per row and at 0.6 m between
previously described by Keller et al. (2020). Briefly, DF represents the rows (Supplementary Table 2). Granular N:P:K (nitrogen,
Frontiers in Plant Science | www.frontiersin.org 3 April 2022 | Volume 13 | Article 830896
Keller et al. Improving Genomic Predictions for Beans
P205 andK2O) fertilizer wasmanually applied at 125 kg/ha. Trials with 55 additional climbing lines from a previous study (Barbosa
were manually weeded two times. Climbing bean plants were et al., 2018).
staked with 2 to 3 m long wooden poles. Five poles were used
per plot. 2.4. Data Analysis
2.4.1. Phenotypic Data
2.2.2. Field Trials for Bush Beans Best linear unbiased estimators (BLUEs) were extracted from
The field trials for the ADP were carried out at five locations phenotypic data in two stages. In the first stage, the field data were
in three countries: Arusha (3◦21′41′′S 36◦37′34′′E, altitude of corrected for spatial effects using the SpATS R package, setting
1,387 m a.s.l.), Morogoro (6◦51′14.2′′S 37◦39′27.6′′E, altitude row and column as random effects (Rodríguez-Álvarez et al.,
of 526 m a.s.l.), and Mbeya (8◦54′52′′S 33◦31′05′′E, altitude of 2018). The number of plants harvested was binned (binwidth
1,780 m a.s.l.) in Tanzania; Chokwe (24◦30′04′′S 33◦00′09′′E, = 5) and added as a random effect in the spatial analysis. The
altitude of 35 m a.s.l.) in Mozambique; and Wilcox (32◦01′44′′N plots with residuals bigger than ±3 times the SD were treated
109◦41′27′′W, altitude of 1,321 m a.s.l.) in the United States as outliers and removed iteratively as described in Keller et al.
(Supplementary Table 3). (2020). The BLUEs were extracted for each line in each trial
The trials were conducted in two to three replications in from the SpATS model (first-stage BLUEs) from all the data sets,
a randomized complete block design (Supplementary Table 3). except for the ADP. In the ADP lines, the first-stage BLUEs were
The plot size differed between trials with lengths of 5–8 m, extracted using replicate (block) as a factor for the fixed effects
1–0.5 m spacing between rows, and 1–4 rows. The plot sizes, since the row and column information was not available for these
soil types, and growth conditions of each trial are described in randomized complete block design trials. In the second stage,
Supplementary Table 3. Irrigated trials were watered about two second-stage BLUEs for each line (Li) were calculated across all
times a month and drought trials were rain fed. The growing trials using the following model:
season of one trial (MZCH15D_heat) coincided with the high
temperature season at the site. The trials in Tanzania and in the yij = Li + Ej + εij (1)
United States were not fertilized. In the trials in Mozambique,
single superphosphate and ammonium sulfate were applied at where yij is the first-stage BLUE of the ith line in the jth trial,
about 30 kg P/ha and 6.3 kg N/ha, respectively. Li is the fixed effect for each line, Ej is the fixed effect for each
trial, and εij the error term. The inverse of the squared standard
2.2.3. Available Phenotypic Data From Previous error⊕(SE) of the mean was used as a weight, i.e., ε ∼ N(0,R);
R n
= j=0 Rj where RStudies j = diag[(SEij)2], and (SEij) is the standard
Of the 290 VEC lines, 43 were phenotypically characterized error of a mean of the ith line in the jth trial (Möhring and
in two trials in a previous study (Barbosa et al., 2018). Piepho, 2009). In addition, BLUEs for yield were scaled (mean
Phenotypic data were also publicly available for the AxM = 0, SD=1 resulting in Yd_scaled) for each panel to compare
(Mayor Duran, 2016), the MIP (Diaz et al., 2021), and the VEF relative differences between the lines beyond panel and growth
(Keller et al., 2020) from ten, one and three trials, respectively type. Kernel density estimates of the BLUEs were calculated using
(Supplementary Table 4). a Gaussian kernel with 1/30 bandwidth of the data range. The
estimates were drawn as a smoothed histogram with the integral
equal to one using the ggplot R package (Wickham, 2016). The
2.3. Genotyping significance of the genotype-by-environment interaction (GxE)
GBS was conducted as previously described by Nay et al. was tested by comparingmodel 1 with and without an interaction
(2019). Briefly, DNA was extracted from leaf tissue and digested term ELij between jth trial and ith line using the likelihood-ratio
with the ApeKI restriction enzyme as described in Elshire test. The test statistic was compared with a chi-square value with
et al. (2011). The sequencing was performed on the Illumina one degree of freedom and the p-value was adjusted following the
HiSeq 2500 platform at the HudsonAlpha Genome Sequencing work by Self and Liang (1987).
Center (Huntsville, AL, United States). The sequence reads were
demultiplexed using NGSEP (v3.3.0) (Tello et al., 2019), trimmed 2.4.2. Linkage Disequilibrium and Population
using Trimmomatic (v0.36) (Bolger et al., 2014), and aligned Structure
to the reference genome of P. vulgaris G19833 v2.1 (Schmutz Pairwise measures of linkage disequilibrium (LD) were calculated
et al., 2014) using Bowtie2 (v2.2.30) (Langmead and Salzberg, for each population in sliding windows of 100 markers. The LD
2012). The variant calling was carried out using NGSEP, filtering measures were corrected for kinship in the population (r2V ) as
SNPs with a genotype quality below 40, minor allele frequency implemented in the R package LDcorSV (v1.3.2) (Mangin et al.,
(MAF) below 0.05, and removing SNPs with less than 60% of 2012). The LD decay was estimated regressing the pairwise r2V
genotype calls, which yielded a matrix with 20% of missing data. values on the physical distance of their markers using the locally
The imputation of the missing data was performed with Beagle estimated scatterplot smoothing implemented in the R function
v.5.0 (Browning et al., 2018) using 100 as an effective population “loess” (v4.1.0), with a span value of 0.5.
size and using the genetic map reported by Diaz et al. (2019). The The phylogenetic tree was constructed using the GGTREE R
SNP calling was carried out once on the VEC separately and once package on the hierarchically clustered SNP matrix (Yu et al.,
on all five panels together. All VEC lines were genotyped together 2017). The population structure was assessed on the SNP matrix
Frontiers in Plant Science | www.frontiersin.org 4 April 2022 | Volume 13 | Article 830896
Keller et al. Improving Genomic Predictions for Beans
using principal components analysis (PCA) implemented in the second-stage BLUEs). In the case of the genotype model using all
FactoMineR R package (Lê et al., 2008). The correlation of the trials separately, j represents always the same environment. The
supplemental phenotypic traits (second-stage BLUEs) with the same applied to the genotype model when using second-stage
principal components (PCs) of the SNP matrix was calculated BLUEs across all trials.
using the same FactoMineR package (Lê et al., 2008).
2.5.2. GxE Model
2.4.3. Genome-Wide Association Studies In the GxE model, the GEBV for each trait were estimated for
To carry out GWAS, the Bayesian-information and Linkage- each environment (i.e., location) by adding an effect for the
disequilibrium Iteratively Nested Keyway (BLINK) algorithm interaction between the jth environment and ith GEBV (gEij) to
implemented in GAPIT was used (Huang et al., 2019). The the model (2). This resulted in the following model:
first five PCs were used to correct for population structure. The
imputed SNP matrix was used for GWAS. Identified QTL tagged yij = gi + Ej + gEij + εij (3)
by SNP markers were labeled with “Trait_Chr_Postion” as QTL
ID, whereas yield was abbreviated with Yd and growth habit with gE ∼ N(0, I ⊗ K 2
σgE), where I is the identity matrix for the
with GH. The SNP position derived from the reference genome environments and⊗ denotes the Kronecker product. This means
G19833 v2.1 (Schmutz et al., 2014) was in Mbp rounded to no correlations between environments were considered.
two digits.
The network of significant marker-trait associations was 2.5.3. Factor Analysis Model
visualized similarly as suggested in Fang et al. (2017) using In the factor analysis (FA) model, the GEBV for each line in each
the ggnetwork R package (Briatte et al., 2020). The distance environment (gij) were estimated for each trait using SNPmarker
between the connected nodes represents the LD between the two information and the covariance of the phenotypes (yij) between
SNPs calculated as the squared Pearson correlation coefficient. the trials. The following equation was used:
The SNPs were connected when LD >0.25. The haplotypes
associated with one trait were assembled using the identified yij = gij + εij (4)
SNPs in the indicated region for that trait. In case there were
with g ∼ Nj(0,G⊗K 2
σg ), where G represents a covariance matrix
less than six SNPs selected, these SNPs were directly assembled
to haplotypes. When more than five SNPs were selected, groups of phenotypes between trials calculated asG = BBT+9 , where B
of haplotypes were constructed by hierarchical clustering of all is a matrix of loadings (regressions of the original random effects
lines based on those SNPs using the stats R package. The optimal into common factors) and9 is a diagonal matrix whose non-null
number of clusters was determined using the average silhouette entries give the variances of factors that are trait-specific. The
width implemented in the factoextra R package (Kassambara and model residuals were assumed to follow a multivariate normal
Mundt, 2020). distribution ε ∼ Nj(0,Rε ⊗ In), where Rε is a covariance matrix
of model residuals and In represents an n-dimensional identity
2.5. Genomic Prediction matrix, where n is the number of phenotypes per environment.
The GEBV were estimated using Bayesian generalized linear Three common factors were selected.
regression (BGLR) implemented in the BGLR R package (Pérez
and de los Campos, 2014). For the factorial models, the BGLR 2.5.4. Training Population Optimization
extension for the Multiple-Trait Model (MTM) was used (http:// For each VEC line in the validation set, the 50 closest related
quantgen.github.io/MTM/vignette.html). The Gibbs sampler ran lines from all five panels were added to the TP (TP optimized).
with 20,000 iterations of which the first 10,000 were burned- Genetic relationships were calculated as the cophenetic distances
in and the remaining were thinned by factor 5. Three different of the hierarchically clustered lines based on the SNP matrix,
model approaches were tested. i.e., as the height of the dendrogram where the two branches
including the considered lines join each other. This means that
2.5.1. Genotype Model Among and Across Trials the TP could consist of between 50 lines (in case the lines to be
For each trait, the phenotypes adjusted per trial (yij) were tested would have all the same 50 closest relatives) and 2,500 lines
modeled as the sum of the GEBV for each line (gi) estimated (if they would not have the same closest relatives). The number
based on SNP marker information, i.e., the additive relationship of related lines (50) was chosen arbitrarily but seemed reasonable
matrix (K), a fixed effect for the trial (Ej), and an error term for a panel of 1,500 lines with some population structure. The
εij ∼ N(0, 2
σ ). The following linear model was used: optimized TP was compared to a TP containing only VEC lines
ε
(TP VEC) and a TP containing all lines from all five panels
yij = gi + Ej + εij (2) (TP extended).
with g ∼ N(0,K 2
σg ) and K was calculated as the normalized cross 2.5.5. Cross-Validation
product of the SNPmatrix using the rrBLUP package (Endelman, The cross-validation was done 100 times, randomly splitting
2011; Endelman and Jannink, 2012). The yij were calculated the dataset into training and validation set, i.e., for each cross-
either among all trials (using first-stage BLUEs), for all trials validation step, 70% of the lines were selected for the training
separately (using first-stage BLUEs), or across all trials (using and 30% for the validation set. The lines in the training and
Frontiers in Plant Science | www.frontiersin.org 5 April 2022 | Volume 13 | Article 830896
Keller et al. Improving Genomic Predictions for Beans
validation sets were referred to as TP and new lines, respectively. 3.2. Comparing Phenotypes of the Five
The Pearson correlation coefficient (r) between predicted and Breeding Panels
observed values was calculated at each cross-validation step to The phenotypic distribution of DF, 100SdW, SdFe, and yield
assess the PA. The prediction accuracy (PAcc) is defined as the were compared among all five panels across all trials (Figure 2A).
quotient of PA and the square root of heritability. The climbing beans of the VEC showed on average 28%, 21%,
and 67% higher DF, SdFe, and yield than the bush type panels,
2.5.6. Genomic BLUPs Without Cross-Validation respectively. The trait 100SdW depended primarily on the gene
The genomic r is defined as the Pearson correlation coefficient of pools, showing lower values in the Mesoamerican MIP and the
modeled vs. observed values when all lines were used in the TP AxM, consisting of inter gene pool crosses. The observation in
to fit the model. The genomic r was derived from the predicted the VEC, that SdFe correlated negatively and DF and 100SdW
genomic BLUPs without cross-validation. positively to yield, was also true in the four bush bean panels
(Supplementary Figure 4).
3. RESULTS 3.3. Structure and Diversity of the Five
Breeding Panels
The traits DF, 100SdW, SdFe, and yield were analyzed in In the joint analysis of all five panels, 14,913 SNP markers
27, 28, 6, and 33 field trials, respectively, using a total distributed over the whole genome were kept from the raw
of 1,869 lines belonging to five different breeding panels 169,087 SNPs after filtering for genotype quality calls, MAF, and
(Supplementary Figure 1). For this study, six and eleven field missingness. The dendrogram of the 1,869 lines showed grouping
trials were newly evaluated, adding data for the VEC and ADP, intoMesoamerican (represented by theMIP) and Andean origin,
respectively. The number of evaluated lines per trial was 290 for whereas the AxM and part of the VEC formed an admixture
the VEC (Supplementary Table 2) and ranged between 41 and branch (Figure 2B). These genetic groups were also visible in the
268 for the ADP (Supplementary Table 3). PCA analysis of all lines: the first PC clearly grouped the lines
according to their Andean and Mesoamerican origin spreading
3.1. Phenotypes of Climbing Beans the admixed AxM lines in-between (Figure 2C). In agreement,
The lines of the VEC showed different phenotypic distributions the first PC was highly correlated with 100SdW (r = −0.67)
for the traits DF, 100SdW, SdFe, and yield among all field which differentiated the two gene pools (Figures 2A,C). The
trials (Figure 1A). Especially, strong environmental effects were second PC explained variation for growth type, separating mostly
observed for yield. The trials carried out in the lowlands the VEC lines from the others, and was correlated to DF (r = 0.58,
and in the warmer climates (Palmira, Kawanda, and Kagera) Figure 2C). The first and second PC explained 28.9 and 3.4%
had less yield compared to the remaining trials in the high- of the genetic variance, respectively. The third PC, explaining
altitude locations Darién and Popayán. Phenotypic variation 2.4% of the variance, captured variationmainly between the AxM
among and across trials was also observed for additional and MIP, whereas the fourth and fifth PC showed the smallest
traits, including important breeding targets such as canning variation for the VEF (Supplementary Figure 5). Finally, the
quality and PHI (Supplementary Figure 2). The phenotypic sixth PC showed no clear pattern among the panels. Therefore,
correlations between trials were evaluated using the first-stage five first PCs were included as fixed effects for GWAS. The LD
BLUEs: positive correlations were revealed across the trials decay observed for all the breeding panels was faster for the
for DF, 100SdW, and SdFe, whereas yield of the Pal19D and combined panel than for the separate panels, enabling higher
TzKg19D trials were mainly negatively correlated to the other detection power for GWAS (Figure 2D). In general, the genetic
trials (Figure 1B). The likelihood-ratio test confirmed significant diversity was bigger between the two gene pools than between the
GxE (p value < 0.001) for all traits, especially for yield, where growth types.
the variance component for GxE was higher than for the lines
(Supplementary Table 5). 3.4. Genome-Wide Association Studies
Across all trials, the correlations between traits were Across All Panels
evaluated using the second-stage BLUEs. Positive correlations Carrying out GWAS using the second-stage BLUEs across all
were revealed between the traits SdFe, SdZn, DPM, and trials and breeding panels, a total of 69 significant marker-trait
DF (Supplementary Figure 3). However, the correlations were associations were identified below the 1% Bonferroni corrected
negative between yield, nitrogen use efficiency, and 100SdW. significance level (Figure 3A and Table 1). The observed p-
Additionally, PHI was negatively correlated to SdFe and SdZn. value distribution showed a clear deviation of the identified
Since yield showed strong GxE, the correlations to other traits significant SNPs from the expected uniform distribution if no
differed across trials, e.g., the yield was negatively correlated with genetic linkage was present (Figure 3B). A p-value inflation
DF in the Pop15B and Pal19D trials but positively correlated was visible mainly in growth habit and 100 SdW. Most striking
with DF in the remaining trials (Supplementary Figure 4). was the region at around 45Mbp on Chr 1 associated with
In summary, strong GxE was observed for yield, while the significant effects for growth habit, DF, DPM, and 100 SdW
remaining traits showed moderate GxE among the different (Supplementary Figure 6). Furthermore, different QTL for SdFe
environments and trials. and scaled yield across all panels were identified.
Frontiers in Plant Science | www.frontiersin.org 6 April 2022 | Volume 13 | Article 830896
Keller et al. Improving Genomic Predictions for Beans
A Dar14B Dar19B Pop15B TzKg19D
Dar18B Pal19D Pop17D UgKw19D
0.0015
0.075
0.2 0.04
0.0010
0.050
0.1 0.02
0.0005
0.025
0.0 0.000 0.00 0.0000
40 50 60 70 0 25 50 75 50 75 100 0 2000 4000 6000
DF [d] 100SdW [g] SdFe [mg/kg] Yield [kg/ha]
BLUEs
B
Pearson
correlation
−1.0 −0.5 0.0 0.5 1.0
DF 100SdW SdFe Yield
Pal19D Pop17D Dar19B
Dar19B
Dar18B Dar14B Dar18B
Pop17D Dar14B
Dar19B Dar18B Pop17D
Dar19B Dar18B UgKw19D
UgKw19D Pop15B Dar14B Pop15B
Pop15B Pal19D Pal19D
Dar14B TzKg19D Pop17D TzKg19D
FIGURE 1 | Phenotypes of the climbing bean panel (VEC). (A) Density diagrams for days to flowering (DF), 100 seed weight (100SdW), seed iron concentration
(SdFe), and seed yield of up to 290 VEC lines in eight trials are shown. Best linear unbiased predictors (BLUEs) were calculated from each trial corrected for spatial
effects in the field. (B) Pearson correlation coefficients were calculated across all trials for each trait based on the BLUEs. Trials were abbreviated based on the
location Darién (Dar), Palmira (Pal), Popayán (Pop) in Colombia, Kagera in Tanzania (TzKg), or Kawanda in Uganda (UgKw), the year and the planting season
(sequentially A to D). For a detailed description of each trial refer to Supplementary Table 2.
3.4.1. Growth Habit and Pleiotropy contrast, two QTL on Chr 4 (GH_4_1.42 and GH_4_40.05)
For growth habit, 38 significant marker-trait associations and on QTL on Chr 5 (GH_5_0.74) differentiated bush
were detected (Table 1). These QTL determined growth from climbing types. The minor SNP variant on Chr 2
habit, e.g., bush (type I and II) from climbing beans (type (GH_2_40.05) was the only one exclusively associated with the
III and IV) or they distinguished growth determinacy, e.g., two climbing growth types (type III and IV). However, this
separated determinate types I from indeterminate type II, III, association is to interpret cautiously since this SNP variant was
and IV (Supplementary Figure 7A). The SNPs associated rare (MAF= 0.05).
with growth habits on Chr 1, 3, and 6 (GH_1_43.71, The region at the end of Chr 1 showed significant pleiotropic
GH_3_1.29, and GH_6_23.87) differentiated mainly the effects on different traits (Supplementary Figure 6). In that
determinate growth type I from the other three types. In region, significant SNPs were identified not only for growth
Frontiers in Plant Science | www.frontiersin.org 7 April 2022 | Volume 13 | Article 830896
Density
Dar14B
Pop15B
UgKw19D
Dar19B
Pop17D
Dar18B
Pal19D
TzKg19D
Pal19D
Pop15B
Dar18B
Dar19B
Dar14B
Pop17D
Pop17D
Dar14B
Dar18B
Dar19B
TzKg19D
Pal19D
Pop15B
UgKw19D
Pop17D
Dar14B
Dar18B
Dar19B
Keller et al. Improving Genomic Predictions for Beans
A ADP AxM MIP VEC VEF
0.100 0.06
0.3 0.0020
0.075
0.2 0.04 0.0015
0.050
0.02 0.0010
0.1 0.025 0.0005
0.0 0.000 0.00 0.0000
30 40 50 60 20 40 60 20 40 60 80 100 0 1000 2000 3000
DF [d] 100SdW [g] SdFe [mg/kg] Yield [kg/ha]
B
M
M CA
R IRC 7
R C __
A _ 0042 _1
6
D 0 45_ 0 M 1
E 0 8 09
N 5 AX 0 422
F 1 P−_0
E _ _00
N 0 M 2
F 2 AD
_ 8 AX _02
CC 2 M 57
E GG 2
N AA 6 AX _0032
_ 527
F __
_ 00 M _002
17
86
M 21 03 AX
XM__0
007
N 1 A XM _ 7
C A M 18
_ _
1 AX
EN 26 AXM
7
M F
A _ 13
A 0
C _ M_ 32
_1 7
MEE
5
GNN 0 9
0 AX 1442
N A 9 XM 7
CF
F A M__110
__
_ 001 AX _
3 350 AXM
1 46 AXM
0
AD
D P 012
−
M A 0
A A 0 M_
26 AX
MN C _
D
AD _ 1
AZZ 1 5
A _072
A __ 1 5
_010 M
DP 310
2
29 AX
BN −07 127
D A 9
A _0 1 AXM_
B 0 _102
_ 7 12127
AB 6 AXM__1
U 1
_ 6
0 AXM
CAG 3 XM_112
B 7 A
UA__0064
_037
10 AXM 8729
AXDMP−_013
61
ADP−04
G1759
DDAABB__ G 40GR_0 713
N 56 ADP−0
U 8059
A_326 G233 65
NMU _2
AC_3148
MN 4
C_468 MNC__45923
CMB_073 AXM_062
AXM_152
MMACCA__004477
DDAA MNC_359
BB__99
A 32
D 54
P−0587
PAN_127 ENF_
ENF_0
_19
012
820
ADDAPN−_0016163 ENF_085
NCC_092 ABU_044
ABDNPA−_0060525
ADP−0603
ADP−0096 D
ADP−0671 DA
AA
B_
_13
2 0
31
DAN_007
DAB_633 DDAAAB__001602
RMA__0072
01150
A DDAB_
ACC_ DA 5
A A 41
A _
_ 1
035 0 0
NUA_ D 9 1
027 A 6
B_9
ABU_
32 D 2
A 9
A_07
ABU_0 3
029
DAN_
D
09 AS
A
DE
B
_022 PQ
_
D −_
5
ABU_0
01
7
A 10
4
B_5 127 97
AA 9
5 DD
−057
A PP−−
6 D 00
ADP 50
A P 90
D − 21
P 0
_001−1
− 04
0 6
63
SAP 0
14 D
U__0
0
9 A A
AB 01
A_ 9 C D BA
54 O P __
− 31
BN 1 D
−0 S 0 93
55 _ 5 85A
A_ 11 CA A_ 0 6
ADP 97
LP _0
_1 AB 1 1 1
U_ 2 6
L_ 0
ALXP
MA 3
14 1 LL 1 4
PP 44
A_ 05 AA 39
__
DA N_
2 D 74L 60
DA 1190 P
A__5 528636 AD
A 06A
CA
N_
DA _0
_0 D CN
_407
DAB
AN_ A _ 11
DAA A 0 2
DD
0 AA _ 2
BB
0
0 22
56 6 AA __ 0
B_ 59
7 DD 99 5
14
SA B_
62 A PP 90
DA 22 D −−−
B_
_6 A P 000055
DA 8684 F − 324
DAB
001
0
12 R 0 680
_ 3
P− 7 9
P−0
P−0 0 5
8
AD
AAD
C 150 DF D
0.4
AxM
ENF_098
100 0.3 ADP
Yield MIP
50 0.2
VEF
DOR_500 0.1
0 VEC
DAB_606 SdW
SEC_035 Combined
0.0
−50 0 50 100 150 0.0 0.5 1.0 1.5 2.0
PC1: 28.9% Distance (Mbp)
FIGURE 2 | Phenotypic and genotypic characterization of five common bean breeding panels including lines with bush and climbing growth habit originating from the
Andean and Mesoamerican gene pools. (A) Density diagrams of best linear unbiased estimators among five breeding bean panels are shown for days to
(Continued)
Frontiers in Plant Science | www.frontiersin.org 8 April 2022 | Volume 13 | Article 830896
Density
PC 2:  3.4 %
A
AA DP
AA D
D −
DD P
P 0
PP −
−0 52
−−− 0 0 3000 0 7000 6587 1
1
116
D
SD A
R MA
A
AD A C
B _
P A _
_ 1
11 0
− _ 20
7
0 0 90 2 4AD 6 1P 6
AA −D 0D
A P
3
P−−
5
0 4
AA D 01A DD P 40
A D PP−− 672
A D P 0
0
A 6
DDP −
− 664
P 00 60
P 66
AD −−− 6
5
A DA 000 0
0
D P 6
66
PA− 7
48
−_0 6
57
006D 28
1
A 13
0
B 3_56
A 9
DP−01
A 8D 3P−05
A 3D 4
A PD −P 0
A −
0
D 0
5
P 0 6
− 1
0 6
011
ADP−003
AA 1
DD
PP
−−
00602006
AADDPP−−00224721
ADP−0021
ADP−0
A 4
D 80
AA P
DD −
PP 0
−− 2
00
0
53
6
6093
MCGBCA__005096
EENNFF__010281
ENF_001
SMC_199
SMC_213
SEC
BS
_
NE
07
C 1
A__00S 1726
MC_017
SSMMCC__211597
SSMMRR__
M 0074
SS
S 2
6
RR_
SS
MM __16
MM
R 11
NN 5
385
0 8__0
2 6750
LD ( S
r EN_116
V) SSCERC
S _
_
0044E _ 5
0
SE FF 4
S _ 0
0 7
E 5F 5SE _0
S FE _
68
SE F
0
_ 31F_ 00 41 75
ALB_034
MSSMIEBCC___702882066S
S M
MSEM C
IBNN
_
___
0
0 3
F 4122
7
E 9937B
DSD
BTG
_
AEO 2A ARR 2L ___ 6S A
5
M T
1A 013
B _34M914
16
284
_ 7A0 7NA 0 C
S _ 6 A
I 6NE 8
A BA 6X __M 800
_ 4 5
1 1
72
34 0
_0 4
2
520490
DAN
P− 18
UA_
D _
0 3
N _
0022536A
DAB
DANN_
_0
DA −0
4
36
BNA
ADP
037
B_
0
_602
CM C_
DA
B
AC 3
5
−0
6
AD
P
1
04
−0 851
−0 15
0
AD
P 5 0
AD
P
P−
004
P− 1
ADAD 0
3
A_
DA 5
2
_0
DA
A
040
N_
DA _
004
DA
Z 057
A_
DA _
021
DA
A
A_0
68
DA _
062
DA
A __15
4184
DDAAAB
_01
7
SAA
04
NDZ
_0 81
− 5
ADP
0
09
DAB
_9064579DSAAAB
__ 0502
5
ACC_DAB_
2
DAB_
366
DDAAAB
__0289958
DAA_0
7
DAB_914
DAB_904
LAPBAU__604536
ADP−0682
CMB_012
ABADRDPBP−_−10095485360
AXM_067
AXM_082AXM_087
AXM_092
DAN_046
AAXXM_077M_177
AXM_04
7
A _1
62
ADXXA
MMA__01
1567
MVAB_
091
AN_04
1421
NUC_
70
ENF
_0
23
ENF
_0
A__
4021
092
MCNGC
ENF
_20
90
EN
F_1
188
C_ 015103
NU C
_
C_0
07
NCNCC
__03
5
NCNC
C
78
_1 69
MN
C _1 817
EN
FF_
EN 73
_3
NU
C
70
_0 19
NC 8M _
01136
VC
M
BC
__1
M 9
5
EN
F
_0
2
NUC
03 33
_1
3
07979
BC_
60
MNCAC_
_0_0
M _1
7
4
AC
M AI
CR 5
2
MC _
09
_1
IIRC
Keller et al. Improving Genomic Predictions for Beans
FIGURE 2 | flowering (DF), 100 seed weight (100SdW), seed iron concentration (SdFe), and seed yield. (B) Dendrogram of 1,869 lines characterized by 14,913 SNPs
shows the hierarchical relationships between lines and panels [following the same color code for panels as in (A)]. (C) Principal components (PC) 1 and 2 visualize the
genetic similarity across all five breeding panels. The arrows show quantitative supplementary phenotypic traits. Their cosines indicate the correlation with PC axes
and their length approximate the SD of the variable. The extreme lines on the PC 1 axis are labeled. (D) Linkage disequilibrium (LD) decay is shown for all panels
separately and combined. The LD was calculated in sliding windows of 100 markers and corrected for kinship in the population (r2V ).
A B
30 30
20 20
10 10
0 0
30 30
20 20
10 10
0 0
12 12
8 8
4 4
0 0
10.0 10.0
7.5 7.5
5.0
5.0
2.5
2.5
0.0
0.0
6
6
4
4
2
2
0
0 0 1 2 3 4
Chr1 Chr2 Chr3 Chr4 Chr5 Chr6 Chr7 Chr8 Chr9 Chr10 Chr11 −log10(expected p)
SNP marker
FIGURE 3 | Genome-wide association studies among five breeding panels differing in growth habits. (A) The Manhattan plots show the genetic associations with
climbing growth habit, days to flowering (DF), 100 seed weight (100SdW), seed iron concentration (SdFe), and seed yield. Seed yield was scaled among panels to
allow comparison between them. The horizontal black lines show the Bonferroni corrected significance threshold at the 1% level. The vertical lines indicate the
position of the two most significant markers for each trait. (B) Quantile distribution plots show the deviation of expected to observed p-values of SNP to trait
associations for each trait.
habits but also for DF, DPM, and 100SdW. Additionally, the on Chr 1 (Supplementary Figure 7B). Interestingly, the QTL
significant SNPs for growth habit on Chr 3 and 6 showed a on Chr 4 (GH_4_1.42 and GH_4_40.05) showed again a
tendency toward pleiotropic effects as observed for the QTL different pattern than the others: these SNPs increased 100SdW
Frontiers in Plant Science | www.frontiersin.org 9 April 2022 | Volume 13 | Article 830896
−log10(p)
Yield_scaled SdFe 100SdW DF Growth habit
−log10(observed p)
Keller et al. Improving Genomic Predictions for Beans
with increasing DF while, surprisingly, yield (scaled across 3.5. GWAS Within the Climbing Bean
populations) decreased. In summary, the 38 SNPs significantly Germplasm
associated with growth habit determined the four different For the genetic analyses within the climbing bean
growth types in different proportions, whereas only a few of these germplasm, a total of 15,589 SNPs were identified in
SNPs expressed pleiotropic effects. the VEC. The population structure was moderate with
PC1 and PC2 explaining 19.1 and 5.8 % of the genetic
variance, respectively (Supplementary Figure 9). In total,
3.4.2. SNP Effects Among Breeding Panels 22 significant marker-trait associations were identified in
In general, the significant marker-trait associations identified the VEC (Supplementary Table 6). Two QTL for PHI
in the whole population of 1,869 lines showed effects also in and DF were identified on Chr 5 in proximity at 38.67
the panels separately (Figure 4). An exception was the QTL and 39.34 Mbp, respectively, indicating tight linkage
DF_9_29.38 which revealed significant associations only in the (Supplementary Figure 10A). Several SNPs significantly
VEF and ADP, the two Andean panels. An interesting breeding associated with SdFe and SdZn were identified on Chr 2, 4, 6, 7,
target for yield is Yd_7_4.86, expressing a clear effect in all five and 10. Furthermore, a QTL for canning quality was identified
panels. This QTL was physically linked to DF_7_4.82 and in LD on Chr 7 at 2.67 Mbp. No clear p-value inflation or deflation
with several other QTL on different chromosomes (Figure 5A). was observed compared to the expected p-value distribution
Finally, two QTL for SdFe (SdFe_2_46.07 and SdFe_6_22.37) (Supplementary Figure 10B). In summary, further SNPs
showed an effect in all three panels which had data for SdFe were identified in the VEC separately to specifically improve
(Figure 4). Except for DF, major QTL for 100SdW, SdFe, and climbing beans.
yield were identified showing effects in all five breeding panels.
3.6. Genomic Prediction
The GEBV for new VEC climbing bean lines (validation set) were
3.4.3. Haplotype Effects Across and Among Growth calculated either with parts of the VEC or with parts of all five
Types panels as TP.
For traits with complex genetic architecture, single SNPs poorly
explain the variance caused by the associated genetic region. 3.6.1. Genomic Prediction Among Environments
Therefore, haplotypes were constructed on Chr 1 between The PAs for new VEC lines within each trial and across all trials
43.71 and 45.47 Mbp for SNPs significantly associated either differed between the traits as well as the used model approaches
with DF or growth type. The haplotypes for DF on Chr 1 (Figure 6). In general, PAs followed the heritability and the
explained 44.6% of the variance for DF across all growth types genomic r calculated using all available lines as TP. PAs for yield
(Supplementary Figure 8A). In contrast, the best SNP for DF reached about r ≈ 0.5 in the high-altitude locations Darién
on Chr 1 explained only 8% of the variance. The haplotypes and Popayán, where climbing beans are better adapted. The PAs
for growth type on Chr 1 differentiated type I from the other for yield in other locations were lower. The trials Dar14B and
types (Supplementary Figure 8B). The first haplotype (“101”) Pop15B showed lower PAs with higher variability because these
was almost exclusively associated with the determinate growth trials comprised only 100 lines.
type I. The second haplotype (“001”) was mainly associated On average, the FA model showed the best performance for
with growth types II and IV. The remaining three haplotypes DF. The genotype model for single trials performed best for
were associated with climbing growth habit (types III and IV). 100SdW, SdFe, and yield. The FA andGxEmodel reached slightly
An important breeding goal is to shorten DF in all growth lower average PAs for SdFe and yield, respectively. The strong
types while maintaining other agronomic traits. Therefore, all GxE for yield was reflected in the different marker effects among
SNPs in the region from 44.60 to 45.47 Mbp on Chr 1 were the locations (Supplementary Figure 11). The PAcc reached the
clustered, resulting in ten distinct haplotypes according to the highest values for 100SdW using the genotype model among
average silhouette width (Figure 5B). As expected, the haplotypes trials (77.4%), DF using the FA model (67.4%), SdFe using the
showed strong effects on DF and explained almost 50% of the genotype model for single trials (63.5%), and yield using the
variation for DF. However, these haplotype effects were not genotype model for single trials(72.5%, Supplementary Table 7).
consistent among the growth types, explaining 3.7, 36.9, 0.0, In summary, the PAs differed among models and traits, i.e.,
and 5.3% of the DF variation in growth types I, II, III, and when predicting for a single trial, the genotype model performed
IV, respectively. In addition, the effect on the remaining traits best, except for DF. When predicting for multiple years in one
varied substantially across the haplotypes. Since only a few location, the FA model was promising for DF and SdFe while for
SNPs expressed significant pleiotropic effects, some trade-offs yield and 100SdW, the GxE models performed best.
between traits could be removed by breaking the LD of QTL on
different chromosomes (Figure 5A). In summary, the QTL on 3.6.2. Genomic Prediction Across Environments With
Chr 1 between 44.60 to 45.47 Mbp controlled major processes Optimization of the Training Population
across the growth types but showed minor effects among them. To increase PAs for new VEC lines, three different approaches
Furthermore, the QTL exhibited varying pleiotropic effects on were tested: when only VEC lines were in the TP (TP VEC),
SdFe, 100SdW, and yield which can be decreased by breaking the when all lines of the five panels were in the TP (TP extended),
LD between this and further QTL (e.g., Yd_7_4.86). or when distantly related lines were excluded (TP optimized;
Frontiers in Plant Science | www.frontiersin.org 10 April 2022 | Volume 13 | Article 830896
Keller et al. Improving Genomic Predictions for Beans
TABLE 1 | Significant marker-trait associations across five bean breeding panels below the 1% significance level according to genome-wide association studies.
QTL ID Trait Chr Pos p-value Marker MAF
GH_1_0.53 Growth habit 1 526025 1.18E-10 Pv2.1_01_526025_A/T 0.24
GH_1_2.85 Growth habit 1 2850955 2.78E-10 Pv2.1_01_2850955 C/G 0.38
GH_1_43.71 Growth habit 1 43707108 1.93E-35 Pv2.1_01_43707108 A/G 0.41
GH_1_45.04 Growth habit 1 45044047 3.90E-36 Pv2.1_01_45044047 C/T 0.13
GH_1_45.37 Growth habit 1 45374662 4.01E-24 Pv2.1_01_45374662 T/C 0.08
GH_1_47.44 Growth habit 1 47439968 1.23E-08 Pv2.1_01_47439968 G/C 0.16
GH_2_32.21 Growth habit 2 32208182 9.73E-08 Pv2.1_02_32208182 C/G 0.32
GH_2_40.05 Growth habit 2 40045968 6.88E-18 Pv2.1_02_40045968 C/A 0.05
GH_3_1.29 Growth habit 3 1289075 2.35E-08 Pv2.1_03_1289075 T/A 0.45
GH_3_1.92 Growth habit 3 1921492 1.22E-09 Pv2.1_03_1921492 A/G 0.03
GH_3_42.49 Growth habit 3 42488118 9.24E-09 Pv2.1_03_42488118 G/C 0.3
GH_3_44.06 Growth habit 3 44056005 2.62E-09 Pv2.1_03_44056005 A/T 0.12
GH_4_0.45 Growth habit 4 448766 2.45E-07 Pv2.1_04_448766 A/T 0.03
GH_4_1.42 Growth habit 4 1421295 2.52E-09 Pv2.1_04_1421295 A/G 0.47
GH_4_1.92 Growth habit 4 1919859 1.29E-09 Pv2.1_04_1919859 A/G 0.26
GH_4_2.16 Growth habit 4 2164168 9.82E-08 Pv2.1_04_2164168 T/C 0.08
GH_4_2.56 Growth habit 4 2559941 1.04E-11 Pv2.1_04_2559941 A/T 0.25
GH_4_47.17 Growth habit 4 47174835 1.95E-07 Pv2.1_04_47174835 T/A 0.26
GH_5_0.39 Growth habit 5 394641 7.38E-10 Pv2.1_05_394641 T/G 0.13
GH_5_0.74 Growth habit 5 739798 1.33E-16 Pv2.1_05_739798 G/A 0.28
GH_5_10.7 Growth habit 5 10696009 1.63E-07 Pv2.1_05_10696009 G/C 0.02
GH_6_22.3 Growth habit 6 22301003 3.03E-15 Pv2.1_06_22301003 A/G 0.03
GH_6_22.51 Growth habit 6 22508433 3.94E-16 Pv2.1_06_22508433 T/G 0.16
GH_6_23.87 Growth habit 6 23868436 1.66E-08 Pv2.1_06_23868436 G/C 0.37
GH_6_26.05 Growth habit 6 26054074 2.23E-07 Pv2.1_06_26054074 C/T 0.1
GH_6_29.43 Growth habit 6 29429040 1.53E-08 Pv2.1_06_29429040 G/A 0.02
GH_7_3.05 Growth habit 7 3047903 1.18E-09 Pv2.1_07_3047903 G/A 0.28
GH_7_7.15 Growth habit 7 7150019 6.01E-10 Pv2.1_07_7150019 C/T 0.1
GH_7_38.99 Growth habit 7 38987037 3.73E-10 Pv2.1_07_38987037 A/G 0.42
GH_8_2.11 Growth habit 8 2107245 2.42E-07 Pv2.1_08_2107245 C/T 0.23
GH_9_0.86 Growth habit 9 860918 6.16E-07 Pv2.1_09_860918 G/A 0.17
GH_9_13.92 Growth habit 9 13924731 3.14E-08 Pv2.1_09_13924731 T/C 0.49
GH_9_34.74 Growth habit 9 34742076 1.48E-07 Pv2.1_09_34742076 A/G 0.23
GH_9_36.11 Growth habit 9 36110952 1.03E-07 Pv2.1_09_36110952 A/T 0.01
GH_10_3.19 Growth habit 10 3191949 6.26E-08 Pv2.1_10_3191949 G/T 0.21
GH_10_6.5 Growth habit 10 6495216 2.67E-09 Pv2.1_10_6495216 G/T 0.27
GH_11_1.01 Growth habit 11 1010422 1.14E-08 Pv2.1_11_1010422 T/C 0.49
GH_11_2.78 Growth habit 11 2775768 9.58E-11 Pv2.1_11_2775768 T/G 0.16
DF_1_41.08 DF 1 41082526 7.72E-09 Pv2.1_01_41082526 G/A 0.22
DF_1_44.6 DF 1 44604072 1.88E-08 Pv2.1_01_44604072 A/C 0.35
DF_1_44.93 DF 1 44927394 5.64E-09 Pv2.1_01_44927394 C/T 0.25
DF_1_45.04 DF 1 45044047 5.11E-08 Pv2.1_01_45044047 C/T 0.13
DF_1_45.23 DF 1 45233651 3.51E-09 Pv2.1_01_45233651 A/G 0.08
DF_1_45.47 DF 1 45469012 1.21E-33 Pv2.1_01_45469012 G/A 0.1
DF_2_7.9 DF 2 7900892 1.22E-08 Pv2.1_02_7900892 A/C 0.26
DF_2_31.53 DF 2 31525478 4.99E-09 Pv2.1_02_31525478 A/C 0.16
DF_4_45.98 DF 4 45975326 5.51E-08 Pv2.1_04_45975326 T/C 0.38
DF_7_4.82 DF 7 4824001 5.29E-09 Pv2.1_07_4824001 C/G 0.26
DF_8_15.48 DF 8 15481886 1.12E-08 Pv2.1_08_15481886 T/A 0.01
(Continued)
Frontiers in Plant Science | www.frontiersin.org 11 April 2022 | Volume 13 | Article 830896
Keller et al. Improving Genomic Predictions for Beans
TABLE 1 | Continued
QTL ID Trait Chr Pos p-value Marker MAF
DF_9_23.12 DF 9 23116201 9.41E-09 Pv2.1_09_23116201 C/T 0.28
DF_9_29.38 DF 9 29378513 1.25E-10 Pv2.1_09_29378513 C/T 0.07
DF_9_37.83 DF 9 37825106 3.55E-09 Pv2.1_09_37825106 T/G 0.21
SdW_1_17.16 100SdW 1 17164005 1.19E-07 Pv2.1_01_17164005 A/G 0.16
SdW_1_42.52 100SdW 1 42515535 5.49E-10 Pv2.1_01_42515535 A/G 0.07
SdW_1_47.26 100SdW 1 47260087 7.66E-09 Pv2.1_01_47260087 A/G 0.28
SdW_2_2.23 100SdW 2 2234763 3.64E-08 Pv2.1_02_2234763 T/G 0.45
SdW_4_46.54 100SdW 4 46540184 1.51E-08 Pv2.1_04_46540184 C/G 0.36
SdW_5_1.06 100SdW 5 1062441 1.09E-07 Pv2.1_05_1062441 C/A 0.44
SdW_5_4.07 100SdW 5 4065359 8.94E-11 Pv2.1_05_4065359 C/T 0.28
SdW_6_18.46 100SdW 6 18456447 1.48E-08 Pv2.1_06_18456447 G/A 0.2
SdW_6_25.25 100SdW 6 25248245 1.34E-08 Pv2.1_06_25248245 G/C 0.21
SdW_6_28.9 100SdW 6 28898246 3.91E-13 Pv2.1_06_28898246 G/A 0.16
SdW_7_28.77 100SdW 7 28765036 1.85E-11 Pv2.1_07_28765036 C/T 0.02
SdW_9_28.29 100SdW 9 28287352 6.16E-07 Pv2.1_09_28287352 G/A 0.12
SdW_11_1.55 100SdW 11 1547189 1.39E-09 Pv2.1_11_1547189 C/T 0.12
SdFe_2_46.07 SdFe 2 46068228 3.20E-08 Pv2.1_02_46068228 G/C 0.19
SdFe_6_22.37 SdFe 6 22365971 3.95E-11 Pv2.1_06_22365971 A/C 0.32
SdFe_9_36 8 SdFe 9 36804490 6.26E-08 Pv2.1_09_36804490 C/T 0.16
Yd_7_4.86 Yield scaled 7 4856975 3.28E-08 Pv2.1_07_4856975 C/T 0.23
The QTL ID, trait, chromosome (Chr), physical position (Pos) in bp, association strength (p value), name of the marker, and minor allele frequency (MAF) of the significantly associated
SNPs is reported. The panels included bush and climbing growth habits and were evaluated for days to flowering (DF), 100 seed weight (100SdW), seed iron concentration (SdFe), and
seed yield (Yield_scaled, scaled for each panel).
Supplementary Figure 12). The optimized TP increased PAs for joint group of panels with fast LD decay enabled the identification
DF, SdFe, and yield (scaled among panels) when adding bush type of SNPs tightly linked to the causal loci while controlling for
lines of other panels and reached a PAcc of 66.8, 66.6, and 22.7% population structure (Sul et al., 2018; Huang et al., 2019).
corresponding to a 0.7, 1.8, and 8.8% increase on the averaged PA, On the one hand, several QTL for growth habit and DF
respectively (Figure 6). Regarding 100SdW, the TP optimization were confirmed from previous studies. On the other hand,
and extension decreased PAs slightly compared to the TP with especially for SdFe and 100SdW, new QTL and candidate genes
only VEC lines. In summary, in complex traits such as SdFe and were identified.
yield, the PA can be improved by adding related lines from other
panels which are not in the TP even though they were tested in 4.1.1. Previously Described and New QTL for DF and
different trials. Growth Habit
Considering significant marker-trait associations less than 1
4. DISCUSSION Mbp away from previously reported positions, QTL for DF
and growth habit were confirmed on Chr 1, 4, 9, and 11. In
Based on the largest assembly of phenotypic and genotypic addition, on Chr 6, 7, and 8, significantly associated SNPs
common bean data, we showed increased PAs for important traits were mapped to a distance of 1 to 4 Mbp from previously
of climbing bean by the addition of related bush type beans reported QTL. The terminal flowering gene PvTFL1y (fin
from other trials to the TP. In addition, the extended pool of locus identified as Phvul.001G189200 on Chr 1 at 44.85
lines, including 1,869 genotypes from distinct breeding panels, Mbp) (Norton, 1915; Koinange et al., 1996; Kwak et al.,
was useful to predict growth type and to increase power in the 2008; Repinski et al., 2012; González et al., 2016) was
detection of QTL using GWAS (Spindel and McCouch, 2016). confirmed by DF_1_44.60 (Table 1). The phytochrome A
Hence, this comprehensive study provides a solid basis to harness gene (Ppd locus identified as Phvul.001G221100 on Chr 1
the large genetic diversity of common bean germplasm and to at 47.64 Mbp) conferring photoperiod sensitivity (Coyne
implement marker-assisted and genomic selection strategies for and Schuster, 1974; Gu et al., 1998; Kamfwa et al., 2015;
more efficient climbing bean breeding. Weller et al., 2019) was tightly linked to GH_1_47.43. In
agreement, the haplotype “101” constructed in the region
4.1. QTL Across Breeding Panels from 43.71 to 45.37 Mbp almost exclusively differentiated
For all studied traits, QTL with clear effects on the phenotypes between determinate and indeterminate growth types
in all five breeding panels were detected (Figure 4). This diverse (Supplementary Figure 8B). Thus, multiple SNPs are required
Frontiers in Plant Science | www.frontiersin.org 12 April 2022 | Volume 13 | Article 830896
Keller et al. Improving Genomic Predictions for Beans
Allele dosage 0 2 Panel ADP AxM MIP VEC VEF
60 60
50 50
40 40
30 30
60 60
50 50
40 40
30 30
20 20
10 10
100 100
80 80
60 60
40 40
2.5
0.0
−2.5
−5.0
FIGURE 4 | Boxplots for allele dosage effect (0 or 2 alternative alleles) of the most significant SNPs associated with days to flowering (DF), 100 seed weight
(100SdW), seed iron concentration (SdFe), and seed yield for five breeding panels. The panels included lines with bush and climbing growth habits originating from the
Andean and Mesoamerican gene pools. The trait, the name of the associated marker and the QTL ID is given.
Frontiers in Plant Science | www.frontiersin.org 13 April 2022 | Volume 13 | Article 830896
Yield [scaled] SdFe [mg/kg] 100SdW [g] DF [d]
Yd_7_4.86 SdFe_2_46.07 SdW_6_28.9 DF_1_45.47
Pv2.1_07_4856975_CT Pv2.1_02_46068228_GC Pv2.1_06_28898246_GA Pv2.1_01_45469012_GA
SdFe [mg/kg] 100SdW [g] DF [d]
SdFe_6_22.37 SdW_7_28.77 DF_9_29.38
Pv2.1_06_22365971_AC Pv2.1_07_28765036_CT Pv2.1_09_29378513_CT
Keller et al. Improving Genomic Predictions for Beans
A
GH_1_45.37 SdW_6_28.9
−log10(p)
DF_9_29.38 DF_1_45.47
10
15
DF_1_45.23
20
25 GH_1_45.04
30
SdW_7_28.77
35
DF_1_44.6 GH_2_40.05
Trait
Growth habit GH_1_43.71 SdFe_2_46.07
DF
100SdW
SdFe
Yield_scaled DF_1_44.93
SdFe_6_22.37
Yd_7_4.86
B Panel ADP AxM MIP VEC VEF
Growth type I Growth type II Growth type III Growth type IV
50
40
60
50
40
30
20
100
80
60
40
2000
1000
Haplotype
FIGURE 5 | Significant marker-trait associations were analyzed for growth habit, days to flowering (DF), 100 seed weight (100SdW), seed iron concentration (SdFe),
and seed yield scaled among the five breeding panels (Yield_scaled). (A) A network of the SNPs significantly (below the 1% significance level) associated with each of
the four traits forming clusters according to their linkage disequilibrium is shown. Each dot represents a significant SNP and its size the associated −log10 p-value.
The SNPs on Chr 1 between 43.71 and 45.47 as well as the two most significant SNPs per trait are labeled with the QTL ID. (B) Haplotypes including all SNPs
between 44.60 and 45.47 Mbp on chromosome 1 were constructed using hierarchical clustering. The averaged SNP effects of the haplotypes were evaluated in all
traits among all growth types, i.e., growth type I (determinate bush type), type II (indeterminate bush), type III (determinate climber), and type IV (indeterminate climber).
The error bars show the SD.
Frontiers in Plant Science | www.frontiersin.org 14 April 2022 | Volume 13 | Article 830896
Phenotype
Yield [kg/ha] SdFe [mg/kg] 100SdW [g] DF [d]
1
2
3
4
5
6
7
8
9
10
1
2
3
4
5
6
7
8
9
10
1
2
3
4
5
6
7
8
9
10
1
2
3
4
5
6
7
8
9
10
Keller et al. Improving Genomic Predictions for Beans
Genotype (single trials) GxE TP extended TP VEC Genomic r
Genotype (among trials) Factor analysis TP optimized
1.00
0.75
0.50
0.25
1.00
0.75
0.50
0.25
0.00
1.00
0.75
0.50
0.25
0.00
0.5
0.0
−0.5
Dar14B Dar18B Dar19B Pop15B Pop17D Pal19D TzKg19D UgKw19D BLUEs
FIGURE 6 | Prediction abilities for different traits and models for new lines from the climbing bean panel (VEC). The first-stage best linear unbiased estimators (BLUEs)
for each trial or second-stage BLUEs across trials were used. The first-stage BLUEs models were tested accounting for either genotypic effects only, i.e., for each trial
separately (single trials) or all together (among trials), genotypic x environment interaction (GxE), or correlation between trials (Factor analysis). The second-stage
BLUEs were used for models which take into account genotypic effects based on the VEC (TP VEC), on all five panels (TP extended), and all five panels with
optimization of the training population (TP optimized). The predicted traits were days to flowering (DF), 100 seed weight (100SdW), seed iron concentration (SdFe),
and seed yield. Seed yield was scaled among panels for the models based on second-stage BLUEs to allow comparison between the different growth habits. The
horizontal line is the square root of heritability indicating the heritable variance of the trait in each trial. Trials were abbreviated based on the location Darién (Dar),
Palmira (Pal), Popayán (Pop) in Colombia, Kagera in Tanzania (TzKg), or Kawanda in Uganda (UgKw), the year, and the planting season (sequentially A to D). For a
detailed description of each trial refer to Supplementary Table 2.
to determine growth type including the allelic version of the In our study, GH_7_38.99 was detected proximal to PvTFL1z. On
PvTFL1y gene. the upper arm of Chr 8 at 4.9 Mbp, another QTL for DF was
In the region of the two identified SNPs for growth habit reported previously (Raggi et al., 2019). Similarly, we detected
on Chr 4 (GH_4_1.42 and GH_4_1.92), a QTL associated with GH_8_2.11 at less than 3 Mbp distance. A second fin locus
climbing ability and plant height was previously reported (linked (fin’) on Chr 9 was tagged by molecular markers at 13.39Mbp
to Pvctt001 marker at 0.51 Mbp) (Checa and Blair, 2008). (de Campos et al., 2011) and at around 20 cM (González et al.,
Furthermore, in proximity to GH_6_29.43, a QTL for DF was 2016). This fin’ locus is probably tagged by GH_9_13.92. A QTL
previously reported at 31.6Mbp (Raggi et al., 2019). The terminal for DF on Chr 11 at around 9 cM reported in Bhakta et al. (2017)
flowering gene PvTFL1z (Phvul.007G229300 a homolog of was linked to GH_11_1.01 (at around 8.8 cM). A new QTL for
PvTFL1y) was located on Chr 7 at 35.31 Mbp (Kwak et al., 2008). growth habit, GH_2_40.05 was identified with one SNP variant
Frontiers in Plant Science | www.frontiersin.org 15 April 2022 | Volume 13 | Article 830896
Prediction ability (r)
Yield SdFe 100SdW DF
Keller et al. Improving Genomic Predictions for Beans
exclusively associated to climbing beans, however, with a low major effects in all breeding panels (Figure 4). The asparagine
MAF of 5%. Newly identified SNPs for growth habit with lower synthetase remobilizes nitrogen from sources to sinks and was
LOD scores have to be interpreted carefully due to the observed reported to increase seed weight and soluble protein content in
p-value inflation, indicating remaining population structure. In Arabidopsis seed (Gaufichon et al., 2016, 2017). In conclusion,
conclusion, the joint analysis consisting of diverse common bean for the major QTL for SdFe and 100SdW, plausible candidate
populations showed high detection power for DF and growth genes were identified whose putative functions remain to be
habit QTL. Four to potentially seven QTL known from previous further validated.
studies and several new QTL for DF and growth habit were
identified and tagged by tightly linked SNP markers. 4.2. QTL Detected Within Breeding Panels
Several QTL were identified only in the VEC, e.g., a pleiotropic
4.1.2. QTL for SdFe and Their Independence From QTL for DF and PHI on Chr 5 between 38.68 and 39.34 Mbp.
Growth Habit The QTL for PHI differed from QTL previously identified in bi-
Four meta-QTL were previously reported for SdFe on Chr 1 parental bush type populations (Mukeshimana et al., 2014; Diaz
(between 43.3 and 48.5 Mbp), on Chr 6 (between 28.2 and 29.5 et al., 2018). The PHI was weakly and positively correlated to
Mbp), on Chr 9 (between 11.7 and 13.5 Mbp), and Chr 11 DF and yield (Supplementary Figure 3). This pleiotropic effect
(between 2.3 and 5.3 Mbp) (Izquierdo et al., 2018). All four of the identified QTL might be related to the time from flowering
QTL fall right onto or next to QTL for growth habit identified until harvest which could affect DF and PHI. Regarding DF, the
in the current study. Since climbing beans in general exhibit a major QTL tagged by DF_1_44.60 was mapped more closely to
higher SdFe (Blair et al., 2010; Petry et al., 2015), the previously PvTFL1y in the current joint analysis than in a separate analysis of
reported QTL is probably confounded with population structure. the ADP and VEF at 48.34 and 49.72 Mbp, respectively (Kamfwa
The effects of these QTL on SdFe were also detected in our et al., 2015; Keller et al., 2020). In agreement, the LD decay of
analysis but the associations did not exceed the 5% significance the combined panel was faster than that of the separate panels.
threshold (Supplementary Figure 7B). A QTL for SdZn was Furthermore, we concluded that the genetic control of flowering
recently reported on Chr 1 at 49.37 Mbp in European landraces is different on the single SNP level (Figure 4) and haplotype level
next to PvTFL1y, suggesting a genetic linkage between DF and among the growth types (Figure 5B). The major QTL for DF
SdZn (Caproni et al., 2020). From the three major QTL for SdFe on Chr 1 showed no effect within growth type III lines. One
detected in our study, two were identified for the first time and possible explanation is that DF within the climbing growth habit
SdFe_6_22.37 was previously reported in proximity at 22.8 Mbp is regulated by additional genes or growth habit specific alleles.
(Diaz et al., 2020). The identified SNPs for SdFe on Chr 2 and 6 Additionally, haplotype 1 showed increased yield for growth
showed an effect in all three evaluated breeding panels and are types I and II while the effect on DF differed. It demonstrates that
thus independent of the growth habit (Figure 4). The QTL at the DF and yield are loosely linked as suggested earlier byWhite et al.
end of Chr 2 (SdFe_2_46.07) is of particular interest because it (1992). Similarly, the QTL for SdFe on Chr 2 at 2.91 and 48.86
was linked to a QTL for DPM in this study and a pleiotropic QTL Mbp were valid only for the VEC and were not detected in the
affecting DF, DPM, 100SdW, and yield in the VEF (Keller et al., joint analysis suggesting some specific alleles were present only
2020). In conclusion, various QTL for SdFe reported previously in the climbing bean panel.
seemed to be confounded with growth habit, while our joint
analysis allowed us to detect new SNPs associated with SdFe 4.3. Genomic Predictions Across and
across different breeding panels. Within Breeding Panels
Using only VEC lines as TP in the different modeling approaches,
4.1.3. Candidate Gene Identification for SdFe and the FA model performed well for DF and SdFe, showing
100SdW the importance of covariance between trials for those traits
Four new candidate genes were identified for SdFe and 100SdW: (Figure 6). Such FAmodels work well for different environments,
on Chr 6 at 22.18 Mbp, in a distance of less than 200,000 bp from where some lines were already tested, were previously reported
the SNP most significantly associated to SdFe (SdFe_6_22.37), for sweet cherry (Prunus avium L.) (Hardner et al., 2019).
the Phvul.006G113100 gene was annotated as a homologous The PAcc for yield of the genotype model using the first-
to a ferric-chelate reductase, reported to be involved in iron and second-stage BLUEs among all environments were lower
uptake from the soil (Robinson et al., 1999; Wu et al., 2005; in comparison to the FA model and GxE model which took
Jeong et al., 2008; Asard et al., 2013). In proximity, less than environmental differences into account. In agreement, the high
125,000 bp away from SdFe_2_46.07 and SdFe_9_36.80, the GxE was visible in the correlations between the trials which
genes Phvul.002G292900 and Phvul.009G247600, respectively, reached values from –0.41 to 0.64 (Figure 1B). Those interactions
were annotated. These two genes putatively expressAtox1-related were additionally reflected in the relatively high PA of the
copper transport proteins, which are involved in copper and iron GxE model in yield and in the differing marker effects among
homeostasis (Himelblau et al., 1998; Puig et al., 2007). A QTL locations (Figure 6 and Supplementary Figure 11). Comparing
for copper and iron uptake was shown previously to have close different traits, the correlations between yield and DPM in
genetic linkage (Waters and Grusak, 2008). Regarding 100SdW, bush type beans altered when changing from low- to high-
a putative asparagine synthetase (Phvul.006G188400) on Chr 6 at altitude locations (Diaz et al., 2018). Therefore, the single-
28.87 Mbp, less than 30,000 bp away from SdW_6_28.90, showed trial genotype model and the GxE model, which account for
Frontiers in Plant Science | www.frontiersin.org 16 April 2022 | Volume 13 | Article 830896
Keller et al. Improving Genomic Predictions for Beans
location effects such as high- and lowland, worked best for the AUTHOR CONTRIBUTIONS
predictions of yield. Similarly, high GxE for yield was shown
before in the VEF under different environmental conditions BR and BS conceived the study. CM, AP-B, HB, WA, SM,
(Keller et al., 2020). Across all trials, the TP optimization turned JM, TP, JB, and PM conducted the experiments and provided
out to be a promising strategy to predict the performance of the phenotypic data. BK, DA-S, and JA performed modeling,
new climbing bean lines. For complex traits such as SdFe and data analysis, and interpretation. DA-S prepared the genotypic
yield, optimizing the TP with related lines from different panels data. BK drafted the manuscript, which was improved by
tested under different conditions improved the PA by 1.8 and DA-S, BS, PM, TP, and BR. All authors approved to the
8.8%, respectively. Especially for smaller breeding programs manuscript submission.
or new breeding panels, it is, therefore, advisable to optimize
the TP even when the added lines were tested in different FUNDING
field trials.
The authors would like to thank the COOP Research Program
on Sustainability in Food Value Chains of the ETH Zurich
5. CONCLUSION World Food System Center and the ETH Zurich Foundation
for supporting this project. The COOP Research Program
The benefits of introducing genes conferring climbing ability into was supported by the COOP Sustainability Fund. The work
new breeding lines are limited because growth habits are fixed in was additionally founded by Tropical Legumes III-Improving
specific production systems and their modification would have Livelihoods for Smallholder Farmers: Enhanced Grain Legume
pleiotropic effects, most likely affecting traits like DF and SdFe Productivity and Production in Sub-Saharan Africa and South
negatively. The most significant QTL for growth habit at the end Asia (OPP1114827), and by the AVISA-Accelerated varietal
of Chr 1 was pleiotropic. However, this QTL was in LD with improvement and seed delivery of legumes and cereals in
several other QTL which can be selected separately since they Africa (OPP1198373) projects funded by the Bill & Melinda
were located on different chromosomes (Figure 5A). In addition, Gates Foundation. The Norman Borlaug Cooperative Research
this QTL had no effect on common beans of the growth type Initiative (NBCRI) Grain Legumes Project of the US Agency
III. Regarding other QTL, we detected stable SdFe and yield for International Development, Project No. 0210-22310-004-96R
QTL which showed effects in all tested panels without significant supported the Tanzanian trials in Arusha and Mbeya. We would
pleiotropic effects. The identified markers were validated in very like to thank USAID for their contributions through the CGIAR
diverse germplasm including all growth types and both gene Research Program on Grain Legumes and Dryland Cereals. Open
pools. The resulting fast LD decay allowed mapping of QTL access publication is kindly covered by ETH Zurich.
more precisely than was achieved in separated panels. Genomic
prediction models were established across populations enabling ACKNOWLEDGMENTS
the selection and hybridization of the best lines of all populations
to combine favorable alleles. Especially for yield, GxE needs to We very much thank the CIAT bean teams for supporting and
be modeled when breeding for different environments. The large carrying out the experiments. Special thanks to Brenda Nakyanzi,
common bean diversity presented in this study was used to Eninka Mndolwa, and others for managing the field trials. Many
identify markers across populations and to establish and improve thanks to Lucy Milena Diaz Martinez and Eliana Macea for DNA
prediction models. The joint population will provide a basis to extraction and library preparation.We further acknowledge Prof.
exploit this genetic diversity and will contribute to the quick and Achim Walter from ETH Zurich for sharing infrastructure with
targeted development of new (climbing bean) lines. the Molecular Plant Breeding group.
DATA AVAILABILITY STATEMENT SUPPLEMENTARY MATERIAL
The original contributions presented in the study are publicly The Supplementary Material for this article can be found
available. This data can be found here: https://doi.org/10.7910/ online at: https://www.frontiersin.org/articles/10.3389/fpls.2022.
DVN/RLAWYN. 830896/full#supplementary-material
REFERENCES yield of climbing bean. Plant Soil. 428, 223–239. doi: 10.1007/s11104-018-
3665-y
Akdemir, D., and Isidro-Sánchez, J. (2019). Design of training populations Bhakta, M. S., Gezan, S. A., Michelangeli, J. A. C., Carvalho, M., Zhang, L.,
for selective phenotyping in genomic prediction. Sci. Rep. 9, 1446. Jones, J. W., et al. (2017). A predictive model for time-to-flowering in the
doi: 10.1038/s41598-018-38081-6 common bean based on QTL and environmental variables. G3 7, 3901–3912.
Asard, H., Barbaro, R., Trost, P., and Bérczi, A. (2013). Cytochromes b561: doi: 10.1534/g3.117.300229
ascorbate-mediated trans-membrane electron transport. Antioxidants Redox Bitocchi, E., Bellucci, E., Giardini, A., Rau, D., Rodriguez, M., Biagetti, E., et al.
Signal. 19, 1026–1035. doi: 10.1089/ars.2012.5065 (2013). Molecular analysis of the parallel domestication of the common bean
Barbosa, N., Portilla, E., Buendia, H. F., Raatz, B., Beebe, S., and Rao, I. (Phaseolus vulgaris) in Mesoamerica and the Andes.New Phytol. 197, 300–313.
(2018). Genotypic differences in symbiotic nitrogen fixation ability and seed doi: 10.1111/j.1469-8137.2012.04377.x
Frontiers in Plant Science | www.frontiersin.org 17 April 2022 | Volume 13 | Article 830896
Keller et al. Improving Genomic Predictions for Beans
Black, R. E., Allen, L. H., Bhutta, Z. A., Caulfield, L. E., Onis, M., d., et al. (2008). a MAGIC population of common bean (Phaseolus vulgaris L.) under drought
Maternal and child undernutrition: global and regional exposures and health conditions. BMC Genomics 21, 799. doi: 10.1186/s12864-020-07213-6
consequences. Lancet 371, 243–260. doi: 10.1016/S0140-6736(07)61690-0 Elshire, R. J., Glaubitz, J. C., Sun, Q., Poland, J. A., Kawamoto, K., Buckler, E. S.,
Blair, M.W. (2013). Mineral biofortification strategies for food staples: the example et al. (2011). A robust, simple genotyping-by-sequencing (GBS) approach for
of common bean. J. Agric. Food Chem. 61, 8287–8294. doi: 10.1021/jf400774y high diversity species. PLoS ONE 6, e19379. doi: 10.1371/journal.pone.0019379
Blair, M.W., Prieto, S., Diaz, L.M., Buendia, H. F., and Cardona, C. (2010). Linkage Endelman, J. B. (2011). Ridge Regression and other kernels for
disequilibrium at the APA insecticidal seed protein locus of common bean genomic selection with r package rrBLUP. Plant Genome 4, 250–255.
(Phaseolus vulgaris L.). BMC Plant Biol. 10, 15. doi: 10.1186/1471-2229-10-79 doi: 10.3835/plantgenome2011.08.0024
Bliss, F. A. (1993). “Breeding common bean for improved biological nitrogen Endelman, J. B., and Jannink, J.-L. (2012). Shrinkage estimation of the realized
fixation,” in Enhancement of Biological Nitrogen Fixation of Common Bean in relationship matrix. G3 2, 1405–1413. doi: 10.1534/g3.112.004259
Latin America: Results from an FAO/IAEA Co-ordinated Research Programme, Fang, C., Ma, Y., Wu, S., Liu, Z., Wang, Z., Yang, R., et al. (2017). Genome-wide
1986–1991, Developments in Plant and Soil Sciences, eds F. A. Bliss and G. association studies dissect the genetic networks underlying agronomical traits
Hardarson (Dordrecht: Springer Netherlands), 71–79. in soybean. Genome Biol. 18, 161. doi: 10.1186/s13059-017-1289-9
Boccio, J. R., and Iyengar, V. (2003). Iron deficiency. Biol. Trace Elem. Res. 94, 1–31. Gaufichon, L., Marmagne, A., Belcram, K., Yoneyama, T., Sakakibara, Y.,
doi: 10.1385/BTER:94:1:1 Hase, T., et al. (2017). ASN1-encoded asparagine synthetase in floral organs
Bolger, A. M., Lohse, M., and Usadel, B. (2014). Trimmomatic: a flexible contributes to nitrogen filling in Arabidopsis seeds. Plant J. 91, 371–393.
trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120. doi: 10.1111/tpj.13567
doi: 10.1093/bioinformatics/btu170 Gaufichon, L., Rothstein, S. J., and Suzuki, A. (2016). Asparagine
Briatte, F., Bojanowski, M., Canouil, M., Charlop-Powers, Z., Fisher, J. C., Johnson, metabolic pathways in Arabidopsis. Plant Cell Physiol. 57, 675–689.
K., et al. (2020). ggnetwork: Geometries to Plot Networks with ‘ggplot2’. doi: 10.1093/pcp/pcv184
Available online at: https://CRAN.R-project.org/package=ggnetwork Gepts, P., Osborn, T. C., Rashka, K., and Bliss, F. A. (1986). Phaseolin-protein
Browning, B. L., Zhou, Y., and Browning, S. R. (2018). A one-penny imputed variability in wild forms and landraces of the common bean (Phaseolus
genome from next-generation reference panels. Am. J. Hum. Genet. 103, vulgaris): evidence for multiple centers of domestication. Econ. Bot. 40,
338–348. doi: 10.1016/j.ajhg.2018.07.015 451–468. doi: 10.1007/BF02859659
Caproni, L., Raggi, L., Talsma, E. F., Wenzl, P., and Negri, V. (2020). González, A. M., Yuste-Lisbona, F. J., Saburido, S., Bretones, S., De Ron, A.
European landrace diversity for common bean biofortification: a genome-wide M., Lozano, R., et al. (2016). Major contribution of flowering time and
association study. Sci. Rep. 10, 19775. doi: 10.1038/s41598-020-76417-3 vegetative growth to plant production in common bean as deduced from a
Checa, O., Ceballos, H., and Blair, M. W. (2006). Generation means analysis comparative genetic mapping. Front. Plant Sci. 7, 1940. doi: 10.3389/fpls.2016.
of climbing ability in common bean (Phaseolus vulgaris L.). J. Heredity 97, 01940
456–465. doi: 10.1093/jhered/esl025 Graham, P. H. (1981). Some problems of nodulation and symbiotic nitrogen
Checa, O. E., and Blair, M. W. (2008). Mapping QTL for climbing ability and fixation in Phaseolus vulgaris L.: a review. Field Crops Res. 4, 93–112.
component traits in common bean (Phaseolus vulgaris L.). Mol. Breed. 22, doi: 10.1016/0378-4290(81)90060-5
201–215. doi: 10.1007/s11032-008-9167-5 Gu, W., Zhu, J., Wallace, D. H., Singh, S. P., and Weeden, N. F. (1998). Analysis of
Cichy, K. A., Fernandez, A., Kilian, A., Kelly, J. D., Galeano, C. H., genes controlling photoperiod sensitivity in common bean usingDNAmarkers.
Shaw, S., et al. (2014). QTL analysis of canning quality and color Euphytica 102, 125–132. doi: 10.1023/A:1018340514388
retention in black beans (Phaseolus vulgaris L.). Mol. Breed. 33, 139–154. Guild, G. E., Paltridge, N. G., Andersson, M. S., and Stangoulis, J. C. R. (2017).
doi: 10.1007/s11032-013-9940-y An energy-dispersive X-ray fluorescence method for analysing Fe and Zn in
Cichy, K. A., Porch, T. G., Beaver, J. S., Cregan, P., Fourie, D., Glahn, R. P., et al. common bean, maize and cowpea biofortification programs. Plant Soil. 419,
(2015). A Phaseolus vulgaris diversity panel for andean bean improvement. 457–466. doi: 10.1007/s11104-017-3352-4
Crop Sci. 55, 2149–2160. doi: 10.2135/cropsci2014.09.0653 Hardner, C. M., Hayes, B. J., Kumar, S., Vanderzande, S., Cai, L., Piaskowski,
Cortes, L. T., Zhang, Z., and Yu, J. (2021). Status and prospects of genome-wide J., et al. (2019). Prediction of genetic value for sweet cherry fruit
association studies in plants. Plant Genome 14, e20077. doi: 10.1002/tpg2.20077 maturity among environments using a 6K SNP array. Hortic. Res. 6, 1–15.
Coyne, D. P., and Schuster, M. L. (1974). Inheritance and linkage relations of doi: 10.1038/s41438-018-0081-7
reaction toXanthomonas phaseoli (E. F. Smith) Dowson (common blight), stage HarvestPlus (2022). BIOFORTIFICATION A food-systems approach to ensuring
of plant development and plant habit in Phaseolus vulgaris L. Euphytica 23, healthy diets globally. Technical brief, HarvestPlus and the World Food
195–204. doi: 10.1007/BF00035858 Programme.
Crossa, J., Perez-Rodriguez, P., Cuevas, J., Montesinos-Lopez, O., Jarquin, Himelblau, E., Mira, H., Lin, S.-J., Cizewski Culotta, V., Peñarrubia, L., and
D., de los Campos, G., et al. (2017). Genomic selection in plant Amasino, R. M. (1998). Identification of a functional homolog of the
breeding: methods, models, and perspectives. Trends Plant Sci. 22, 961–975. yeast copper homeostasis gene ATX1 from Arabidopsis. Plant Physiol. 117,
doi: 10.1016/j.tplants.2017.08.011 1227–1234. doi: 10.1104/pp.117.4.1227
de Campos, T., Oblessuc, P. R., Sforça, D. A., Cardoso, J. M. K., Baroni, R. M., de Huang, M., Liu, X., Zhou, Y., Summers, R. M., and Zhang, Z. (2019).
Sousa, A. C. B., et al. (2011). Inheritance of growth habit detected by genetic BLINK: a package for the next level of genome-wide association studies
linkage analysis using microsatellites in the common bean (Phaseolus vulgaris with both individuals and markers in the millions. Gigascience 8, 1–12.
L.).Mol. Breed. 27, 549–560. doi: 10.1007/s11032-010-9453-x doi: 10.1093/gigascience/giy154
Diaz, L. M., Ricaurte, J., Tovar, E., Cajiao, C., Terán, H., Grajales, M., et al. Izquierdo, P., Astudillo, C., Blair, M. W., Iqbal, A. M., Raatz, B., and Cichy, K.
(2018). QTL analyses for tolerance to abiotic stresses in a common A. (2018). Meta-QTL analysis of seed iron and zinc concentration and content
bean (Phaseolus vulgaris L.) population. PLoS ONE 13, e0202342. in common bean (Phaseolus vulgaris L.). Theor. Appl. Genet. 131, 1645–1658.
doi: 10.1371/journal.pone.0202342 doi: 10.1007/s00122-018-3104-8
Diaz, S., Ariza-Suarez, D., Izquierdo, P., Lobaton, J. D., de la Hoz, J., Acevedo, Jeong, J., Cohu, C., Kerkeb, L., Pilon, M., Connolly, E. L., and Guerinot, M. L.
F., et al. (2019). Replication Data for: Genetic mapping for agronomic traits in (2008). Chloroplast Fe(III) chelate reductase activity is essential for seedling
a MAGIC population of common bean (Phaseolus vulgaris L.) under drought viability under iron limiting conditions. Proc. Natl. Acad. Sci. U.S.A. 105,
conditions’. BMC Genomics. 21, 799. doi: 10.7910/DVN/JR4X4C 10619–10624. doi: 10.1073/pnas.0708367105
Diaz, S., Ariza-Suarez, D., Ramdeen, R., Aparicio, J., Arunachalam, N., Hernandez, Kamfwa, K., Cichy, K. A., and Kelly, J. D. (2015). Genome-wide association
C., et al. (2021). Genetic architecture and genomic prediction of cooking study of agronomic traits in common bean. Plant Genome 8, 1–12.
time in common bean (Phaseolus vulgaris L.). Front. Plant Sci. 11, 622213. doi: 10.3835/plantgenome2014.09.0059
doi: 10.3389/fpls.2020.622213 Kassambara, A., and Mundt, F. (2020). factoextra: extract and visualize the results
Diaz, S., Ariza-Suarez, D., Izquierdo, P., Lobaton, J. D., de la Hoz, J., Acevedo, of multivariate data analyses. Available online at: https://CRAN.R-project.org/
F., et al. (2020). Replication data for: genetic mapping for agronomic traits in package=factoextra
Frontiers in Plant Science | www.frontiersin.org 18 April 2022 | Volume 13 | Article 830896
Keller et al. Improving Genomic Predictions for Beans
Katungi, E., Larochelle, C., Mugabo, J., and Buruchara, R. (2019). Climbing Petry, N., Boy, E., Wirth, J. P., and Hurrell, R. F. (2015). Review: the potential
bean as a solution to increase productivity in land-constrained environments: of the common bean (Phaseolus vulgaris) as a vehicle for iron biofortification.
evidence from Rwanda. Outlook Agric. 48, 28–36. doi: 10.1177/00307270188 Nutrients 7, 1144–1173. doi: 10.3390/nu7021144
13698 Petry, N., Egli, I., Gahutu, J. B., Tugirimana, P. L., Boy, E., and Hurrell, R. (2014).
Keller, B., Ariza-Suarez, D., de la Hoz, J., Aparicio, J. S., Portilla-Benavides, A. Phytic acid concentration influences iron bioavailability from biofortified
E., Buendia, H. F., et al. (2020). Genomic prediction of agronomic traits in beans in Rwandese women with low iron status. J. Nutr. 144, 1681–1687.
common bean (Phaseolus vulgaris L.) under environmental stress. Front. Plant doi: 10.3945/jn.114.192989
Sci. 11, 1001. doi: 10.3389/fpls.2020.01001 Puig, S., Andrés-Colás, N., García–Molina, A., and Peñarrubia, L. (2007).
Kelly, J. D., and Bornowski, N. (2018). “Marker-assisted breeding for economic Copper and iron homeostasis in Arabidopsis: responses to metal deficiencies,
traits in common bean,” in Biotechnologies of Crop Improvement, Volume interactions and biotechnological applications. Plant, Cell Environ. 30,
3: Genomic Approaches, eds S. S. Gosal and S. H. Wani (Cham: Springer 271–290. doi: 10.1111/j.1365-3040.2007.01642.x
International Publishing), 211–238. Raggi, L., Caproni, L., Carboni, A., and Negri, V. (2019). Genome-wide association
Koinange, E. M. K., Singh, S. P., and Gepts, P. (1996). Genetic control of study reveals candidate genes for flowering time variation in common bean
the domestication syndrome in common bean. Crop Sci. 36, 1037–1045. (Phaseolus vulgaris L.). Front. Plant Sci. 10, 962. doi: 10.3389/fpls.2019.00962
doi: 10.2135/cropsci1996.0011183X003600040037x Rehman, H. M., Cooper, J. W., Lam, H.-M., and Yang, S. H. (2019). Legume
Kornegay, J., White, J. W., and de la Cruz, O. O. (1992). Growth habit and gene biofortification is an underexploited strategy for combatting hidden hunger.
pool effects on inheritance of yield in common bean. Euphytica 62, 171–180. Plant Cell Environ. 42, 52–70. doi: 10.1111/pce.13368
doi: 10.1007/BF00041751 Repinski, S. L., Kwak, M., and Gepts, P. (2012). The common bean growth habit
Kwak, M., Toro, O., Debouck, D. G., and Gepts, P. (2012). Multiple origins of the gene PvTFL1y is a functional homolog of Arabidopsis TFL1.Theor. Appl. Genet.
determinate growth habit in domesticated common bean (Phaseolus vulgaris). 124, 1539–1547. doi: 10.1007/s00122-012-1808-8
Ann. Bot. 110, 1573–1580. doi: 10.1093/aob/mcs207 Robinson, N. J., Procter, C. M., Connolly, E. L., and Guerinot, M. L. (1999).
Kwak, M., Velasco, D., and Gepts, P. (2008). Mapping homologous sequences for A ferric-chelate reductase for iron uptake from soils. Nature 397, 694–697.
determinacy and photoperiod sensitivity in common bean (Phaseolus vulgaris). doi: 10.1038/17800
J. Heredity 99, 283–291. doi: 10.1093/jhered/esn.005 Rodriguez, M., Rau, D., Bitocchi, E., Bellucci, E., Biagetti, E., Carboni, A., et al.
Langmead, B., and Salzberg, S. L. (2012). Fast gapped-read alignment with Bowtie (2016). Landscape genetics, adaptive diversity and population structure in
2. Nat. Methods 9, 357–359. doi: 10.1038/nmeth.1923 Phaseolus vulgaris. New Phytol. 209, 1781–1794. doi: 10.1111/nph.13713
Lê, S., Josse, J., and Husson, F. (2008). FactoMineR: an r package for multivariate Rodríguez-Álvarez, M. X., Boer, M. P., van Eeuwijk, F. A., and Eilers, P. H.
analysis. J. Stat. Softw. 25, 1–18. doi: 10.18637/jss.v025.i01 (2018). Correcting for spatial heterogeneity in plant breeding experiments with
Lehermeier, C., Schön, C.-C., and Campos, G. d. l. (2015). Assessment P-splines. Spatial Stat. 23, 52–71. doi: 10.1016/j.spasta.2017.10.003
of genetic heterogeneity in structured plant populations using Ronner, E., Descheemaeker, K., Marinus, W., Almekinders, C. J. M., Ebanyat, P.,
multivariate whole-genome regression models. Genetics 201, 323–337. and Giller, K. E. (2018). How do climbing beans fit in farming systems of the
doi: 10.1534/genetics.115.177394 eastern highlands of Uganda? Understanding opportunities and constraints at
Mangin, B., Siberchicot, A., Nicolas, S., Doligez, A., This, P., and Cierco- farm level. Agric. Syst. 165, 97–110. doi: 10.1016/j.agsy.2018.05.014
Ayrolles, C. (2012). Novel measures of linkage disequilibrium that correct Rosales-Serna, R., Kohashi-Shibata, J., Acosta-Gallegos, J. A., Trejo-López, C.,
the bias due to population structure and relatedness. Heredity 108, 285–291. Ortiz-Cereceres, J., and Kelly, J. D. (2004). Biomass distribution, maturity
doi: 10.1038/hdy.2011.73 acceleration and yield in drought-stressed common bean cultivars. Field Crops
Mayor Duran, V. M. (2016). Identificación de QTLs de frijol común (Phaseolus Res. 85, 203–211. doi: 10.1016/S0378-4290(03)00161-8
vulgaris) asociados a tolerancia a sequía (Ph.D. thesis). Universidad Nacional Sarinelli, J. M., Murphy, J. P., Tyagi, P., Holland, J. B., Johnson, J. W., Mergoum,
de Colombia, Palmira. M., et al. (2019). Training population selection and use of fixed effects to
Mayor Duran, V. M., Raatz, B., and Blair, M. W. (2016). Desarrollo de líneas de optimize genomic predictions in a historical USA winter wheat panel. Theor.
frijol (Phaseolus vulgaris L.) tolerante a sequía a partir de cruces inter acervo Appl. Genet. 132, 1247–1261. doi: 10.1007/s00122-019-03276-6
con genotipos procedentes de diferentes orígenes (Mesoamericano y Andino). Schmutz, J., McClean, P. E., Mamidi, S., Wu, G. A., Cannon, S. B., Grimwood, J.,
Acta Agronómica 65, 431–438. doi: 10.15446/acag.v65n4.48680 et al. (2014). A reference genome for common bean and genome-wide analysis
Meuwissen, T. H., Hayes, B. J., and Goddard, M. E. (2001). Prediction of total of dual domestications. Nat. Genet. 46, 707–713. doi: 10.1038/ng.3008
genetic value using genome-wide dense marker maps.Genetics 157, 1819–1829. Self, S. G., and Liang, K.-Y. (1987). Asymptotic properties of maximum likelihood
doi: 10.1093/genetics/157.4.1819 estimators and likelihood ratio tests under nonstandard conditions. J. Am. Stat.
Miklas, P. N., Kelly, J. D., Beebe, S. E., and Blair, M. W. (2006). Common bean Assoc. 82, 605–610. doi: 10.1080/01621459.1987.10478472
breeding for resistance against biotic and abiotic stresses: from classical to MAS Singh, S. P. (1981). A key for identification of different growth habits of Phaseolus
breeding. Euphytica 147, 105–131. doi: 10.1007/s10681-006-4600-5 vulgaris L. Technical report, International Center for Tropical Agriculture
Möhring, J., and Piepho, H.-P. (2009). Comparison of weighting in (CIAT), Cali, Colombia.
two-stage analysis of plant breeding trials. Crop Sci. 49, 1977–1988. Spindel, J. E., and McCouch, S. R. (2016). When more is better: how data sharing
doi: 10.2135/cropsci2009.02.0083 would accelerate genomic selection of crop plants. New Phytol. 212, 814–826.
Mukamuhirwa, F., and Rurangwa, E. (2018). Evaluation for high iron and zinc doi: 10.1111/nph.14174
content among selected climbing bean genotypes in rwanda. Adv. Crop Sci. Stangoulis, J., and Sison, C. (2008). Crop sampling protocols for micronutrient
Technol. 6, 1–6. doi: 10.4172/2329-8863.1000344 analysis. Harvest Plus Tech. Monogr. Ser. 7, 1–20.
Mukeshimana, G., Butare, L., Cregan, P. B., Blair, M. W., and Kelly, J. Stein, A. J. (2010). Global impacts of human mineral malnutrition. Plant Soil. 335,
D. (2014). Quantitative trait loci associated with drought tolerance 133–154. doi: 10.1007/s11104-009-0228-2
in common bean. Crop Sci. 54, 923–938. doi: 10.2135/cropsci2013. Sul, J. H., Martin, L. S., and Eskin, E. (2018). Population structure in genetic
06.0427 studies: confounding factors and mixed models. PLoS Genet. 14, e1007309.
Nay, M. M., Mukankusi, C. M., Studer, B., and Raatz, B. (2019). Haplotypes doi: 10.1371/journal.pgen.1007309
at the Phg-2 locus are determining pathotype-specificity of angular Teixeira, F. F., Ramalho, M. A. P., and Abreu, N. D. F. B. (1999). Genetic control of
leaf spot resistance in common bean. Front. Plant Sci. 10, 1–11. plant architecture in the common bean (Phaseolus vulgaris L.). Genet Mol Biol.
doi: 10.3389/fpls.2019.01126 22, 577–582. doi: 10.1590/S1415-47571999000400019
Norton, J. (1915). Inheritance of habit in the common bean (Masters Theses). Tello, D., Gil, J., Loaiza, C. D., Riascos, J. J., Cardozo, N., and Duitama, J. (2019).
1911-February 2014. NGSEP3: accurate variant calling across species and sequencing protocols.
Pérez, P., and de los Campos, G. (2014). Genome-wide regression and Bioinformatics 35, 4716–4723. doi: 10.1093/bioinformatics/btz275
prediction with the BGLR statistical package. Genetics 198, 483. Waters, B. M., and Grusak, M. A. (2008). Quantitative trait locus
doi: 10.1534/genetics.114.164442 mapping for seed mineral concentrations in two Arabidopsis thaliana
Frontiers in Plant Science | www.frontiersin.org 19 April 2022 | Volume 13 | Article 830896
Keller et al. Improving Genomic Predictions for Beans
recombinant inbred populations. New Phytol. 179, 1033–1047. Yu, G., Smith, D. K., Zhu, H., Guan, Y., and Lam, T. T.-Y. (2017). ggtree:
doi: 10.1111/j.1469-8137.2008.02544.x an r package for visualization and annotation of phylogenetic trees with
Weller, J. L., Vander Schoor, J. K., Perez-Wright, E. C., Hecht, V., González, their covariates and other associated data. Methods Ecol. Evolut. 8, 28–36.
A. M., Capel, C., et al. (2019). Parallel origins of photoperiod adaptation doi: 10.1111/2041-210X.12628
following dual domestications of common bean. J. Exp. Bot. 70, 1209–1219.
doi: 10.1093/jxb/ery455 Conflict of Interest: The authors declare that the research was conducted in the
White, J., Kornegay, J., Castillo, J., Molano, C., Cajiao, C., and Tejada, G. absence of any commercial or financial relationships that could be construed as a
(1992). Effect of growth habit on yield of large-seeded bush cultivars of potential conflict of interest.
common bean. Field Crops Res. 29, 151–161. doi: 10.1016/0378-4290(92)9
0084-M Publisher’s Note: All claims expressed in this article are solely those of the authors
White, J. W., and Izquierdo, J. (1989). “Dry bean: physiology of yield potential and and do not necessarily represent those of their affiliated organizations, or those of
stress tolerance,” in Common Beans: Research for Crop Improvement (Santiago: the publisher, the editors and the reviewers. Any product that may be evaluated in
CAB International).
this article, or claim that may be made by its manufacturer, is not guaranteed or
Wickham, H. (2016). “Toolbox,” in ggplot2: Elegant Graphics for Data
Analysis, Use R!, ed H. Wickham (Cham: Springer International Publishing), endorsed by the publisher.
33–74.
Wu, H., Li, L., Du, J., Yuan, Y., Cheng, X., and Ling, H.-Q. (2005). Molecular Copyright © 2022 Keller, Ariza-Suarez, Portilla-Benavides, Buendia, Aparicio,
and biochemical characterization of the Fe(III) chelate reductase gene family Amongi, Mbiu, Msolla, Miklas, Porch, Burridge, Mukankusi, Studer and Raatz.
in Arabidopsis thaliana. Plant Cell Physiol. 46, 1505–1514. doi: 10.1093/pcp/ This is an open-access article distributed under the terms of the Creative Commons
pci163 Attribution License (CC BY). The use, distribution or reproduction in other forums
Wu, J., Wang, L., Fu, J., Chen, J., Wei, S., Zhang, S., et al. (2020). is permitted, provided the original author(s) and the copyright owner(s) are credited
Resequencing of 683 common bean genotypes identifies yield component and that the original publication in this journal is cited, in accordance with accepted
trait associations across a north–south cline. Nat. Genet. 52, 118–125. academic practice. No use, distribution or reproduction is permitted which does not
doi: 10.1038/s41588-019-0546-0 comply with these terms.
Frontiers in Plant Science | www.frontiersin.org 20 April 2022 | Volume 13 | Article 830896