QTL analyses for tolerance to abiotic stresses in a common bean (Phaseolus vulgaris L.) population

Common bean productivity is reduced by several abiotic stress factors like drought and low soil fertility, leading to yield losses particularly in low input smallholder farming systems in the tropics. To understand the genetics of stress tolerance, and to improve adaptation of common bean to adverse environments, the BAT 881 x G21212 population of 95 recombinant inbred lines (RILs) was evaluated under different abiotic stress conditions in 15 trials across four locations in Colombia, representing two higher altitude (Darién, Popayán) and two lower altitude (Palmira, Quilichao) locations. Stress vs non-stress treatments showed that yields were reduced in drought trials in Palmira by 13 and 31%, respectively, and observed yield reductions in low phosphorus stress were 39% in Quilichao, 16% in Popayán, and 71% in Darién, respectively. Yield components and biomass traits were also reduced. Traits linked to dry matter redistribution from stems, leaves and pods to seed, such as pod harvest index and total non-structural carbohydrates, were found to be important factors contributing to yield in all conditions. In contrast, early maturity was correlated with improved yield only in lower altitude locations, whereas in higher altitudes delayed maturity promoted yield. Superior RILs that combine stress tolerance and high cross-location productivity were identified. Lines that showed good yield under strong stress conditions also performed well under non-stress conditions, indicating that breeder’s selection can be applied for both conditions at the same time. Quantitative trait loci (QTL) analyses revealed a stable yield QTL on chromosome Pv04, detected individually in all locations, several stress treatments and in best linear unbiased predictions (BLUPs) across all trials. Furthermore, two QTL hotspots for maturity traits were identified on Pv01 and Pv08, which are the most stable QTL. The constitutive yield QTL could serve as a good candidate for marker development and could be used in marker assisted selection. Increased understanding of the physiology of abiotic stress tolerance, combined with the availability of superior germplasm and molecular tools, will aid breeding efforts for further improvement of these plant traits.

Increased understanding of the physiology of abiotic stress tolerance, combined with the availability of superior germplasm and molecular tools, will aid breeding efforts for further improvement of these plant traits.
Several years of field evaluation for drought resistance and acid soil tolerance identified the landrace G21212 to be adapted to these stresses. In contrast the breeding line BAT 881 is sensitive to both stresses [24][25][26]. Both lines belong to race Mesoamerica. The authors concluded that productivity under low P is heritable but is based on multiple mechanisms, hence, there may be opportunities to recombine genes and mechanisms and obtain higher levels of stress tolerance.
In recent years new genotyping technologies have been developed, utilizing several thousands of Single Nucleotide Polymorphisms (SNPs) that have become available through next generation sequencing (NGS). SNP based markers have become the preferred choice for genetic studies, due to their abundance and uniform distribution throughout genomes [27]. They have been applied to study diversity, genome mapping and QTLs detection for biotic and abiotic stresses in common bean [28][29][30]. However, adoption into breeding processes has been more limited to biotic stress resistance. Kompetitive Allele Specific PCR (KASP) genotyping technology offers high precision and robustness in assays at low cost [31]. Approximatelỹ 1500 KASP markers have been established for the common bean community by the Generation Challenge Program, to support genetic studies and breeding programs.
The main objectives of the current study were to increase the understanding of adaptation of common bean to different stress environments, and to perform multilocational evaluations of yield stability and stress tolerance. The Mesoamerican RIL population BAT 881 x G21212 (B x G) was evaluated for yield and yield components, plant vigor, dry matter redistribution and physiological traits under different conditions of abiotic stress (drought conditions, low soil fertility and medium-high aluminum levels). The objective of this study was to identify QTLs associated with agronomic traits and to provide markers for use in marker assisted selection (MAS) employing a saturated genetic B x G map.

Plant materials
A total of 100 common bean (P. vulgaris) genotypes were used in this study. These included 95 recombinant inbred lines (RILs) [32], the two parents (BAT 881 and G21212) used as controls and three other lines used in the bean breeding program of the International Center for Tropical Agriculture (CIAT) as additional controls in each trial. Darién, Popayán 1999 and CIAT-Palmira trials used G 3513, G 4017 (tolerant to low fertility) and DOR 364 (sensitive to low fertility) as controls; Quilichao trials used DOR 364, BAT 477 and VAX 1 (tolerant to aluminum and common bacterial blight). RILs were developed by single seed descent (SSD) from the F2 to the F5 generation, followed by seed advance to the F7 generation. BAT 881 is an advanced breeding line of CIAT, with indeterminate bush growth habit (IIa) and small light brown seed that is resistant to bean common mosaic virus (BCMV) and intermediately resistant to the leafhopper (Empoasca kraemeri) and angular leaf spot [32]. G21212 is a small black seeded Colombian landrace, growth habit IIIb, tolerant to low soil fertility, drought and with excellent yield, but is susceptible to BCMV [5,24,32,33]. Both genotypes showed intermediate levels of resistance to the melon thrips [32] and belong to the same race Mesoamerica of the Mesoamerican genepool.

Description of experimental locations
Phenotypic data were obtained in four different locations in Colombia under field conditions with 10 x 10 balanced lattice design with three replicates, except three trials at Quilichao (Q02a_MAl-LP, Q03a_HAL_MP, Q03b_HAL_MP) that had two replicates. A trial refers to an evaluation of all 100 genotypes in a specific location, season, and management, a replicate represents 100 experimental plots harboring all genotypes and a block refers to a sub section of 10 plots within a replicate.
Description of locations: Darién (Dar, Inceptisol, 3˚55'N latitude, 76˚28'W longitude, 1457 masl, low P stress), Popayán (Pop, Inceptisol, 2˚25'N latitude, 76˚40'W longitude, 1730 masl, low P stress), Palmira (Pal, Mollisol, 3˚30'N latitude, 76˚21'W longitude, 965 masl, drought stress) and Santander de Quilichao (Qui, Oxisol, 3˚06'N latitude, 76˚31'W longitude, 990 masl, high Al stress). Specific soil conditions are shown in S1 Table. Two trials were established in Darién: 1997 (September-December) with high P supply (20 kg ha -1 supplied, Dar97_HP) and in 1998 (September-December) with low and high P supply (9 and 90 kg ha -1 , Dar98_LP and Dar98_HP respectively). Experimental units (4 rows of 4 m length with row-to-row distance of 0.6 m) were sown manually with 10 plants per m row length. Two trials were sown in Popayán: 1999 (April-September) with low and high P (0 and 40 kg ha -1 , Pop99_LP and Pop99_HP,) and 2005 (June to August) with high P supply (39 kg ha -1 , Pop05_HP). Experimental units (4 rows of 4.9 m length with row-to-row distance of 0.6 m) were tractor sown with 14 seeds per m. Two trials were established in Palmira: 2000 and 2002 in the same semester (June to September), with two levels of water supply (drought and irrigated, Pal00_D, Pal00_I, Pal02_D, and Pal02_I). Both environments were irrigated to ensure good plant establishment at early growth as needed until early pod filling. Thereafter, irrigation was discontinued in the drought treatment while irrigated or non-stress (NS) plots continued to receive supplemental irrigation of 20 mm of water twice a week for a total of additional 188 and 211 mm of water in 2000 and 2002, respectively (S1 Fig). Experimental units (4 rows of 5 m length with row-to-row distance of 0.6 m) were tractor sown with 14 seeds per m row length. Finally, four trials were sown in Santander de Quilichao: 2002 (April to June) in moderate Al saturation (40 and 64% in top 0-10 and 10-20 cm, respectively), low available P content (3.3 and 1.9 μg g -1 ) and high organic matter content (10.7 and 9.1%), supplementing with P 15 kg ha -1 (Qui02a_MAl_LP). A second trial in 2002 (October to January 2004), in moderate Al saturation (40 and 45%), moderate available P content (17.1 and 6.9 μg g -1 ) and moderate organic matter content (5.6 and 5.2%) soil with low P supply (10 kg ha -1 , Qui02b_MAl_MP). In 2003 (March to June), with high Al saturations (61 and 62%), moderate available P content (23.5 and 13.1 μg g -1 ) and moderate organic matter content (5.9 and 5.4%) with low P supply (10 kg ha -1 , Qui03_HAl_MP). Lastly, in 2003 (October to January 2004) two trials were sown with high Al (61-61% and 64-67%), high and moderate available soil P content (20.3-17.5 μg g -1 and 13.2-5.5 μg g -1 ), moderate organic matter content (5.6-5.4% and 6.1-5.4%), with low P supply (10 kg ha -1 , Qui03_HAl_HP and Qui03_HAl_MP, respectively). Experimental units (2 rows of 3.72 m length in Qui02a_MAl_LP, 4 rows of 4.9 m in Qui02b_MAl_MP and 3 rows of 4.9 m length in the other Quilichao trials) were tractor sown with 14 seeds per m row length and rows separated with 0.6 m.
In the supplemental material S1 Fig the conditions of precipitation and evaporation, minimum and maximum temperature, and total rainfall where available are detailed. In all locations recommended agronomic procedures were followed.

Phenotypic traits
Five groups of plant traits were evaluated in this study: yield components, plant vigor, dry matter redistribution, phenological traits, and mineral nutrients. Not all variables were measured in every trial.
Yield components: Yield (Yd, kg ha -1 ) and 100 seed weight (100SdW, g 100 seeds -1 ). At the time of harvest, 0.3 m 2 per plot were harvested to measure seed number per area (SNA, seeds m -2 ), pod number per area (PNA, pods m -2 ). Plant vigor: Canopy biomass (CB, kg ha -1 ), leaf area index (LAI, m 2 m -2 ), and leaf biomass (LB), pod biomass (PBMP) and stem biomass (SBMP, kg ha -1 ) were evaluated at mid-pod filling. Plants from 0.3 m 2 per plot were harvested, leaves were separated from plants to measure leaf area (LI 3000; LI-CORm Inc., Lincoln, NE). Separated plant parts were oven-dried at 60˚C for two days until establishing constant weight, to determine dry weight. Dry matter redistribution: Indices such as harvest index (HI), pod partitioning index (PPI), pod harvest index (PHI) were calculated. Seed and shoot total nonstructural carbohydrates content (TNC mg g -1 ) were measured to explore genotypic differences in dry matter redistribution and grain filling [14,16,18,33,34]. These dry matter redistribution traits were determined as HI (%): (seed biomass as dry weight at harvest) / (total canopy biomass as dry weight at mid-pod filling) x 100; PHI (%): (seed biomass as dry weight at harvest) / (pod biomass as dry weight at harvest) x 100; PPI (%): (pod biomass as dry weight at harvest) / (total canopy biomass as dry weight at mid pod filling) x 100; seed and shoot TNC were determined by comparison with glucose standards according to the method described by [35]. Because of leaf fall in common bean during seed filling, dry weight of total canopy biomass was determined at mid-pod filling [18]. Use of dry weight of canopy biomass and seed biomass at harvest time overestimates HI and PPI values in common bean. Phenological traits: Days to flowering (DF, days); days to harvest (DH, days). DF is defined as the number of days after planting until 50% of the plants have at least one open flower and DH is the number of days after planting until harvest. Mineral nutrients: One plant per plot was harvested at mid pod filling, washed with deionized water, oven-dried at 60˚C for two days and ground in a 1 mm sieve. For seed evaluations seeds from 0.3 m 2 plot area were oven-dried at 60˚C for two days to be ground in a 1 mm sieve. Nutrient (N, P) composition of plant parts was determined as described by [36]. Shoot nutrient uptake was calculated as the product of shoot dry matter yield and nutrient concentration. Nutrient use efficiency (NUE_Sh, PUE_Sh) were estimated as grams of shoot biomass produced per gram of shoot nutrient uptake [37] and nutrient use efficiency for seed yield (NUE_Sd, PUE_Sd) were estimated as grams of seed biomass produced at harvest per gram of shoot nutrient uptake at mid pod filling. Additional information on description of phenotypic traits can be found in "Trait Dictionaries for Fieldbook Development" at http://mbp.generationcp.org and http://www.cropontologycurationtool.org/.

Statistical analysis
Variance and covariance analyses were carried out, and the results of the covariance are presented only for those traits that showed significant co-variation between genotypes and treatments. Pearson's correlations between all traits were determined for all trials using data of individual replications. Data shown in this work are adjusted means for each genotype and the environment, that were obtained using SAS (v 9.3) PROC MIXED and PROC CORR [38], considering the effects of replication and blocks within replications as random and genotypes as fixed. Also least significant difference (LSD) was calculated with the same software to compare trait means. Variability and frequency distributions were represented as violin plots for each environment evaluated, generated by the statistics software R [39] with the package ggplot2 [40]. To estimate the overall adjusted means and adjusted means for each of the four locations, Best Linear Unbiased Predictions (BLUP) were used, according to Resende [41] for the four major traits (yield, 100SdW, DF, DH). We used the package lme4 [42] for analysis in the statistical software R. The following mixed linear model was defined: Y = B0 + B1 trl + B2 rep + B3 blk + B4 + err where Y is the vector of observed values, B0 is intersect, B1 is the effect of each trial, B2 the effect of each repetition, B3 the effect of each block within trials and replicates, B4 the effect of each genotype or line and err the residual error. Treatments and replicates were considered fixed effects while blocks and genotypes were considered random. The new BLUP variables were named after their respective traits (yield, 100SdW, DF, DH) across each location, e.g. Yd_D, Yd_P, Yd_C and Yd_Q, and for overall adjusted means Yd_all. We also calculated the differential phenotypic response under stress for yield, expressed as percent yield loss under irrigation vs drought (%YdL_D) and high vs low phosphorus (%YdL_LP). They were calculated: %YdL_LP = ((BLUP_higher P -BLUP_lower P) / BLUP higher P) Ã 100 (combining the trials Dar98, Pop99, Qui03b); and %YdL_D = ((BLUP irrigated -BLUP_drought) / BLUP irrigated) Ã 100 (Pal00 and Pal02 trials). For these BLUPs only data sets were used where stress vs non-stress conditions were evaluated in the same season. The same formula was applied for the calculation of percent loss in 100SdW, DF and DH.

DNA, markers and QTL analysis
DNA was extracted from RILs and parents according to the method described before [43]. DNA quality was evaluated on 0.8% agarose gels followed by quantification on a Hoefer DyNA fluorometer (DNA Quant™ 200). The genetic map of the B x G population was constructed from 339 markers, from which 115 markers correspond to previous analyses [32] and we added 53 AFLP, 2 RAPD, 42 SSRs [44][45][46][47][48]. PCR amplification conditions were described before [44,49,50] and PCRs were conducted in 96-well plates using a PTC-100 thermal cycler (MJ Research, Watertown, USA). We also added 127 SNPs markers genotyped at LGC Genomics (www.lgcgenomics.com, Hoddesdon, UK) with KASP genotyping technology. These markers had been established by the Generation Challenge Program (https://www. generationcp.org/), based on the BARKBean chip [30,51]. The genetic map was constructed using MapDisto v. 1.7.5 (http://mapdisto.free.fr/), using Kosambi mapping function, with a minimum LOD score of 3 and maximum recombination fraction of 0.35; SARF and seriation as criteria for ordering and ripple. All genetic mapping was confirmed with the command for ''best order" in MapDisto. Linkage groups were named based on the chromosome information of Phaseolus vulgaris genome v2.1 https:// phytozome.jgi.doe.gov/pz/portal.html (consulted July 2018). QTL analysis was performed through QTL IciMapping software V4.0 based on the inclusive composite interval mapping of additive (ICIM-Add) [52]. The mapping parameters of each step for ICIM-ADD were set at 1.0 cM. LOD threshold value of 2.96 was applied for all traits, determined through completing 1000 permutations with the Type I error level set at α = 0.05 for each mapping method.

Phenotypic data from four locations with abiotic stress treatments
The BAT 881 x G21212 RIL population was evaluated at four locations, totaling 15 trials with independent location/treatment combinations. Two trials were grown under drought stress, four in low P stress, and another three in high Al stress environments. A total of 31 traits were evaluated, only yield and DH were evaluated in all 15 trials (Fig 1 and Table 1), resulting in 225 data sets, including the averages for each locality for the major traits (Yd, 100SdW, DH, DF) and their overall averages. RILs presented continuous population distributions of traits by treatment (Fig 1 and S2 Fig). Transgressive segregation was also evident in B x G lines with lower or higher values compared to parents.

Drought and low P stress
Drought stress and non-stress treatments were compared in two seasons at Palmira (Pal) location. Drought reduced yield by 13% and 31%, respectively (Fig 1 and Table 1). Lowest precipitation among drought trials during pod filling stage was recorded in 2002 (88 mm vs 147 mm, S1 Fig), leading to stronger yield reduction and lower seed weight (Fig 1). Similarly, yield components PNo and SdNo were diminished in 2002, by 40% by and 41%, respectively, whereas seed size was not affected (S2 Fig). Biomass traits were significantly reduced by drought in the Pal 2002 trials (CB = 42%, LAI = 26% LB = 28%, SBMP = 38%, PBMP = 54%). In contrast, in Pal 2000 where the drought was less intense and more rainfall coincided with the pod filling phase no major reductions in biomass were observed. Dry matter redistribution trait values were elevated under drought, suggesting better allocation of dry matter resources to seed under stress (Table 1, S2 Fig). The lack of this data in the Palmira 2000 trial did not allow us to make a comparison between these traits between the two years. In most trials shoot traits LAI and LB showed positive yield correlations, however, in drought conditions negative yield correlations were observed (S2 Table). SdN was higher under drought conditions, whereas SdP was reduced.
Effects of low P stress were evaluated comparing stress and non-stress trials at three locations. Yield reductions were more severe in Darién (Dar) where low P reduced yields by 71% QTL detection for abiotic stress tolerance in common bean (Phaseolus vulgaris L.) (Fig 1, S2 Fig, Table 1). In Popayán (Pop) yield was lowered by 16%, in a season probably affected by excessive rainfall (S1 Fig). In Pop99_HP exclusively positive yield correlations with DF and DH were found, in contrast, in Pop99_LP all correlations were negative (S2 Table). This condition-effect was not observed in Dar or Quilichao (Qui). In trials at Qui yields under moderate P levels (MP) were reduced by 39% compared to high P (HP) conditions. Effects of reduced P on yield components, biomass and dry matter redistribution traits were only evaluated in Qui. All were found to be diminished in moderate P conditions, including all dry matter redistribution traits which was unexpected. Even though shoot biomass production was reduced, lower values of HI, PHI and PPI indicate that seed formation and filling were even more compromised. SdP was reduced by 9% under low P conditions. Qui 2002 trials under moderate Al stress did not show the same tendencies as the trial in MP conditions out-yielded the HP treatment, however, these were sown in different seasons, where climatic differences likely having a larger effect on traits than soil P levels. In addition, good performance of Qui02_MAl_LP trial may be due to higher organic matter content in the soil (S1 Table).
The differential phenotypic expression under stress was evaluated for the four major traits (Fig 2). Yield was reduced more under LP than under drought across these experiments. Most significant correlations between yield loss under stress (%YdL_D and %YdL_LP) and yield were as expected negative correlations in the respective stress trials (drought or LP, S2 Table). Furthermore, several positive significant yield correlations were found, so lines with high yields are more likely to see yields reduced under stress. Negative correlations were observed with pod and seed number as well as pod biomass, suggesting that these yield component traits aid stress tolerance. Positive correlations with maturity traits indicate that late maturing lines loose more yield under stress. Maturity traits saw more reductions in drought trials, under low P DF values were even extended. Accordingly, many more significant trait correlations were found for %DHL_D compared to %DHL_LP, generally contrasting the DH correlations.
Qui 2002 and Qui 2003 trials were planted on soils that were evaluated to have moderate and high Al levels, respectively. No significant effect attributed to differences in soil Al on yield and other traits could be detected across seasons.
Superior RIL lines that performed well under stress and across trials were identified. Yields under most severe drought and low P stress trials were compared to average yields across all 15 trials (Fig 3A and 3B). High correlations indicate that most stress tolerant lines are among the most productive in all trials, displaying high general yield stability. This means that selection of best lines under stress conditions are also generally well performing under optimal conditions and vice versa. RILs R9, 47, 154, are the best lines to combine drought and cross-location yield stability, and another seven lines are excellent for both. For low P tolerance the correlation was also positive, but less clear as the top 5% of lines for each condition do not overlap. The best lines that combine cross-location yield stability and low P tolerance are R30, 45, 80, and the parental line G21212, which is excellent in low P. Taken together, drought and low P stress treatments had significant effects on productivity, yield components and biomass, however, variability between trials due to climatic and other conditions was also large. RILs 1, 80, and 154 can be useful as parents to develop multiple stress resistant genotypes.
Among the parental lines, G21212 parent showed higher yield and earliness than BAT 881 under LP stress in Pop and Qui, drought and combinations of high Al and moderate P (Fig 1). BAT 881 on the other hand, was superior under favorable conditions underlining the contrasting behavior of the parental lines. The additional controls G 3513 and G 4017 used in the Dar and Pop localities were similar in yield, their performance being generally intermediate between the two parents under low P. G 3513 performs poorly in Pal. The additional checks DOR 364, BAT 477 and VAX 1 selected for contrasting aluminum stress response are mostly inferior to G21212.

Effect of locations on traits
Comparisons of phenotypic data sets showed strong correlations among trials. DH showed significant positive correlations between nearly all data sets (Fig 4A), demonstrating that the data set is useful for comparative analysis and variable response is highly stable. Yield correlations between locations were expectedly lower, all significant correlations were expectedly positive ( Fig 4A). Best correlations were found between trials at Dar and Pop, which are at higher altitudes (1457-1730 masl, respectively) compared to lower altitude sites Pal and Qui (965-990 masl). Fewer significant yield correlations were observed between Pal and Qui locations, probably due to their contrasting soil fertility conditions.
Phenotypic data for major traits varied widely among trials (Fig 1). Highest yields were observed in Dar location under non-stress environments with average yields of 3.7 t/ha (Table 1). Pop was the most variable. Lower yields were observed in Pal and Qui. High altitude location Pop recorded most delayed flowering and maturity with highest DF and DH values. Highest seed sizes were observed at Qui location, a trait not evaluated at Dar.
Site regression analysis based on yield data (Fig 4B) resulted in Pop05_HP being the most distinct trial, Dar98_HP and Dar97_HP appear clustered to each other, different from all other data sets that appear in one large group. This result mainly mirrors average yield levels. Site regression did not generally result in distinct clusters for locations. Furthermore, no clustering of trials was observed that were subjected to the same conditions like drought, irrigation, high and low P. Variation between seasons due to climate or other factors seem to have larger effects on yield and other traits than locations or stress treatments.

Yield, yield components and dry matter redistribution traits
Correlations between 224 trait data sets that were investigated are shown in S2 Table. Adjusted mean yield across all trials (Yd_all) was significantly positively correlated with all individual yield trials, underlining the general correlation among all experiments (Fig 4A). Location means for Dar and Pop highly correlated with trials from these sites, but not significantly with Pal and Qui trials, showing the differentness of the two higher altitude locations. For traits 100SdW and DF all trials and location means were highly correlated to each other. Taken together, this confirms that phenotypic data of all trials is comparable showing no drastic differences among trials.
Yield component traits PdNo and SdNo were, as expected, positively correlated with yield, as well as virtually all of 100SdW data sets (S2 Table). Regarding biomass traits, CB and PBMP show predominantly positive yield correlations. Within all 15 yield data sets, the Pal02_I trial had the highest number of correlations with other traits and also produced the highest number of QTLs. We interpret that as a high quality trial that can be used as a representative data set. Fig 3C shows positive yield correlations with CB and 100SdW in this most representative data set Pal02_I.
An evaluation of phenology traits revealed no connections between DF and yield in five trials in Dar and Pop, but significant negative yield correlations in all five Qui trials (not evaluated in Pal). Commitment to early initiation of reproductive development appears constitutively important in this population at lower altitude Qui location. This trend is equally clear in the DH trait, where all significant DH-Yd correlations were positive in Dar and Pop and negative in Pal (Fig 3D) and Qui. This shows different mechanisms for success, where extended seed fill is useful in cooler climates in high altitude locations, whereas early maturity is preferable in lower altitudes.
Dry matter redistribution and transport plays an important role in yield accumulation. PHI showed good yield correlations in Pal, Qui and Pop locations (not evaluated in Dar), reaching up to 0.71 in trial Pal02_I (Fig 3D). HI and PPI had positive yield correlations of lesser significance. Results underline utility of the PHI trait to select high yielding lines, holding up across seasons and locations and showing no particular preference for stress conditions. TNC is an interesting trait, with positive yield correlations if evaluated in shoots (Fig 3D, Table 2), but hardly in seed (S2 Table). Investigating further the relationship with biomass, reveals positive correlations with CB, but not with biomass of vegetative tissues, leaves and stem, where all significant correlations resulted negative (Table 2). Positive correlations were again observed with PBMP, hence the positive effect on CB (at mid pod fill) is solely due to pod mass. TNC_sh correlations were very negative with phenology traits DF/DH. Taken together, high TNC_sh levels lead to pod and seed formation and are not linked with vegetative biomass.
All significant SdN-Yd correlations are negative within trials where both traits were scored (Fig 3E), SdP mainly follows the same tendency, but has one positive yield correlation in Qui. The negative relationship is likely a dilution effect from carbohydrate formation which has a higher variability due to genetic and environmental factors. Interestingly, also ShN and ShP showed only negative correlations with yield, and similarly yield components and TNC_sh. These traits were only measured in Pal and Qui locations. Furthermore, ShN evaluations were strongly positively correlated with nearly all DF and DH data sets, suggesting a link between late maturity and ShN. Hence higher protein/ammonia content of leaves may cause (or be caused by) increased vegetative development and later maturity, which has the aforementioned negative effect on grain productivity in these locations. Surprisingly this was not found for ShP, the only two significant correlations with DF/DH were negative, suggesting a different mechanism for P.
RILs R8, R53 and R67 combined good SdN with above average yield (Fig 3E), though, a solid negative correlation did not allow to identify lines that excel in both traits. Line R53 is above average on general, stress yield and SdN content and may be used for breeding to combine these traits. Taken together, combining high productivity and N content is not easy but to a certain extent possible. Yield components and dry matter redistribution traits PHI and TNC_sh are generally useful proxies for yield, whereas, effects of maturity and biomass traits are site and condition dependent.

Identification of QTLs
An improved genetic map of the B x G population was generated that is suitable for QTL analysis. Genetic map B x G with 8 linkage groups, using a total of 115 markers covering 611 cM was previously reported [32]. We added AFLP (53), RAPD (2), microsatellite (42) and SNP (127) markers for a total of 339 markers, to improve map density and enable high precision QTL mapping in the B x G population (S3 Table). The total length of the genetic map is 1068 cM, which represents an increase of 457 cM over the previous map, with an average distance between markers of 3.15 cM. A total of 143 QTL were identified on all chromosomes for 25 of the 31 evaluated traits. 115 QTLs located in chromosomes Pv01, Pv04, Pv07 and Pv08 are shown in Fig 5 and Table 3 (complete QTL list in S4 Table).

listed in
A consistent yield QTL was identified at the start of chromosome 4 (Fig 5). QTL Yd4.1 was detected in 6 out of 15 trials in all four locations, including Al and low P stress treatments, and means across Dar, Pop, and all sites, indicating a very stable condition-independent QTL. The G21212 (G) derived allele explained up to 30% of phenotypic variability and boosted yield by an average of 112 kg/ha. Yd4.1 is co-localized with yield component QTLs for PBMP, PNo and 100SdW, positively correlated traits also supported by the G allele. Another yield QTL Yd7.1 was observed in two trials in Dar and Pop next to Dar location and across-locations yield means, explaining up to 23% also powered by the G allele. QTL for differential response under stress %YdL_D4.1 and %YdL_LP7.1 were discovered in these same regions. Also % DHL_LP1.1 and %DHL_D8.1 were found in similar regions as their corresponding trait, however, two QTL (%YdL_LP7.1 and %DHL_D8.1) were supported by the opposing allele. In some cases the positive allele that corresponds to higher values for the trait leads to more loss under stress, in the other cases the positive allele leads to more stress tolerance. %SdWL_LP8.1 is the only QTL that was solely identified for differential response under stress but not for the trait itself. QTL hotspots for maturity were found on chromosome 1 and 8 (Fig 5). DF scored 31 QTL, the most of all traits even though not measured in all trials. On chromosome 1, 18 DF QTLs and 4 DH QTLs were detected, all supported by BAT 881 (B) parent; DF1.2 is strongest, explaining 12-52% of variation, data suggests three separate loci on this chromosome. On chromosome 8, 10 DF QTLs and 14 DH QTLs were detected, explaining up to 17% and 32% of variation, respectively. In this chromosome all positive alleles originate from G, data suggest 5 independent loci. These two QTL hotspots for maturity cover the majority of the most significant genotype-phenotype associations.
For TNC_Sh several yield adjacent QTL were identified. TNC_sh6.1 was observed in two Qui trials, and co-localizes with Yd6.1 for Yd_all. The B allele explained 9-13% of variability

Bean productivity is affected by abiotic stress factors
This project produced a rich source of data, collected over 4 locations, 6 years, and three stress treatments. Nine out of fifteen trials had average grain yields lower than 1000 kg/ha, including most stress trials. This resembles typical farmer field yields in the tropics, hence, the data set is comparable and outcomes transferable to breeding efforts aimed at tropical smallholder farming systems. Drought is the major abiotic stress affecting common bean and other legumes [53]. While approximately 67% of globally cultivated lands are affected by P deficits [54], low P tolerance in legumes has been mainly studied in common bean [53]. A deeper understanding of these traits is required to improve productivity under stress. In the B x G population drought and low P significantly reduced yields within evaluated location-season combinations, reductions of 13-31% at Palmira, 16% at Popayán, 71% at Darién and 39% at Quilichao were observed, in line with many field studies reviewed before in common bean, [18,23]. Drought stress caused yield reduction between 47 to 69% compared to the non-stress conditions [55], similarly, a reduction of 31% in the same conditions was reported [56]. For low P conditions compared to high P supply in same location reported yield reductions range between 50-60%. In other legumes, yield reductions due to stress have also been described in soybean, observing 36% yield differences under drought [57,58]. In drought trials, significant correlations with yield provided support for the important role that adaptation plays in selecting for drought resistance across environments, even though drought trials are inevitably noisier than non-stress trials. Yield correlations with vegetative biomass traits were found to be mostly negative, while correlations with overall biomass remained positive as previously reported [11,13]. Yield correlations with dry matter redistribution indices under drought were strong. Early maturity was strongly linked to increased productivity in trials in Qui and Pal. Next to drought avoidance and tolerance, drought escape of shorter duration varieties has been found to be an important trait in chickpea [59], cowpea [60] and lentil [61]. This indicates that focus on dry matter redistribution towards seed production has to be emphasized in water limiting conditions, as indicated before [62].
Similar to drought stress, low P conditions also reduced yields, yield components and biomass traits, in accordance with previous reports [63,64]. Interestingly dry matter redistribution trait values were reduced in the lower P treatment in Qui, the only low P trial in which these traits were evaluated. This contrasts observations under drought where the efficiency of dry matter distribution to grain was increased. In line with these findings, positive Yd-PHI correlations in low P along with reduced mean PHI values in low P were previously reported [64]. Retaining higher PHI under low P may represent a strategy for improved low P tolerance. QTL for differential trait response under stress were found mostly in the same regions as the corresponding traits. However, the alleles that lead to yield loss under stress are not always the same as those that lead to higher trait values indicating different mechanisms: Some alleles support high trait values under non-stress conditions and are thereby lined to higher losses under stress, and in other cases the allele supports the trait even more under stress. A significant effect of high Al stress could not be found as variability between seasons was high and contrasting treatments were only compared across seasons.
While drought and low P conditions significantly reduced yields, these impacts were overshadowed by cross-season and location effects that caused even more variation in yield and other traits. Late maturity is generally an advantage in higher altitudes, the opposite is observed in locations with lower altitude, overall, altitude was the major location effect observed. Soybean was found to have longer growth cycles but lower yields in higher altitudes [65] but in general few studies exist on this topic that deserves more attention.
Site regression analysis based on yield did not result in clear location clustering, unlike maturity traits that were found to be very location dependent. This suggests that even though yield correlations are more significant among locations in high or low altitudes, selection in all locations can be utilized to improve lines for all environments. Lines were identified that displayed good performance under strong stress and good cross-location yield stability, at least for the range of environments investigated here. This shows that improvement through breeding for these traits can be achieved and that selection under extreme stress does not carry a penalty for productivity under optimal conditions [3,26]. This is an important finding to plan selection in breeding programs, allowing identification of lines with good cross-location productivity, rather than requiring site specific breeding activities.

Productivity dependence on yield components and dry matter redistribution
Yield showed expected solid correlations with yield component traits and biomass traits, in line with many previous reports [11,16,62,64]. Several groups reported positive yield correlations of 100SdW in stress and non-stress conditions [14,62,66], other reports showed positive yield correlations only in stress conditions [13,15,64], occasionally even negative in non-stress and positive under stress [67]. In this population 100SdW was always positively linked to yield, seed fill and seed size are required for yield and not replaced by a larger number of smaller seed. In this work no major QTL hotspots for 100SdW appeared which are often found in such projects [50,64,68], indicating that there may be no major gene variability for 100SdW in this population.
Previous reports on yield correlations with maturity traits have been diverging. Several negative yield correlations were reported [14,64,69], suggesting that earliness aids productivity in stress as well as non-stress conditions. Other reports observed significant negative correlations only in stress conditions [11], or even contrary results between environments, with positive correlations under non-stress conditions [13,70], up to strong positive yield correlations of maturity in all conditions [62]. In this work in Pal and Qui locations DF/DH were always negatively correlated with yield. On the other hand, enhanced vegetative development contributed to increased productivity at Pop and Dar sites, hence this is rather an environment dependent effect and not due to the specific genetic background. Multi-location comparisons of effects of maturity as shown in this data set have not been published before with several major abiotic stress factors, this helps to put previous reports in perspective. Diverging correlations between yield and maturity are rather due to differences between evaluated environments than evaluated germplasm.
PHI continues to be a good proxy for yield, displaying positive yield correlations in all evaluated conditions. This confirms that effective seed fill is a useful selection criteria. Some previous publications reported correlations of PHI with seed yield in all conditions [13,14,16,62,66], while in others PHI only had significant yield correlations under stress conditions [15,64,67].
In this work the effect does not appear to be stress specific. It also appears to be useful in cooler climates in higher-altitudes.
TNC_sh proved to be an interesting trait that is linked to yield formation. TNC levels were reported to be positively correlated with starch in storage tissue of legume roots [71]. In soybean studies drought was found to reduce pod number and yield accompanied by an increase of %TNC in pods, showing that lower sink strength increases %TNC [72]. In this study positive correlations with yield, total and pod biomass, but not with vegetative plant tissue suggest two possible source-sink relationships: the active vegetative growth tissues are competing for TNCs thereby reducing yield, or the variability in TNC levels lies in the variability in photosynthetic activity of source tissues (leaves) to provide enough energy for seed production. Bush beans usually cease vegetative development at some point, which reduces sink (young) tissues, thereby increasing TNC levels and TNC uptake by reproductive sink tissues. Increased TNCs through photosynthetic activity should correlate with both increased yield and leaf and stem biomass. As the latter is not observed, the first relationship is more likely. Vegetative tissues dilute and use sugar and starch not fully committing to their source role at crucial stages and thereby deprive pods and seed from energy for generating yield. TNC_sh represents a powerful selection trait, but with current technology may be too cumbersome to phenotype to be utilized in large scale breeding.
In any environment biomass is still a prerequisite for yield, a large, strong plant is needed as a basis, hence, the resource allocation into vegetative and reproductive development needs to be balanced. Plants need to grow large quickly, and then fully commit to reproductive development without unnecessary sinks to fill seed before senescence sets in. In cool conditions at Dar and Pop sites the process is not so confined in time, and an extended seed fill phase with larger DH, not DF, aids productivity.

Major QTL affect yield and maturity
G21212 was earlier maturing than BAT 881, observed in 25 of 26 DF and DH evaluations in the 15 trials (Fig 1). QTL analysis reveals the genetic basis of this observation, which is likely the maturity QTL hotspot in chromosome 1, with 22 identified maturity QTLs over all data sets. DF1.2 BG is the strongest QTL found in this study, the allele of BAT 881 increased flowering time on average by 1.2 days. QTL for DF on chromosome 1 were previously reported showing location stable DF QTL on the start and end of chromosome 1, next to a days to harvest maturity QTL [68,73]. In the same region of DF1.2, near marker BM201 the fin locus was reported, that determines the determinate growth habit. The fin locus has been related to changes in architecture and growth period that may help the plant to tolerate stresses [74], however, both parental lines have an indeterminate growth habit. Data in this study suggests several loci (4-5) ( Table 3) with similar effects and source, for both maturity hotspots on chromosome 1 and 8. Genetically this is not very credible, the authors assume that there is only one locus per chromosome, QTL mapping software may not pick up the precise location range. In summary, one of the clearest physiological differences observed between parental lines can be explained with the most significant QTL.
The most exciting discovery of the QTL analysis is the identification of a fairly stable yield QTL. Yd4.1 BG was found in all four locations, high P, low P, as well as high Al and irrigated conditions. G21212 is higher yielding in most trials, and its allele at Yd4.1 BG explains up to 30% of variation (average of 112 kg/ha yield boost). Stable yield QTL are rare to find. A yield QTL was reported on Pv02 in two out of four trials only in non-stress conditions [50], and similarly a QTL on Pv03 was reported in two out of seven trials [66]. A condition stable QTL was published on Pv07 but observed in only two out of two trials [64]. An interesting seed yield QTL SY1.1 BR was discovered in five out of 12 trials over three locations, but not across stress treatments in any location [68]. The unusually strong and condition-independent yield QTL Yd4.1 BG found in this study is not clearly associated to another trait or a mechanism, unlike e.g. the reported QTL SY1.1 BR [68] which is a QTL for delayed flowering time. Other examples are the yield and flowering time/ growth habit QTL in tomato [75], or DRO1 QTL that modifies root architecture in rice [76]. The mechanism of action of Yd4.1 BG is unknown, but the absence of obvious physiological effects that may be undesirable for breeders suggests that it is easier to introgress into elite breeding material. Dry matter redistribution traits PHI and TNC_sh were found to correlate with yield in several trials. Three TNC_sh QTL co-localized with yield QTL, providing a mechanistic explanation for these QTL, however, the effects are not as strong and stable as Yd4.1 BG or the maturity hotspots. Direct molecular markers linked to yield have been sought after for several years with considerable investments, but their application in breeding is scarce in any crop and not reported in common bean. Yd4.1 BG represents the best candidate for use in MAS. Sequencing data from a panel will be used to identify G21212 specific SNPs tagging this locus aiming to employ it in diverse genetic backgrounds.

Summary and conclusion
This comprehensive study produced a rich source of data from 15 field trials, quantifying significant yield losses under different abiotic stresses. The study of dry matter redistribution to grain proved to be important under all conditions, and PHI and TNC_sh allowed indirect selection of higher yielding lines. Early maturity aids productivity in lower altitudes, in contrast, in higher altitudes an extended growth season is generally of advantage. Selection under strong stress conditions identifies lines that are also performing well under non-stress conditions, an important information that can be applied in selection process. QTL analysis revealed a stable yield QTL, Yd4.1 that showed significant association in all locations, next to two strong maturity QTL hotspots. Yd4.1 is the best candidate for marker assisted selection (MAS).
Drought and low P are among the major limitations affecting the production of common bean and many other crops in the tropics. Understanding which traits are associated with productivity under expected increase of drought and extreme temperature conditions due to climate change can guide research and germplasm development to increase yields in stress and non-stress conditions. Data obtained in this study will help bean breeding efforts through information on germplasm and development of molecular tools.  Table 1.