Multienvironment trials (METs) are widely used to assess the performance of promising crop germplasm. Though seldom designed to elucidate genetic mechanisms, MET data sets are often much larger than could be duplicated for genetic research and, given proper interpretation, may offer valuable insights into the genetics of adaptation across time and space. The Cooperative Dry Bean Nursery (CDBN) is a MET for common bean (Phaseolus vulgaris) grown for > 70 years in the United States and Canada, consisting of 20–50 entries each year at 10–20 locations. The CDBN provides a rich source of phenotypic data across entries, years, and locations that is amenable to genetic analysis. To study stable genetic effects segregating in this MET, we conducted genome-wide association studies (GWAS) using best linear unbiased predictions derived across years and locations for 21 CDBN phenotypes and genotypic data (1.2 million SNPs) for 327 CDBN genotypes. The value of this approach was confirmed by the discovery of three candidate genes and genomic regions previously identified in balanced GWAS. Multivariate adaptive shrinkage (mash) analysis, which increased our power to detect significant correlated effects, found significant effects for all phenotypes. Mash found two large genomic regions with effects on multiple phenotypes, supporting a hypothesis of pleiotropic or linked effects that were likely selected on in pursuit of a crop ideotype. Overall, our results demonstrate that statistical genomics approaches can be used on MET phenotypic data to discover significant genetic effects and to define genomic regions associated with crop improvement.
Almost every crop improvement program assesses the performance of promising germplasm and breeding material via multienvironment trials (METs). The phenotypic data produced by these trials are extremely important guides to growers, private seed companies, and public institutions involved in crop improvement, because combining trial data from multiple years and locations increases the probability of identifying genotypes that perform well, or show especially desirable traits (Bowman 1998). Many cooperative testing networks conduct METs to enable cooperators and other interested parties to observe performance over a wider range of environments than if they were only tested locally (Annicchiarico 2002). This supports the identification of advanced lines with stable, high performance in multiple production environments. Among many others, crop testing networks that conduct METs include the U.S. cooperative regional performance testing program, the University Crop Testing Alliance, and the Cooperative Dry Bean Nursery (CDBN) (Singh 2000).
Longstanding METs such as the CDBN have often focused on breeding for crop ideotypes, in addition to breeding to eliminate defects and to select for yield. Donald (1968) defined a crop ideotype as an idealized plant with trait combinations expected to produce a greater yield quantity or quality. In contrast, approaches that eliminate defects or select for yield do not consider desirable combinations of traits; thus, these approaches only produce desirable combinations by chance. Selection for an ideotype involves selection for correlated traits, and could lead to substantial pleiotropy, where a single gene affects multiple traits. METs like the CDBN, which were used to select for specific crop ideotypes, could provide insight into the genetics of trait correlations in crop genomes.
Though METs are often used to measure genetic gain over time (Graybosch and Peterson 2010; Vandemark et al. 2014), the vast majority of METs are designed to measure phenotypic responses to a broad set of targeted growing environments. The experimental designs of METs can pose substantial analytical challenges to additional, unplanned genetic analyses. METs typically produce sparse data matrices of phenotypes across germplasm entries, locations, and years (Figure 1). The frequency of different germplasm entries may vary as part of the normal selection process. Thus, entries with good performance are often tested in more locations and years than those with poor performance. With the exception of few standard checks, the set of genotypes tested each year typically varies, with most genotypes tested in only 1 or 2 years. In addition, the total number of genotypes tested each year can vary substantially, and this number is typically too small for genome-wide association on any 1 year’s data alone. Over the years, MET cooperators can also join or leave the network, and add or drop MET sites or phenotypes due to changes in research focus, personnel, or funding. All of these variations make METs into large unbalanced data sets that need to be handled properly for genetic work. Genetic analyses of MET germplasm can also be hampered by the difficulty of obtaining and genotyping previously evaluated entries, particularly entries with poor trial performance that were not tested further. This difficulty may bias or prevent studies that require genetic diversity to explain phenotypic variation, such as genome-wide association studies (GWAS). In contrast, field experiments designed for genetic studies assess complete, balanced designs, and produce data matrices of phenotypes across genotypes and environments with few or no missing cells. Ideally, the number of genotypes is identical across all environments, and a minimum of a few hundred genotypes are tested in each environment. Each genotype is also tested an equivalent number of times across sites and years.
Despite these analytical issues, METs often produce decades of phenotypic data, which gives them substantial appeal for use in genetic analyses of phenotypic variation. Genetic analyses of MET data sets have recently been implemented in several crop species (Hamblin et al. 2010; Rife et al. 2018; Sukumaran et al. 2018). Common bean has nutritional and agronomic importance, a long history of METs, and many emerging genomic tools, making it an outstanding species in which to assess METs that might support the genetic analysis of phenotypic variation. Common bean is the most consumed plant protein source worldwide and is a particularly important source of protein in the developing world (Jones 1999). In North America, common bean improvement efforts remain mostly in the public sector, and over the past 70 years, the CDBN has been a major testing platform for these improvement efforts. The CDBN is the largest MET for common bean in the United States and Canada (Myers 1988; Singh 2000), and CDBN cooperators have collected phenotypic data on > 150 traits for hundreds of advanced breeding lines and released cultivars (hereafter entries) of common bean at > 70 locations (Figure 1), which produced up to 18,000 recorded data points per trait (Figure 2A). The traits are of economic and/or agronomic importance to bean producers, and include seed yield, growth habit, seed size, phenology, and disease responses, among others (Figure 2A, Supplemental Material, Figure S1).
More than 500 CDBN entries have been grown since the 1980s (Figure 1). These entries include released cultivars and unreleased advanced breeding lines, which represent most bean types grown in North America. Entries include > 13 market classes of common bean that group into three major races from two independent domestication events (Mamidi et al. 2011) (Figure 1). Therefore, the CDBN can be used as a representative sample of the genetic diversity being used by North American bean breeders in their programs throughout the last 70 years. However, phenotypic data from the CDBN are sparse and unevenly distributed: the average CDBN entry was grown at only 19 of the 70 locations and in 2 of the 34 years, with substantial variation in these numbers. CDBN cooperators grew between 16 and 61 of the > 500 entries each year, and used 10–28 of the > 70 locations per year (Figure 1). Individual CDBN locations grew between 8 and 514 entries, with a median of 74 entries. Locations were used in the CDBN for as few as 1 to as many as 34 years, with a median of 5 years of participation. Though genotypes are present only intermittently over CDBN locations and years, the vast phenotyping effort on this interrelated set of bean germplasm, when combined with genomic data, offers an excellent opportunity to identify genomic regions affecting phenotypic variation in this species.GWAS have elucidated candidate genes and genomic regions that affect trait variation in many other crop species (Atwell et al. 2010; Kirby et al. 2010; Mackay et al. 2012; Lin et al. 2014; McCouch et al. 2016; MacArthur et al. 2017; Xiao et al. 2017; Togninalli et al. 2018), and have recently been implemented in common bean (Cichy et al. 2015; Kamfwa et al. 2015a,b; Moghaddam et al. 2016; Soltani et al. 2017, 2018; Tock et al. 2017; Nascimento et al. 2018; Oladzad et al. 2019a,b; Raggi et al. 2019). Combining sparse phenotypic data in agricultural data sets to look for pleiotropic or linked effects across conditions has parallels in human biomedical GWAS. In these trials, individual clinics can assess only a subset of human genotypes, and patients are evaluated using institution-specific criteria (Lotta et al. 2017; Visscher et al. 2017). Human GWAS often look for common variants for common diseases and correct phenotypes for effects of age, sex, and location (Schork et al. 2009; Mefford and Witte 2012; Zaitlen et al. 2012). Analogously, we seek common, genetically stable variants for important phenotypes evaluated in a MET, corrected for effects of location, year, kinship, and assessment criteria. In human biomedical GWAS, pleiotropic effects of SNPs on multiple diseases have frequently been observed (Sivakumaran et al. 2011). Selection for a common bean crop ideotype—with a long hypocotyl, many nodes carrying long pods, and without side branches, small leaves, and determinate growth (Adams 1982; Kelly 2001)—is known to have led to pleiotropic or linked effects on multiple traits, such as seed yield, biomass, lodging, and plant height (Soltani et al. 2016). To study the genetic effects of this aspect of the CDBN selection framework, we used multivariate adaptive shrinkage (mash) to find genomic associations with significant effects on one or more CDBN phenotype (Urbut et al. 2019). Mash is a flexible, data-driven method that shares information on patterns of effect size and sign in any data set where effects can be estimated on a condition-by-condition basis for many conditions (here, phenotypes) across many units (here, SNPs). It first learns patterns of covariance between SNPs and phenotypes from SNPs without strong effects, then combines these data-driven covariances with the original condition-by-condition results to produce improved effect estimates. In this way, mash shares information between conditions to increase the power to detect shared patterns of effects. Mash was originally used for analyses of human biomedical data (Urbut et al. 2019) and has yet to be used in an agricultural setting. This analysis method could be used with the rich phenotypic resources of crop METs to understand genetic effects across multiple phenotypes or across multiple locations and years.
Here, we demonstrate that the CDBN MET data set can be used to make genetic discoveries, despite the sparse nature of the data, by using best linear unbiased predictions (BLUPs) for entries phenotyped in the CDBN. We explore whether this approach can find genomic regions significantly associated with phenotypic variation, and compare associations found with this approach to published GWAS results obtained from more balanced trials. We also explore patterns of genomic associations with significant effects on more than one CDBN phenotype using mash. Our results demonstrate the value of adding a genetic component to data sets such as the CDBN and provide a starting point for future work that explores the genetics of phenotypes evaluated in METs.
MET data sets represent substantial phenotypic resources that can aid in the genetic study of important agronomic phenotypes. Several important steps in preparing the CDBN data for analysis fall under the remit of data science, and specifically involve the data processing steps outlined here. First, when available only from printed reports, the data were rendered machine-readable. Processing of the digitized data next involved cleaning the data to remove inconsistencies and spurious data, then filtering to retain only the relevant data. The data were stored in a consistent form where the semantics of the data set matched the way it was stored. Then, various data scales for individual traits such as growth habit were standardized to create phenotypes that were more consistent across locations and years. The phenotypic data were next enriched with additional attributes that made subsequent analyses more meaningful, such as germplasm, environment, and crop management information. Then, the data were aggregated to create summary data by estimating BLUPs for each phenotype. We next used a GWAS modeling approach to determine the genomic regions associated with these data summaries. Finally, we used mash to examine the patterns of overlap between genomic associations with significant effects on one or more phenotype (Urbut et al. 2019).
Phenotypic data for entries grown in the CDBN were available mainly as hard-copy reports providing plot averages at named locations. Some reports were available in the National Agricultural Library from the 1950s onward; however, reports from 1981 onward had substantial additional available genetic material and were the focus for this analysis (Table S1). Reports from 1981 to 2015 were scanned if not in digital format, digitized using optical image recognition as required, and then reformatted using custom SAS (SAS System, version 9.4, SAS Institute, Cary, NC) scripts that also standardized nomenclature and units of measurement.
Much of the phenotypic data required additional processing to allow comparisons across locations and years. The long time span and large number of testing locations led to the scoring of 152 traits. Many of these traits represented distinct methods for scoring similar phenotypes; for example, lodging was scored on a percent scale, a 1–5 scale, a 0–9 scale, and a 1–9 scale at different locations and in different years; for this analysis, these lodging traits were standardized to one lodging phenotype on a 1–5 scale. From 152 traits reported, 22 phenotypes were standardized for use in GWAS, including 8 quantitative phenotypes and 14 qualitative phenotypes created from visual scores and/or specific measurements (Figure 2A). The output from the R script used to standardize the phenotypes across locations and years can be found online at http://rpubs.com/alice_macqueen/CDBN_Phenotype_Standardization.
We generated phenotypes associated with location code, year, and genotype information. A total of 70 location codes were created as four-letter abbreviations with the U.S. state or Canadian province abbreviation as the first two letters, and the specific site abbreviation as the second two letters. Five location codes ending in “2” corresponded to a second trial grown at that location and year, usually with a treatment such as drought or disease applied. Location codes were associated with latitude, longitude, elevation, and other location-specific metadata (Table S2), while genotypes were associated with market class and race, as well as the availability of seed from the holdings of CDBN cooperators and SNP data, where available (Table S2).
In general, location by year (L*Y) combinations with outlier phenotypic values [values above the third quartile or below the first quartile by 1.5 times the interquartile range (IQR)] were removed for every entry in that L*Y combination. Removing outlier L*Y combinations prevented possible bias from linear models using a biased sample of data points for a L*Y, while still removing points that, by IQR measures and by knowledge of reasonable ranges for common bean quantitative phenotypes, were likely due to mismeasurement or data entry errors. The specifics of phenotype standardization for all 22 phenotypes are given in the Supplementary Note and the code is available on GitHub at https://github.com/Alice-MacQueen/CDBNgenomics/tree/master/analysis-paper.
Selection and breeding strategies to generate new bean entries for the CDBN varied across years and among breeding programs. However, in general, new advanced lines were selected from either single, triple, or double crosses among advanced breeding material and released cultivars, which in most cases were already tested within the CDBN in previous years. These lines were bulked to increase seed supply, then field tested to ensure consistency of phenotypic responses in the advanced lines. Entries with favorable characteristics were often entered into the CDBN to be phenotyped in multiple environments. Consequently, most CDBN entries were members of a complex pedigree that has had novel, favorable alleles recombined or introgressed into it over time.
It is clear that the CDBN is not a randomly mating, homogeneous population, and the breeding and selection strategy in the CDBN likely impacts GWAS on this material in a number of ways. Presumably, breeders have increased the frequency of alleles that favorably affect phenotypes over time, which should aid in the detection of these genomic regions via GWAS. The multiple generations of inbreeding should reduce allelic heterogeneity, which should also aid GWAS. Indeed, we find few heterozygous regions in our SNP data set and few examples of multiallelic loci. By the same token, the frequent inbreeding may also increase the size of linkage disequilibrium (LD) blocks or cause spurious patterns of LD, which may cause nonsyntenic associations and make candidate gene identification more difficult. In addition, the infrequent crosses between the gene pools from the two independent domestication events, and the assortative mating practiced as part of the breeding strategy, could lead to an inflated false-positive rate and create correlations between previously uncorrelated traits (Li et al. 2017).
To detect genomic regions associated with phenotypic variation in a GWAS framework, it is particularly valuable to have a large amount of heritable phenotypic variation. Thus, it was equally important to include entries from the CDBN with poor seed yields or nonideal phenotypic traits as high-yielding, commercially released varieties. Thus, we went to considerable effort to obtain seed of unreleased, unarchived materials from the holdings of CDBN cooperators. Germplasm from the entries grown in the CDBN was obtained from multiple sources, including the International Center for Tropical Agriculture, the National Plant Germplasm System, and three common bean diversity panels: the Mesoamerican Diversity Panel (MDP) (Moghaddam et al. 2016), the Durango Diversity Panel (DDP) (Soltani et al. 2016), and the Andean Diversity Panel (ADP) (Cichy et al. 2015). Seed was also obtained from holdings of CDBN cooperators, including Mark Brick (Colorado State University), Jim Kelly (Michigan State University), Phil McClean (North Dakota State University), Phil Miklas (U.S. Department of Agriculture–Agricultural Research Service), James Myers (Oregon State University), Juan Osorno (North Dakota State University), and Tom Smith (University of Guelph).
The SNP data set was created from this germplasm in two ways. First, raw sequence data were obtained from the ADP, DDP, and MDP (Cichy et al. 2015; Moghaddam et al. 2016) for CDBN entries and all parents of CDBN entries that had been sequenced as part of these panels. The remainder of the CDBN was genotyped using identical methodology to these previous diversity panels, dual-enzyme genotyping-by-sequencing (Schröder et al. 2016). Unfortunately, 39 of the older, unreleased varieties would no longer germinate. For these varieties, we obtained DNA for sequencing by rehydrating sterilized seeds on wetted Whatman paper in petri plates for 2–3 days, then dissecting the embryo from the seed and extracting DNA from the embryo. The DNA from the remaining entries was extracted from young trifoliates. The enzymes MseI and TaqI were used for digestion following the protocol from Schröder et al. (2016). SNPs were called from this raw sequence data using the pipeline found at https://github.com/Alice-MacQueen/SNP-calling-pipeline-GBS-ApeKI. Briefly, cutadapt was used to trim adapters and barcodes (Martin 2011), sickle adaptive trimming was used to remove ends of reads with quality scores < 20 (Joshi and Fass 2011), bwa mem was used to align reads to V2.0 of the G19833 reference genome found at https://phytozome.jgi.doe.gov/pz/portal.html#!info?alias=Org_Pvulgaris (Li and Durbin 2010; Schmutz et al. 2014), and NGSEP (Next Generation Sequencing Experience Platform) was used to call SNPs for the entire set of CDBN entries and all parents in the CDBN pedigrees (Duitama et al. 2014). SNPs were imputed using FILLIN in TASSEL. This resulted in the creation of a diversity panel of 327 entries with MET data in the CDBN (Table S2), with aligned SNP data available on the University of Texas Libraries data repository at DOI: https://doi.org/10.18738/T8/RTBTIR for use in the CDBNgenomics R package at https://github.com/Alice-MacQueen/CDBNgenomics.
To explore consistent genetic effects that could be compared to balanced genetic trials, analyses were performed on genetic BLUPs for each phenotype. BLUPs were calculated in the rrBLUP package in R using a kinship matrix, and treating location and the interaction between location and year as fixed effects. The R code to generate the BLUPs is available on GitHub at https://github.com/Alice-MacQueen/CDBNgenomics/tree/master/data-raw. The BLUPs are available in Table S2. R code treating CDBN germplasm entry as a fixed effect, and treating location, the interaction between location and year, the interaction between CDBN germplasm entry and location, and the interaction between CDBN germplasm entry and year as random effects is also available on GitHub at https://github.com/Alice-MacQueen/CDBNgenomics/tree/master/data-raw/ASReml; this model did not significantly affect association rank order or subsequent analyses. rrBLUP was also used to calculate narrow-sense heritability (h2), defined as Va / (Va + Ve), where Va is the restricted maximum likelihood (REML) estimate of the additive genetic variance from rrBLUP and Ve is the REML estimate of the error variance.
For GWAS phenotypes, BLUPs were retained only for CDBN entries phenotyped at least one time in the CDBN. The kinship matrix was calculated using default methods in GAPIT. A total of 1,221,540 SNPs with a minor allele frequency > 5% in the CDBN diversity panel were identified and used for the CDBN GWAS. GWAS analyses were performed using compressed mixed linear models (Zhang et al. 2010) implemented in GAPIT with the optimum level of compression (Lipka et al. 2012). These models used a kinship matrix calculated within GAPIT to account for individual relatedness, and some number of principal components (PCs) to account for population structure. The optimum number of PCs to account for population structure was determined using model selection in GAPIT, and by selecting the number of PCs that maximized the Bayesian information criterion. Typically, zero to two PCs were used (Table S3). The final Manhattan plots were created using the ggman R package. Plots of intersecting sets were created using the UpSetR package (Lex et al. 2014). Candidate genes within a 20-kb interval centered on the peak SNP with P-values above a Benjamini–Hochberg false discovery rate (FDR) threshold of 0.1 were examined further. In a balanced GWAS panel that had many identical genotypes, this interval size was suitable for capturing LD in the data (Moghaddam et al. 2016).
Out of the 21 BLUPs estimated from CDBN phenotypes, a group of 13 also had published associations from GWAS on common bean. To compare the major associations in our study to those of published studies on balanced genetic trials, we collected the major associations reported in 11 published GWAS studies of common bean (Cichy et al. 2015; Kamfwa et al. 2015a,b; Moghaddam et al. 2016; Soltani et al. 2017, 2018; Tock et al. 2017; Nascimento et al. 2018; Oladzad et al. 2019a,b; Raggi et al. 2019). We compared these published associations to the associations for the top 10 SNPs for each of the 13 phenotypes in this study, thinned to one SNP per 20-kb region. Unfortunately, these comparisons were likely very conservative, in that most of these publications used panels of common bean that were comprised of material from different gene pools than the CDBN, with the exception of the MDP and DDP (Moghaddam et al. 2016; Soltani et al. 2016; Oladzad et al. 2019a,b). Both Andean and Middle-American gene pools have been observed to have different SNPs underlying domestication traits (Schmutz et al. 2014). Eight of these publications used v1.0 of the Phaseolus vulgaris genome annotation, while our associations were mapped to v2.0. We used the genome browser located at https://legumeinfo.org/genomes/gbrowse/phavu.G19833.gnm2 to convert associations between these two versions of the genome annotation. We then determined the number of overlapping associations meeting two criteria: first, those within 200 kb of one another, and second, within 20 kb of one another and with the same candidate gene. We determined these overlaps for the 80 associations from the 11 published GWAS to find an expected rate of overlap, then compared this to the rate of overlap between this study and the 11 balanced GWAS.
To increase our power to detect associations above an FDR, and to find genomic associations with significant effects on one or more CDBN phenotype, we used a two-step empirical Bayes procedure, mash, to estimate effects of ∼45,000 SNPs on 20 BLUPs determined from CDBN phenotypes (Urbut et al. 2019). Mash has been demonstrated to increase power to detect effects in analyses of human data (Urbut et al. 2019), and while the methods are extensible to any data set with many SNPs/markers and many phenotypes/conditions, it has not yet been used in an agricultural setting. It can be used for any data set where effects can be estimated for many conditions (here, phenotypes) across many units (here, SNPs). Mash first learns patterns of covariance between SNPs and phenotypes from SNPs without strong effects, then combines these data-driven covariances with the original condition-by-condition results to produce improved effect estimates. In this way, mash shares information between conditions to increase the power to detect shared patterns of effects. Importantly, this method does not have restrictive assumptions about the patterns of effects between markers or conditions. In addition, estimates with little uncertainty are not adversely affected by the inclusion of estimates with high uncertainty. Thus, we included 20 phenotypes in the mash analysis, including 12 phenotypes with no signal above the Benjamini–Hochberg FDR threshold in individual GWAS. Two low-signal phenotypes related to bean common mosaic virus presence or absence were not included; inclusion of these phenotypes did not significantly alter the mash results (data not shown). The procedure we used to generate input matrices for mash is captured in the R package gapit2mashr, available at https://github.com/Alice-MacQueen/gapit2mashr. Briefly, the effect of the alternate allele relative to the reference allele was determined for each SNP using GAPIT. To allow mash to converge effectively on effect estimates, the effects for each phenotype were standardized to fall between −1 and 1, with a mean of 0. Because mash does not accept NA values, when GAPIT calculated SE for ≤ 95% of the SNPs in the GWAS, we instead calculated the SE for that phenotype using Hedges’ G (Hedges and Olkin 1985).
Data-driven covariance matrices were estimated using 45,000 randomly selected SNPs from the entire set of 1,221,540 SNPs. These matrices were then used on the top 4000 SNPs for each of the 20 traits, as determined by P-values in the individual GWAS, which produced a matrix of strong effects for 45,000 SNPs. Then, we explored the patterns of significant effects in the mash output. We first determined which SNPs had evidence of significant phenotypic effects by determining SNPs with the largest Bayes factors. In this analysis, the Bayes factor was the ratio of the likelihood of one or more significant phenotypic effects at a SNP to the likelihood that the SNP had only null effects. Here, following Kass and Raftery (1995), a Bayes factor of > 102 is considered decisive evidence in favor of the hypothesis that a SNP has one or more significant phenotypic effects. We also compared the size of significant phenotypic effects, as determined by SNPs with a local false-sign rate of ≤ 0.05 for one or more phenotype. The local false-sign rate is analogous to an FDR, but is more conservative, in that it also reflects the uncertainty in the estimation of the sign of the effect (Stephens 2017).
Genotypic data are available at the Sequence Read Archive (SRA) under submission number SUB6162710. Code for SNP calling is available at https://github.com/Alice-MacQueen/SNP-calling-pipeline-GBS-ApeKI.
Aligned SNP data are available at https://doi.org/10.18738/T8/RTBTIR. Raw phenotypic data are available in the National Agricultural Library: https://www.nal.usda.gov/. Code used to generate data used in this analysis from the raw phenotypic data are available at Rpubs, found at: http://rpubs.com/alice_macqueen/CDBN_Phenotype_Standardization. Code and data necessary to replicate this analysis are available as part of the R package CDBNgenomics, found at: https://github.com/Alice-MacQueen/CDBNgenomics. Supplemental material available at Texas Data Repository Dataverse: https://doi.org/10.18738/T8/KZFZ6K.
The CDBN contains a wealth of data to study the genetics of phenotypes and phenotypic correlations (Figure 1 and Figure 2A). We were able to obtain and genotype 327 germplasm entries from the > 544 entries present in the CDBN trials from 1981 to 2015, including 124 entries that were neither released commercially nor submitted to the National Plant Germplasm System (National Plant Germplasm System 2017), and 39 entries whose seed would not germinate. Most of the remaining entries were grown in the CDBN before 1990 and had seed stocks that, for reasons of practicality, were no longer maintained by breeders (Figure 1). Genotyping-by-sequencing of the available genotypes generated 1.2 million SNPs for analysis of stable effects in the CDBN.
BLUPs of phenotypes from the CDBN, conditioned on location, location by year, and the kinship matrix, are analogous to breeding values for the CDBN entries. These genetic values can be used to determine the narrow-sense heritability, h2, potentially explainable by GWAS. h2 varied between 6 and 73% in the 21 phenotypic BLUPs (Table 1), reflecting the substantial amount of environmental variation in the data set. We then determined the correlations between the BLUPs of CDBN phenotypes or the genetic correlations. Correlation coefficients between BLUPs of CDBN phenotypes varied between −0.75 and 0.81, and most phenotypes were significantly correlated (Figure S1). Two major groups of phenotypes were positively correlated: biomass, days to flowering, plant height, zinc deficiency score, days to maturity, blackroot presence/absence, and early vigor were in the first of these groups, and white mold damage score, growth habit, seed yield, harvest index, lodging, rust damage score, bean common mosaic virus damage score, and halo blight damage score were in the second of these groups. These two groups had negative phenotypic correlations with each other.
|Seed yield||53,173||222,409||0.193||1,727||kg ha−1|
|Days to maturity||12.3||18.3||0.402||17.0||Days|
|Days to flowering||4.68||6.11||0.434||11.7||Days|
|Seed fill duration||5.584||17.812||0.239||12.7||Days|
|Lodging score||0.244||0.466||0.344||2.97||1–5 scale|
|Growth habit||0.160||0.154||0.509||2.19||1–3 scale|
|Seed appearance score||0.017||0.241||0.067||0.346||1–3 scale|
|CBB damage score||0.282||1.127||0.200||2.32||1–9 scale|
|Rust damage score||3.201||1.957||0.621||7.35||1–9 scale|
|Early vigor score||0.064||0.747||0.078||0.957||1–9 scale|
|White mold damage score||0.143||0.631||0.185||2.60||1–5 scale|
|CTV presence/absence||0.025||0.132||0.157||0.432||0–1 scale|
|Halo blight damage score||0.176||0.708||0.199||1.12||1–5 scale|
|BCMV blackroot response||0.049||0.086||0.363||0.843||0–1 scale|
|BCMV presence/absence||0.016||0.097||0.145||0.428||0–1 scale|
|Root rot damage score||0.455||3.193||0.125||2.25||1–9 scale|
|Zinc deficiency damage score||3.831||1.395||0.733||8.96||1–9 scale|
We conducted GWAS on 21 phenotypes using BLUPs calculated using a kinship matrix, location, and an interaction between location and year as fixed effects (for details, see the GWAS section in the Materials and Methods). To determine if any SNP frequencies had changed over the duration of the CDBN, we also conducted GWAS on the earliest year that each germplasm entry was present in the CDBN as a proxy for the age of the entry. This GWAS was analogous to an environmental GWAS that uses climatic variables associated with a genotype’s location of origin (Hancock et al. 2011; Lasky et al. 2015), though this GWAS is fitted to a variable correlated with the age of the genotype rather than with its location of origin.
Given the analytical issues surrounding the use of METs for unplanned genetic analyses, it was unclear whether GWAS on CDBN phenotypes would find significant associations, or if these associations would be reduced or eliminated by environmental noise, or by experimental design biases. Thus, we determined if any GWAS on CDBN phenotypes had significant associations after a Benjamini–Hochberg FDR correction of 10%. With this criterion, significant associations were discovered for 8 of the 21 phenotypes. More than 33 peaks had SNPs with P-values above the FDR, indicating the presence of ≥ 30 distinct, significant associations with these eight CDBN-derived phenotypes. Phenotypes with associations above the FDR generally had more data points in the CDBN (6500 vs. 2400 data points, Wilcoxon rank sum test P = 0.018; Figure 2A). Phenotypes with associations above the FDR also had significantly higher narrow-sense heritabilities estimated from the phenotypic data (h2 of 40.5 vs. 25%, Wilcoxon rank sum test P = 0.038, Table 1). We briefly discuss the associations above the FDR for these eight phenotypes in the order of most to fewest data points in the CDBN. In cases where there were multiple associations for a single phenotype, we discuss only the top associations by P-value.
Seed yield (kg ha−1) had one significant peak after FDR correction, on Pv01 at 42.2 Mb (Figure 2, B and C and Table S4). This association was correlated with a difference in seed yield of 104 kg ha−1 (Figure 2F and Table S4). Median seed yields in the CDBN for the Durango, Mesoamerican, and Nueva Granada races were 2803, 2443, and 2038 kg ha−1, respectively; thus, this genomic region accounts for changes in seed yield of 3.7–5.1%, or 3–4 years of improvement effort at historical rates of bean improvement (Vandemark et al. 2014). This association was 3.7-kb upstream of the gene Phvul.001G167200, a gene that is highly expressed in the shoot and root tips of common bean at the second trifoliate stage of development (O’Rourke et al. 2014; Dash et al. 2016). The Arabidopsis thaliana homolog of this gene, VERNALIZATION INDEPENDENCE 5 (VIP5), affects flowering time by activating Flower Locus C, which is a repressor of flowering (Oh et al. 2004).
Seed weight (mg) had associations on nine chromosomes that were significant after FDR (Figure 2, D and E); the strongest of these were on Pv02 (Figure 2G), Pv03, Pv05, and Pv08, though each explained only 1–2% of the variation in seed weight (Table S4). Because seed weight correlates strongly with population structure in the three bean races and two bean gene pools, seven PCs were used to correct for population structure in this GWAS (Table S3). The association on Pv02 was 5-kb upstream of gene model Phvul.002G150600, a Sel1 repeat protein. Sel1-like repeat proteins are frequently involved in signal transduction pathways and in the assembly of macromolecular complexes (Mittl and Schneider-Brachert 2007). The association on Pv03 was 10-kb upstream of gene model Phvul.003G039900, a jasmonic acid carboxyl methyltransferase. The association on Pv05 was not within 20 kb of any gene. The association on Pv08 fell in the coding sequence of Phvul.008G290600, a choline-phosphate cytidylyltransferase highly expressed in many tissues, including roots and pods, and seeds at the heart stage and stage 2, or seeds 3–4- and 8–10-mm wide (O’Rourke et al. 2014; Dash et al. 2016).
Days to flowering had one significant peak after FDR, on Pv01 between 13.4 and 17.1 Mb (Figure S2A and Table S4). It was correlated with a difference in flowering time of 2–3 days, depending on the population (Figure S3A). A candidate gene model hypothesized to affect days to flowering, Phvul.001G087500, is located at 13.76 Mb in the V2.0 annotation for P. vulgaris. Gene model Phvul.001G087500 is an ortholog of KNUCKLES (KNU), a protein that is part of the Polycomb repressive complex 2, a complex that affects both flowering time and floral meristem development (de Lucas et al. 2016). KNU is activated in the transition to determinate floral meristem development and functions in a feedback loop that promotes determinate development (Payne et al. 2004; Sun et al. 2014).
Lodging score, where higher scores indicated more stem breakage near ground level, had associations on three chromosomes that were significant after FDR; one on Pv04 at 2.8 Mb, one on Pv05 at 0.4 Mb, and one on Pv07 at 34.5 Mb (Figure S2B and Table S4). In total, these three associations explained 8% of the variation in lodging (Figure S4). The signal on Pv04 fell within gene model Phvul.004G025600; the A. thaliana homolog of this gene is involved in the biosynthesis of inositol pyrophosphate, a cellular signaling molecule involved in metabolism and energy sensing (Desai et al. 2014). The signal on Pv05 fell within gene model Phvul.005G005400, a uridine diphosphate glycosyltransferase superfamily protein (Dash et al. 2016). The strongest signal for lodging, explaining 3% of the variation, fell in the promoter region of gene model Phvul.007G221800, which is orthologous to SUPPRESSOR OF AUXIN RESISTANCE 1 (SAR1). In A. thaliana, SAR1 increases plant height and internode distance, and appears to affect stem thickness (Cernac et al. 1997; Parry et al. 2006).
Harvest index, or the ratio of seed yield weight to total above ground biomass, had one significant association on Pv03 at 2.1 Mb (Figure S2C and Table S4). The alternate allele was associated with an increase in harvest index of 1.5–3.5%, and associated with bean race (Figure S3B). This allele was 20 kb from gene model Phvul.003G023000, a cellulose synthase-like protein highly expressed in green mature pods, whole roots, and leaf tissue at the second trifoliate leaf stage of development (O’Rourke et al. 2014; Dash et al. 2016).
Growth habit encompasses both determinate and indeterminate types (I and II/III), as well as upright and prostrate indeterminate types (II and III). Growth habit had significant associations on every chromosome after FDR; the strongest four associations were on Pv01 at 6.2 and 42.2 Mb, on Pv09 at 30.9 Mb, and Pv10 at 42.7 Mb (Figure S2D and Table S4). There are known to be multiple determinacy loci segregating in different gene pools of common bean (Kolkman and Kelly 2003; Kwak et al. 2012), which could complicate associations between growth habit and genomic regions in the CDBN panel. These four associations were associated with variation in determinacy in this panel; however, these four associations were not sufficient to explain all variation in determinacy, in that 13 genotypes had all alleles that were associated with determinacy, but were indeterminate, and one genotype had all alleles that were associated with indeterminacy, but was determinate (Figure S3C). The association at 6 Mb on Pv01 fell in the coding sequence of the gene model Phvul.001G055600, a RING-CH-type zinc finger protein expressed highly in roots and in stem internodes above the cotyledon at the second trifoliate stage (O’Rourke et al. 2014; Dash et al. 2016). The association at 42.2 Mb was 3.7-kb upstream of the gene VIP5; as noted above, this gene and genomic region were also candidate associations for seed yield (kg ha−1). The association on Pv09 was 5-kb upstream of model Phvul.009G204100, which encodes a signal peptide peptidase A highly expressed in pods associated with stage-2 seeds and in stem internodes above the cotyledon at the second trifoliate stage (O’Rourke et al. 2014; Dash et al. 2016). The association on Pv10 was 1-kb upstream of model Phvul.010G146500, a gene from an uncharacterized protein family highly expressed in roots, pods with seeds at the heart stage, and stem internodes above the cotyledon at the second trifoliate stage (O’Rourke et al. 2014; Dash et al. 2016).
Bean rust (Uromyces appendiculatus) causes leaf and pod pustules, and leads to losses in vigor and seed yield. Higher plant damage caused by rust was indicated by a higher rust score. Rust score had significant associations on 10 chromosomes after FDR (Figure S2E and Table S4). However, the strongest association was located on Pv11 at 50.6 Mb and overlapped a major cluster of disease resistance genes containing the rust resistance genes Ur-3, Ur-6, Ur-7, and Ur-11 (Hurtado-Gonzales et al. 2017). This signal fell just upstream of the gene model Phvul.011G193100, which maps in the interval suggested to contain the resistance gene Ur-3 (Hurtado-Gonzales et al. 2017). The alternate allele was present in the early years of our CDBN data within the Mesoamerican race, but was either absent or rare within the Durango race in the CDBN until 1988, when it was observed in the pinto Sierra and the great northern Starlight. The alternate allele was not widely distributed in the Durango race until the mid-1990s (Figure S3D).
Finally, the presence or absence of curly top virus, a virus characterized by plant stunting and deformation of leaves and fruit, had significant associations on seven chromosomes after FDR; however, the strongest associations were on Pv01, Pv05, Pv07, and Pv11 (Figure S2F and Table S4). The association on Pv01 was 0.5-kb upstream of gene model Phvul.001221100, found to be associated with days to flowering (Kamfwa et al. 2015b), and recently identified as the photoperiod sensitivity locus Ppd, or PHYTOCHROME A3 (Kamfwa et al. 2015b; Weller et al. 2019). The association on Pv05 was within 20 kb of gene model Phvul.005G051400, a VQ motif-containing protein highly expressed in leaf tissue. VQ motif-containing proteins are a class of plant-specific transcriptional regulators that regulate photomorphogenesis, and responses to biotic and abiotic stresses (Jing and Lin 2015). The association on Pv07 was 1-kb upstream of gene model Phvul.007G035300, a pH-response regulator protein. The association on Pv11 was 20-kb downstream of gene model Phvul.011G142800, a terpene synthase gene expressed in young trifoliates, flowers, and young pods (O’Rourke et al. 2014; Dash et al. 2016). Terpenoids are a large class of secondary metabolite, which have roles in plant defense against biotic and abiotic stresses (Singh and Sharma 2015).
The presence of many associations above the FDR threshold supports using MET data for genetic analyses. However, the assortative mating employed purposefully by breeders of entries in the CDBN could potentially lead to a high rate of false positives (Li et al. 2017). Overall, it was unclear whether GWAS using phenotypes derived from sparse MET data sets would yield similar genetic associations as published, balanced field trials. Thus, we compared the top associations discovered here to associations from 11 published GWAS papers on common bean. This allowed us to compare association overlaps for 13 phenotypes, 7 of which had associations above the FDR, and gave 34 top associations from this study to compare to 80 published association regions. In addition to these GWAS associations, the bean rust resistance phenotype overlapped with a candidate rust resistance gene, Ur-3, one of the two genes pyramided in the 1980s in bean breeding to provide comprehensive rust resistance.
Three major associations from this study were within 20 kb of, and had the same candidate gene as, top associations from published, balanced GWAS: days to flowering, on Pv01 at 13.7 Mb; growth habit, on Pv01 at 42.2 Mb; and lodging, on Pv07 at 34.2 Mb (Table 2). Interestingly, when considering all 114 associations, each of these three regions had significant effects for three phenotypes: lodging, growth habit, and days to flowering on Pv01 at 13.7 Mb; growth habit, seed yield, and biomass on Pv01 at 42.2 Mb; and plant height, lodging, and growth habit on Pv07 at 34.2 Mb (Table 2). In this study, the top 10 SNPs for harvest index and days to maturity also had the same candidate gene on Pv03 at 36.8 Mb, the gene model Phvul.003G153100. Phvul.003G153100 is an APETALA2-like ethylene-responsive transcription factor that is highly expressed in root tissue and nodules (O’Rourke et al. 2014; Dash et al. 2016). This region of Pv03 has also been found to have a genetically stable effect on seed yield (Kelly 2018), though our GWAS on seed yield did not detect this region at a 10% FDR.
|FDRa||Trait||Study||Chr||Position in v2.0||Candidate geneb|
|Plant height (cm)||This study||1||6.13c||Phvul.001G054800|
|FDR||Growth habit||This study||1||6.28c||Phvul.001G055600|
|Biomass (kg)||This study||1||6.49c||Phvul.001G057100|
|Lodging score||Resende et al. (2018)||1||13.76c||Phvul.001G087900|
|Growth habit||Resende et al. (2018)||1||13.76c||Phvul.001G087900|
|Days to flowering||Moghaddam et al. (2016)||1||13.76c||Phvul.001G087900|
|FDR||Days to flowering||This study||1||13.45–15.36c||Phvul.001G087900|
|Root rot damage||Oladzad et al. (2019b)||1||23.92|
|Days to flower||Oladzad et al. (2019a)||1||27.68|
|Root rot damage||Oladzad et al. (2019b)||1||33.03|
|Halo blight damage score||This study||1||36.72||Phvul.001G132516|
|Root rot damage||Oladzad et al. (2019b)||1||37.20|
|Plant height (cm)||This study||1||38.74||Phvul.001G143800|
|Root rot damage||Oladzad et al. (2019b)||1||40.20|
|FDR||Growth habit||This study||1||42.17c||Phvul.001G167200|
|FDR||Seed yield||This study||1||42.23c||Phvul.001G167200|
|Growth habit||Moghaddam et al. (2016)||1||42.23c||Phvul.001G167200|
|Biomass (kg)||This study||1||42.27c||Phvul.001G167200|
|Growth habit||Cichy et al. (2015)||1||44.80c||Phvul.001G189200|
|Growth habit||Moghaddam et al. (2016)||1||44.80c||Phvul.001G192200|
|Days to flowering||Nascimento et al. (2018)||1||47.07||Phvul.001G214500|
|Days to flowering||Raggi et al. (2019)||1||48.86|
|Days to flowering||Raggi et al. (2019)||1||49.65|
|Halo blight damage score||This study||2||16.17||Phvul.002G091900|
|FDR||Seed weight||This study||2||30.38||Phvul.002G150600|
|Days to flowering||Oladzad et al. (2019a)||2||38.07|
|Halo blight damage score||Tock et al. (2017)||2||49.08||Phvul.002G326200|
|FDR||Harvest index (%)||This study||3||2.16||Phvul.003G023000|
|FDR||Seed weight||This study||3||4.40||Phvul.003G039900|
|Days to maturity||This study||3||15.75|
|Days to maturity||This study||3||32.04||Phvul.003G128400|
|Harvest index (%)||This study||3||36.82c||Phvul.003G153100|
|Days to maturity||This study||3||36.83c||Phvul.003G153100|
|Seed yield||Kamfwa et al. (2015b)||3||37.60||Phvul.001G136600|
|Days to flowering||Oladzad et al. (2019a)||3||40.27|
|Days to flowering||Oladzad et al. (2019a)||3||41.09|
|Harvest index (%)||Kamfwa et al. (2015b)||3||46.70||Phvul.003G233400|
|Harvest index (%)||Kamfwa et al. (2015b)||3||47.17c||Phvul.003G237900|
|Days to flower||Oladzad et al. (2019a)||3||47.35c|
|Seed yield||Resende et al. (2018)||3||49.28–50.33c||Phvul.003G253700|
|Days to flowering and days to maturity||Oladzad et al. (2019a)||3||51.48|
|Days to flowering and days to maturity||Oladzad et al. (2019a)||3||52.32|
|Days to flowering and days to maturity||Oladzad et al. (2019a)||3||52.6|
|Halo blight damage score||Tock et al. (2017)||4||0.55–1.899c||Phvul.004G007600|
|Days to maturity||Moghaddam et al. (2016)||4||1.94c||Phvul.004G011400|
|FDR||Lodging score||This study||4||2.87||Phvul.004G025600|
|Growth habit indeterminate||Moghaddam et al. (2016)||4||3.20||Phvul.004G027800|
|Halo blight damage score||MacQueen et al. (2020)||4||6.79||Phvul.004G051500|
|Days to flowering||Raggi et al. (2019)||4||16.37|
|Days to flowering||Raggi et al. (2019)||4||36.88|
|Halo blight damage score||Tock et al. (2017)||4||46.20c||Phvul.004G158000|
|Days to flowering and days to maturity||Oladzad et al. (2019a)||4||46.33c|
|Days to flowering and days to maturity||Oladzad et al. (2019a)||4||47.06|
|Halo blight damage score||This study||5||13.25||Phvul.005G074200|
|Halo blight damage score||Tock et al. (2017)||5||39.00||Phvul.005G162500|
|Root rot damage||Oladzad et al. (2019b)||6||0.57|
|Root rot damage||Oladzad et al. (2019b)||6||5.75|
|Root rot damage||Oladzad et al. (2019b)||6||6.89|
|Root rot damage||Oladzad et al. (2019b)||6||8.16||Phvul.006G017211|
|Root rot damage||Oladzad et al. (2019b)||6||12.20|
|Root rot damage||Oladzad et al. (2019b)||6||17.85|
|Plant height (cm)||This study||6||20.89||Phvul.006G098300|
|Growth habit indeterminate||Moghaddam et al. (2016)||6||29.92||Phvul.006G203400|
|Days to flowering||Raggi et al. (2019)||6||31.60|
|Days to maturity||This study||7||1.15||Phvul.007G017000|
|Growth habit indeterminate||Moghaddam et al. (2016)||7||34.12c||Phvul.007G246700|
|Lodging score||Moghaddam et al. (2016)||7||34.20c||Phvul.007G218900|
|FDR||Lodging score||This study||7||33.60–34.51c||Phvul.007G218900|
|Plant height (cm)||Moghaddam et al. (2016)||7||34.20c||Phvul.007G218900|
|Growth habit indeterminate||Moghaddam et al. (2016)||7||34.20c||Phvul.007G218900|
|Biomass (kg)||This study||7||35.74||Phvul.007G233700|
|Seed weight||Moghaddam et al. (2016)||8||1.10c||Phvul.008G013300|
|Root rot damage score||This study||8||1.34c|
|Days to flowering||Raggi et al. (2019)||8||4.93|
|Biomass (kg)||Soltani et al. (2017)||8||6.86||Phvul.008G073000|
|Biomass (kg)||Soltani et al. (2017)||8||7.60||Phvul.008G078200|
|Root rot damage||Oladzad et al. (2019b)||8||15.26|
|Root rot damage||Oladzad et al. (2019b)||8||17.72|
|Days to flowering and days to maturity||Oladzad et al. (2019a)||8||24.95|
|Days to flowering||Raggi et al. (2019)||8||26.40|
|Halo blight damage score||Tock et al. (2017)||8||61.34c||Phvul.008G268700|
|Root rot damage score||This study||8||61.98c||Phvul.008G277352|
|Halo blight damage score||This study||9||5.42||Phvul.009G022400|
|Days to maturity||This study||9||5.85|
|Seed yield||Kamfwa et al. (2015b)||9||10.00||Phvul.009G051600|
|Plant height (cm)||This study||9||27.98||Phvul.009G185100|
|FDR||Growth habit||This study||9||30.93||Phvul.009G204100|
|Biomass (kg)||Soltani et al. (2017)||10||0.60||Phvul.010G003600|
|Seed weight||Moghaddam et al. (2016)||10||2.60||Phvul.010G017600|
|Root rot damage||Oladzad et al. (2019b)||10||16.40|
|Root rot damage||Oladzad et al. (2019b)||10||19.51|
|Root rot damage||Oladzad et al. (2019b)||10||25.44|
|Root rot damage||Oladzad et al. (2019b)||10||29.88|
|Biomass (kg)||Soltani et al. (2017)||10||36.45||Phvul.010G099100|
|Halo blight damage score||This study||10||41.52||Phvul.010G133101|
|Days to flowering||Nascimento et al. (2018)||10||42.50c||Phvul.010G142900|
|FDR||Growth habit||This study||10||42.79c||Phvul.010G146500|
|Biomass (kg)||Soltani et al. (2017)||11||1.59||Phvul.011G020500|
|Days to flowering||Oladzad et al. (2019a)||11||4.02|
|Days to maturity||Moghaddam et al. (2016)||11||4.46||Phvul.011G050300|
|Days to flowering||Oladzad et al. (2019a)||11||10.66|
|days to maturity||MacQueen et al. (2020)||11||14.8|
|Root rot damage||Oladzad et al. (2019b)||11||26.41|
|Days to flowering||Oladzad et al. (2019a)||11||27.27|
|Days to flowering||Oladzad et al. (2019a)||11||36.22|
|Root rot damage||Oladzad et al. (2019b)||11||41.36|
|Days to maturity||Moghaddam et al. (2016)||11||45.09c||Phvul.011G158300|
|Days to flowering||Oladzad et al. (2019a)||11||45.29c|
|Growth habit indeterminate||Moghaddam et al. (2016)||11||46.75||Phvul.011G164800|
|Days to flowering||Oladzad et al. (2019a)||11||47.3|
|Root rot damage||Oladzad et al. (2019b)||11||50.59c|
|Biomass (kg)||Soltani et al. (2017)||11||52.55||Phvul.011G207500|
In comparisons involving only the 11 balanced studies, 9 of 80 associations fell into three 20-kb regions, while 15 of the 80 associations fell into six 200-kb regions. When this study was added, seven additional associations fell into four 20-kb regions, while 12 additional associations fell into 14 overlapping 200-kb regions (Table 2). This study did not identify many new overlaps at the 20-kb level, though it did find associations in all three 20-kb overlapping regions found by comparing the 11 balanced studies alone. However, it did find many new overlaps with previously published studies at the 200-kb level, twice as many as expected given the rate of overlap in the 11 balanced studies (χ2 test P = 0.025). However, as the balanced studies often did not conduct GWAS on similar phenotypes, our “expected” rate of overlap is likely to be biased. Thus, we consider the fact that this study found the same three 20-kb regions that overlap in balanced GWAS comparisons to be stronger evidence than the large number of overlaps at the 200-kb level of the conclusion that this panel can yield similar associations to balanced GWAS of common bean diversity panels.
We observed that numerous CDBN phenotypes had overlapping distributions of significantly associated SNPs. These overlaps could be due to pleiotropy (one genetic locus affecting multiple phenotypes) or due to multiple linked genetic loci affecting multiple phenotypes. To formally compare these overlaps, we used mash on 19 sets of 4000 SNPs with the smallest P-values for phenotypes from the CDBN as well as 4000 SNPs for the earliest year an entry was grown in the CDBN (Figure 3). Mash shares information about effect sizes of SNPs across all phenotypes, while accounting for data-driven covariances in the patterns of effects (Urbut et al. 2019). In contrast to phenotype-by-phenotype analyses, where only 8 phenotypes had associations above the FDR, in mash, all 20 phenotypes had SNPs with P-values below the local false-sign rate, an analog for the FDR. In addition, SNPs typically had local false-sign rates below this threshold for 11–14 phenotypes; thus, there was either extensive pleiotropy or frequent linked effects on multiple phenotypes within entries in the CDBN. SNPs with Bayes factors > ∼102, indicative of decisive evidence favoring that SNP having a significant effect on one or more phenotypes, were distributed very unevenly across the genome, with the vast majority of SNPs clustering within two large regions on Pv01 (Figure 3B and Table S5), from 5.3 to 20.5 Mb, and from 34.1 to 45.4 Mb. These associations encompass two regions of intermediate gene density surrounding the centromere on chromosome 1, but do not include the gene-rich edges of this chromosome (5 and 6 Mb in size, respectively). Interestingly, the two largest Bayes factors across all 20 phenotypes were within these two regions, on Pv01 at positions 15.4 and 42.2 Mb. These associations were two that overlapped with top associations from published, balanced GWAS (Table 2). Outside of chromosome Pv01, the most significant Bayes factor was found for a SNP on Pv07 at 14.5 Mb. This SNP was not within 100 kb of any annotated gene.
The alternate allele for the SNP on Pv01 at 15.4 Mb was associated with significant decreases in biomass, days to flowering, days to maturity, plant height, and seed appearance score. It was also associated with increases in common bacterial blight (CBB) damage score, harvest index, root rot damage score, rust damage score, seed fill duration, white mold damage score, and zinc deficiency damage score (Figure 3D). Here, higher damage scores indicate increased levels of damage. The alternate allele for the SNP on Pv01 at 42.2 Mb was associated with significant decreases in biomass, days to flowering, growth habit (as an increased tendency toward determinacy), harvest index, lodging score, plant height, and seed yield, and increases in root rot damage score (Figure 3E). The allele was also significantly associated with earlier “earliest year in the CDBN,” indicating that this allele has been declining in frequency in entries in the CDBN over time. The alternate allele for the SNP on Pv07 at 14.5 Mb was associated with significant decreases in biomass, days to flowering, plant height, and seed appearance score (Figure 3F). Overall, two groups of phenotypes had consistent patterns of effect sign and effect magnitude for most significant SNPs (Figure 3C). Days to maturity, growth habit, seed yield, days to flowering, biomass, and plant height had a large fraction of SNPs with significant effects with similar effects on these phenotypes; in most pairwise comparisons of these six traits, 40–90% of SNPs had the same sign and similar magnitudes of effect (Figure 3C). The same was true for seed fill duration, white mold damage score, zinc deficiency damage score, harvest index, CBB damage score, and rust damage score; in pairwise comparisons of these six traits, 25–80% of SNPs had the same sign and similar magnitudes of effect (Figure 3C). The phenotypes in the first group corresponded to plant architecture and size, while several phenotypes in the second group were related to disease response. Few other SNPs (< 10%) affected these two clusters of phenotypes to a similar extent with the same sign. Interestingly, groups of highly positively correlated phenotypic BLUPs, or genetic values, did not consistently match groups with large fractions of SNP effects of the same sign and similar magnitude (Figure S4). Overall, 90% of SNPs with Bayes factors > 102 affected ≥ 10 phenotypes (Table S5) and typically affected phenotypes in the two groups in similar ways. However, a few exceptions included Pv03 at 10.64 Mb, which affected only plant height; Pv04 at 17.77 Mb, which affected seed weight and varied with earliest year in the CDBN; Pv07 at 13.94 Mb, which affected biomass; and Pv08 at 33.18 Mb, which affected days to flowering, plant height, and seed appearance.
The genes and genomic regions affecting phenotypic variation in common bean are now being narrowed down with the aid of a recently released high-quality reference genome (Schmutz et al. 2014). The use of previously generated phenotypic data for genetic analysis could circumvent the “phenotypic bottleneck” that has previously constrained our understanding of the genotype–phenotype map in this species. The CDBN offers a vast phenotypic data resource for common bean; however, it was unclear whether the sparse phenotypic data matrix from the CDBN, where only 20–30 entries were tested in each location and year, could be used for GWAS. Our results provide evidence supporting the use of METs such as the CDBN for genetic analysis. First, 8 of the 22 phenotypes created using the CDBN data had associations that fell above the Bonferroni–Hochberg FDR threshold, and 5 of these phenotypes had multiple independent peaks that fell above this threshold. Given our FDR of 10%, there were ≥ 30 distinct, significant associations with these CDBN-derived BLUPs for phenotypes, and these associations tended to be found in phenotypes with higher heritabilities. However, it is still surprising that only 8 of the 22 phenotypes had significant associations by the FDR criterion.
We hypothesized that noise caused by environmental or gene-by-environmental variation in phenotypes across years and locations reduced our ability to find significant associations in a condition-by-condition analysis. Supporting this hypothesis, we found that phenotypes with more data points in the CDBN were more likely to have associations above the FDR. We used mash to increase our power to detect significant effects for 20 of these phenotypes, and used an analog of the FDR, the local false-sign rate, to determine whether an effect was significant. Mash has been demonstrated to have higher power to detect significant shared effects in correlated phenotypes by shrinking effect estimates toward covariance matrices informed by patterns in the data (Urbut et al. 2019). Without mash, we found significant associations for 8 phenotypes; using mash, we found significant associations for all 20 phenotypes. Thus, phenotypes derived from CDBN MET data are suitable for analysis using GWAS, and the additional phenotypic data available in this MET can be analyzed in mash to boost the power to detect significant genetic effects for SNPs that affect multiple phenotypes.
Second, associations found in our GWAS colocalized with results of previous GWAS using balanced phenotypic data sets. Three associations from this study overlapped top associations from published, balanced GWAS: Pv01 at 13.7 Mb, Pv01 at 42.2 Mb, and Pv07 at 34.2 Mb (Table 2). The association at 13.7 Mb fell near the candidate gene KNU, a gene that is activated in, and later promotes, the transition to determinate floral meristem development. This peak falls within an association for days to flowering observed previously (Moghaddam et al. 2016). The association at 42.2 Mb fell near the candidate gene VIP5, an important regulator of flowering time in A. thaliana and other species (Huang et al. 2012). Other mapping studies have also colocated VIP5 with QTL for flowering time (Zhou et al. 2014). The association at 34.2 Mb on Pv07 also overlapped the strongest association for the earliest year each entry was grown in the CDBN, a proxy for the age of the CDBN entry. This association fell near the candidate gene SAR1, which increases plant height and internode distance in A. thaliana (Cernac et al. 1997; Parry and Estelle 2006). A genomic region affecting determinacy in common bean has been confirmed on Pv07; a genomic region affecting determinacy, part of growth habit, would also affect the traits that we find to be affected by this region: biomass, days to flowering, plant height, and seed appearance score. The alternate allele for the signal on Pv07 occurred in newer CDBN entries.
Third, our results are consistent with the recent history of breeding efforts in common bean, and provide a map of the genomic regions that have been associated with improvement in the species. Using mash, we found two major genomic regions on Pv01 associated with many CDBN phenotypes (Figure 3B). Though SNPs for mash were chosen to have low LD (r2 < 0.2), significant SNPs nonetheless highlight two, large linked regions on chromosome 1 that have effects on many phenotypes. Regions of this size (15 and 11 Mb) with evidence for significant effects almost certainly require the presence of multiple linked genes with effects on multiple phenotypes. Linked effects do not preclude pleiotropic effects occurring at single loci within this large region, though molecular confirmation of pleiotropy is still necessary. We suggest that the two major genomic regions on Pv01 associated with many CDBN phenotypes were major targets of selection by breeders for entries that matched an “ideotype” for common bean. The original ideotype had a long hypocotyl, many nodes carrying long pods and without side branches, small leaves, and determinate growth (Adams 1982; Kelly 2001). The primary plant architecture change introduced into genotypes tested in the CDBN over the past 30 years was the adoption of upright indeterminate architecture (type II), which replaced upright determinate (type I) architecture in the Mesoamerican race and was introduced into prostrate indeterminate (type III) germplasm (Kelly 2001; Soltani et al. 2016). Generally, entries with type II architecture yielded more than determinate (type I) entries, due to the increased pod set associated with indeterminate growth (Kelly 2001). Entries with type II architecture could yield more than type III entries under grower-preferred direct harvest (Eckert et al. 2011). An association for growth habit on Pv01 at 42.2 Mb fell near the gene VIP5; this SNP and gene were also candidate associations for seed yield in this study and days to flowering in Moghaddam et al. (2016). The Pv01, Pv09, and Pv10 associations for growth habit alter determinacy and segregate in different genotypes, consistent with the known multiple origins of determinacy segregating in this species (Figure S3C). However, these associations were not sufficient to explain all variation in determinacy present in this panel, perhaps due to the relative rarity of some variants controlling determinacy within the CDBN panel.
Bean breeders in North America generally avoided modifying days to flowering over the years of the CDBN, to protect matching of phenology to specific production environments. However, when type II architecture was introduced from Mesoamerica race into the Durango/Jalisco race, the first entries with this architecture showed delayed flowering, which prevented pod set on lower nodes (Vandemark et al. 2014). Our strongest association for days to flowering was near the candidate gene KNU. This gene is a candidate for the gene Higher response (Hr) (Gu et al. 1998), which affects flowering time. A basic local alignment search tool analysis of random amplified polymorphic DNA primers from previous mapping analyses constrains the location of Hr between 1.4 and 21 Mb on Pv01 (Gu et al. 1998). Hr is thus a plausible candidate for the peak at 13 Mb. Hr is known to be in LD with the common bean gene terminal flower 1 (PvTFL1 or fin) on Pv01, a major determinacy gene in common bean, (Repinski et al. 2012). Thus, this gene could plausibly have been introduced during the introduction of type II architecture.
The primary disease resistance phenotype introduced into entries in the CDBN over the past 30 years was bean rust resistance. Bean rust (U. appendiculatus) was a major disease in North America in the 20th century (Zaumeyer 1947). Although the first rust-resistant varieties were released in the 1940s (Zaumeyer 1947), rust was primarily controlled by chemicals prior to the concerted introduction of rust resistance genes in the mid-1980s (Kelly 2001). Our strongest association for rust damage score fell just upstream of the gene model Phvul.011G193100, which maps in the interval suggested to contain the resistance gene Ur-3 (Hurtado-Gonzales et al. 2017). Initially described by Ballantyne (1978), Ur-3 was the first gene aggressively used by US breeders to address bean rust in the mid-1980s (Hurtado-Gonzales et al. 2017). Combining Ur-3 and Ur-11 provides resistance against all known rust races (Pastor-Corrales et al. 2003), and the two genes formed the basis of breeding efforts to pyramid major bean rust resistance genes that led to the release of pinto, great northern, and black bean germplasm currently used in breeding programs. The alternate allele was present in the early years of the CDBN data in the Mesoamerican race, but was either absent or rare in the Durango/Jalisco race in the CDBN until 1988, when it was observed in the pinto Sierra and the great northern variety Starlight. The alternate allele was not widely distributed in the Durango/Jalisco race until the mid-1990s (Figure S3D). These results agree with the known timing of breeding for rust resistance.
Finally, this work allowed us to characterize the patterns of sharing of genetic effects on phenotypes in the CDBN, using an exciting new method to determine if this data set could uncover novel results. Selection for the common bean ideotype is known to have led to pleiotropic effects on other traits, such as seed yield, biomass, and plant height (Soltani et al. 2016). Previous work has indicated that genes responding to photoperiod have a major influence on many traits, including biomass, harvest index, days to maturity, and plant architecture traits such as the number of branches and nodes (Wallace et al. 1993; Gu et al. 1994). Our associations also revealed substantial overlaps in the genomic regions affecting phenotypic variation, suggesting the presence of substantial pleiotropy or linked genes of major effect. The genomic region on Pv01 from 34 to 48 Mb has also been identified in previous QTL mapping studies as one that affects many traits, including seed yield, days to flowering, days to maturity, seed fill duration, seed weight, biomass, and pod wall ratio (Trapp et al. 2015, 2016). Our mash analysis reveals two major groups of phenotypes with commonly shared SNP effects, one corresponding to plant architecture and size, and the other related to disease response. Very few SNPs had similar effects on both groups of traits (Figure 3C). This indicates pleiotropy or correlated effects within each group of phenotypes, and unlinked effects or antagonistic pleiotropy between these groups of phenotypes. In addition, the two groups of phenotypes that had similar genetic effects at the SNP level did not substantially overlap groups of phenotypes that had highly correlated genetic values by BLUP estimation (Figure S4). Though many genomic regions affect multiple phenotypes in the CDBN, the large shared effects detected by mash do not always combine additively into the overall patterns of genetic correlation present in this data set. However, two sets of phenotypes did have shared SNP effects and similar patterns of phenotypic correlations: lodging, seed yield, and growth, and biomass, plant height, days to flowering, and days to maturity. We suggest that these seven phenotypes were the most important when breeders selected for preferred common bean ideotypes. In contrast, many of the remaining phenotypes were related to disease damage; these phenotypes might be more affected by epistatic interactions between genomic regions or by tradeoffs across environments.
Overall, METs such as the CDBN offer a remarkable opportunity to identify candidate genes underlying phenotypic variation and phenotypic plasticity, and to identify how artificial selection has affected crop phenotypes through time. We note that the genomic regions found with this approach are likely to have consistent, stable phenotypic effects across a large range of environments. These genomic regions are thus likely to be generally useful to bean breeding. Detailed mapping and cloning of the causative genes in these regions will provide insight into molecular mechanisms that control these critical phenotypes important for high productivity of common bean. In the future, we also believe that it would be of great value to crop breeding and genetics to archive DNA from all material used in breeding programs and MET trials.
Many crops, both in the U.S. and worldwide, have public trials that could be mined in a manner similar to our approach. This work will require collaborative efforts between crop breeders and bioinformaticians to digitize, clean, and analyze phenotypic data from METs, and to obtain genetic material from successful and unsuccessful trial entries. Phenotypic and genetic data can be combined using genomic selection approaches, or by GWAS using models that adjust BLUPs for effects of kinship, trial location, and trial year (Rife et al. 2018; Sukumaran et al. 2018). If effect estimates for genetic markers can be obtained and some effects are strong, the patterns of significant effects across markers and phenotypes can be determined using a metanalysis approach such as mash (Urbut et al. 2019). A broader effort to collectively mine such extensive phenotypic data could identify conserved genetic factors important for improved productivity for many crops in major production regions. We present the CDBNgenomics R package and analysis as a resource for other researchers to use the CDBN data set, and to apply these techniques to other MET data sets.
The authors thank An Hang, John Kolar, and Shree Singh for substantial contributions over the history of the CDBN; Melody Yeh for assistance in digitizing and harmonizing the CDBN reports; and thank the Joint Genome Institute for prepublication access and use of the P. vulgaris V2.1 genome and annotation. The U.S. Department of Agriculture is an equal opportunity provider and employer. This research was partially funded by support from the National Science Foundation, Plant Genomes Research Program (grant IOS-1612262 to A.H.M.).
Author contributions: A.H.M., J.M.O., J.W.W., P.E.M., and T.E.J. conceived of the analysis. P.E.M., J.M.O., P.N.M., and J.M. provided seed for varieties in the CDBN. A.H.M., R.L., P.E.M., and J.S. sequenced the CDBN. R.L. and P.E.M. provided sequence data from the ADP and MDP. J.W.W. compiled the phenotypic data from the annual CDBN reports. AHM conducted the analyses with input from J.W.W., P.E.M., and T.E.J. A.H.M. wrote the manuscript with contributions from all authors.