Genetic Associations in Four Decades of Multienvironment Trials Reveal Agronomic Trait Evolution in Common Bean.

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 ), 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.

Figure 1

The CDBN is an extensive multienvironment trial with hundreds of germplasm entries for common bean (P. vulgaris) grown at 75 locations over 35 years. The relationships between germplasm entries (first row and column), locations (second row and column), and years (third row and column) used in the CDBN are shown in the grid. Purple represents genotyped germplasm entries from the Durango race of the Mesoamerican gene pool, blue represents genotyped entries in the Mesoamerican race of the Mesoamerican gene pool, green represents genotyped entries in the Nueva Granada race of the Andean gene pool, and yellow represents entries that were not genotyped. Beans of the two most commonly genotyped market classes from each race are pictured to the right of the color key, and the number of entries of each race that were genotyped is displayed to the left of the color key. On the maps of CDBN locations, pie chart size is scaled relative to the total number of entries grown at that location, and circle saturation represents locations with a greater number of years of data available. CDBN, Cooperative Dry Bean Nursery.

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 ; Rife ; Sukumaran ). 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).

Figure 2

Phenotypic data available in the CDBN, and genetic variation within the two of these phenotypes with the most data. (A) Details of 21 phenotypes present in the CDBN; BLUPs from these phenotypes were used for genome-wide association. Black circles left of the y-axis indicate phenotypes with one or more genetic associations that had P-values above the Benjamini–Hochberg false discovery rate. (B) Manhattan plot of BLUPs for seed yield (kg/ha) from the CDBN data. The dashed lines are the cutoff values for peak significance. SNPs above the Benjamini–Hochberg false discovery rate are above the black line, and those in the 0.001 percentile are above the gray line. (C) Quantile–quantile plot of the goodness of fit of the model on seed yield. (D) Manhattan plot of BLUPs for seed weight (mg) from the CDBN data, as in (B). (E) Quantile–quantile plot of the goodness of fit of the model on seed weight. Distribution of the BLUPs for (F) seed yield and (G) seed weight. Black represents the presence of the reference allele for the top SNP labeled in panels (B) and (D), respectively, and gray represents the presence of the alternate allele. BCMV, bean common mosaic virus; BLUP, best linear unbiased predictor; CBB, common bacterial blight; CDBN, Cooperative Dry Bean Nursery; CTV, curly top virus.

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 ) (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 ; Kirby ; Mackay ; Lin ; McCouch ; MacArthur ; Xiao ; Togninalli ), and have recently been implemented in common bean (Cichy ; Kamfwa ,b; Moghaddam ; Soltani , 2018; Tock ; Nascimento ; Oladzad ,b; Raggi ). 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 ; Visscher ). Human GWAS often look for common variants for common diseases and correct phenotypes for effects of age, sex, and location (Schork ; Mefford and Witte 2012; Zaitlen ). 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 ). 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 ). 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 ). 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 ) 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.

Materials and Methods

Background principles: processing, digitization, and genetic analysis of phenotypic data

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 ).

Phenotypic data processing

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.

Germplasm: CDBN breeding strategies, diversity panel, and SNP data set

CDBN selection framework:

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 ).

Germplasm recovery and sequencing:

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 ), the Durango Diversity Panel (DDP) (Soltani ), and the Andean Diversity Panel (ADP) (Cichy ). 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 ; Moghaddam ) 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 ). 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 . 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 ), 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 ). 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.

GWAS:

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 ) implemented in GAPIT with the optimum level of compression (Lipka ). 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 ). 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 ).

Comparison to published genome-wide associations in common bean:

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 ; Kamfwa ,b; Moghaddam ; Soltani , 2018; Tock ; Nascimento ; Oladzad ,b; Raggi ). 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 ; Soltani ; Oladzad ,b). Both Andean and Middle-American gene pools have been observed to have different SNPs underlying domestication traits (Schmutz ). 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.

Analysis of pleiotropy or linked effects on multiple phenotypes:

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 ). Mash has been demonstrated to increase power to detect effects in analyses of human data (Urbut ), 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).

Data availability

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.

Results

Phenotypic correlations in the CDBN

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.

Table 1

BLUP statistics, with phenotypes ordered as in Figure 2

Phenotype	Vaa	Veb	h2c	BLUP range	Units
Seed yield	53,173	222,409	0.193	1,727	kg ha−1
Seed weight	2,340	1,076	0.685	443	mg
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
Harvest index	13.5	31.5	0.299	20.1	%
Plant height	13.8	42.5	0.245	18.1	cm
Biomass	2.46E+05	6.72E+05	0.268	2,990	kg
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

BCMV, bean common mosaic virus; BLUP, best linear unbiased predictor; CBB, common bacterial blight; CTV, curly top virus; REML, restricted maximum likelihood.

Va is the REML estimate of the additive genetic variance from rrBLUP.

Ve is the REML estimate of the error variance.

h2 is narrow-sense heritability, defined as Va / (Va + Ve).

Eight CDBN phenotypes have genetic associations above the FDR

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 ; Lasky ), 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 ). 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 ; Dash ). 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 ).

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 ; Dash ).

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 ). KNU is activated in the transition to determinate floral meristem development and functions in a feedback loop that promotes determinate development (Payne ; Sun ).

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 ). The signal on Pv05 fell within gene model Phvul.005G005400, a uridine diphosphate glycosyltransferase superfamily protein (Dash ). 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 ; Parry ).

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 ; Dash ).

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 ), 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 ; Dash ). 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 ; Dash ). 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 ; Dash ).

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 ). 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 ). 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 ), and recently identified as the photoperiod sensitivity locus Ppd, or PHYTOCHROME A3 (Kamfwa ; Weller ). 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 ; Dash ). Terpenoids are a large class of secondary metabolite, which have roles in plant defense against biotic and abiotic stresses (Singh and Sharma 2015).

Three CDBN genetic associations overlap genetic associations from balanced genetic field trials

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 ). 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 ; Dash ). 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.

Table 2

Major associations in GWAS from phenotypes from the CDBN and from previously published GWAS

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
FDR	Rust	This study	11	50.67c	Phvul.011G193100
	Rust	Phil McClean	11	50.67c	Phvul.011G193100
	Biomass (kg)	Soltani et al. (2017)	11	52.55	Phvul.011G207500

Chr, chromosome; FDR, false discovery rate; GWAS, genome-wide association study; CDBN, common dry bean nursery.

FDR indicates associations from this paper, which were above the Benjamini–Hochberg FDR correction.

Positions (in Mb) indicate associations in more than one published GWAS that are within 200 kb of one another.

Associations in more than one published GWAS, which fall within 20 kb of one another and have the same candidate gene.

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.

Extensive pleiotropy or linked effects within CDBN genetic associations

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 ). 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.

Figure 3

Patterns of phenotypic effects of genetic associations for 22 phenotypes from the CDBN, determined using multivariate adaptive shrinkage. (A) SNPs with significant effects on 1 or more of the 22 phenotypes in the CDBN. (B) Manhattan plot of the Bayes factor (log10) comparing the model likelihood that the SNP has significant effects to the likelihood that it has no significant effects. Bayes factors of > 102 are considered decisive evidence in favor of the alternate model. Point color represents the number of phenotypes for which the SNP has a local false-sign rate < 0.05. Squares represent even chromosomes, while circles represent odd chromosomes. The top associations for three regions of the genome are highlighted. (C) Correlation in the sign and magnitude of significant effects in all pairwise comparisons of the 22 CDBN phenotypes. Circle size and color indicate the fraction of all significant SNPs that have the same effect sign and similar effect magnitude. (D–F) Effect estimates and standard errors for 22 phenotypes for the top associations from three regions of the genome, (D) P. vulgaris chromosome 1 (Pv01) at 15.4 Mb, (E) Pv01 at 42.2 Mb, and (F) Pv07 at 14.5 Mb. Genomic locations are based on the P. vulgaris v2.1 genome annotation. Standard error bars are colored by the six groups present in (C). CBB, common bacterial blight; CDBN, Cooperative Dry Bean Nursery; CTV, curly top virus.

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.

Discussion

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 ). 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 ). 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 ). 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 ). Other mapping studies have also colocated VIP5 with QTL for flowering time (Zhou ). 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 ; 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 ). 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 ). 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 . 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 ). 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 ), 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 ). 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 ). 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 ). 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 ). Combining Ur-3 and Ur-11 provides resistance against all known rust races (Pastor-Corrales ), 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 ). 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 ; Gu ). 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 , 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 ; Sukumaran ). 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 ). 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.