Effect of drought stress on the genetic architecture of photosynthate allocation and remobilization in pods of common bean (Phaseolus vulgaris L.), a key species for food security

Background Common bean is the most important staple grain legume for direct human consumption and nutrition. It complements major sources of carbohydrates, including cereals, root crop, or plantain, as a source of dietary proteins. It is also a significant source of vitamins and minerals like iron and zinc. To fully play its nutritional role, however, its robustness against stresses needs to be strengthened. Foremost among these is drought, which commonly affects its productivity and seed quality. Previous studies have shown that photosynthate remobilization and partitioning is one of the main mechanisms of drought tolerance and overall productivity in common bean. Results In this study, we sought to determine the inheritance of pod harvest index (PHI), a measure of the partitioning of pod biomass to seed biomass, relative to that of grain yield. We evaluated a recombinant inbred population of the cross of ICA Bunsi and SXB405, both from the Mesoamerican gene pool, to determine the effects of intermittent and terminal drought stresses on the genetic architecture of photosynthate allocation and remobilization in pods of common bean. The population was grown for two seasons, under well-watered conditions and terminal and intermittent drought stress in one year, and well-watered conditions and terminal drought stress in the second year. There was a significant effect of the water regime and year on all the traits, at both the phenotypic and QTL levels. We found nine QTLs for pod harvest index, including a major (17% of variation explained), stable QTL on linkage group Pv07. We also found eight QTLs for yield, three of which clustered with PHI QTLs, underscoring the importance of photosynthate remobilization in productivity. We also found evidence for substantial epistasis, explaining a considerable part of the variation for yield and PHI. Conclusion Our results highlight the genetic relationship between PHI and yield and confirm the role of PHI in selection of both additive and epistatic effects controlling drought tolerance. These results are a key component to strengthen the robustness of common bean against drought stresses. Electronic supplementary material The online version of this article (10.1186/s12870-019-1774-2) contains supplementary material, which is available to authorized users.


Background
Yield potential is genetically determined by the efficiencies in light interception, conversion of light to photosynthates, and translocation and partitioning of photosynthates to biomass of the harvestable organ [1,2]. For biomass partitioning to seeds, breeding efforts in major crops have achieved significant gains [3]. Nevertheless, partitioning remains relatively low in most pulses [4]. Furthermore, photosynthate partitioning and remobilization are important mechanisms of drought tolerance in common bean [5,6]. Other mechanisms, like conservative transpiration [7,8] and deep rooting [9,10] can help maintain high water status at pod filling, thus enhancing remobilization, increasing the harvest index and grain yield.
Understanding the genetic basis of yield productivity, physiological efficiency (rates of biomass and seed yield accumulation), and drought resistance mechanisms is an essential component for the design and implementation of efficient breeding strategies and marker-assisted selection for drought tolerance [11][12][13]. The timing and severity of the drought stress can significantly alter the balance between vegetative and reproductive growth [4]. For example, stress during flowering can induce flower drop and, in turn, cause more vigorous, compensatory vegetative growth, while drought stress during pod filling would cause more biomass production and decreased sink strength, resulting in a low harvest index [5]. In common bean, pod walls serve as a temporary storage of carbohydrates before rapid seed growth and efficient translocation of carbohydrates during grain filling. Thus, carbohydrate storage in pod walls and their translocation to seeds is of great importance for productivity [14].
Common bean (Phaseolus vulgaris L.; 2n = 2x = 22) is the most important legume for direct human consumption, because of its wide consumption and cultivation, and nutritional and agronomic role, which complements that of cereals and other sources of carbohydrates like root crops and plantain [15,16]. Its role as a complementary source of nutrition is limited, nevertheless, because more than half of its area of production is affected by drought, which causes yield reduction or even crop failure as the production of common bean is mostly rainfed [17,18]. Terminal (TD) and intermittent (ID) drought are two of the patterns of drought most common in bean production. The specific growth and development stage affected by drought stress, as well as its duration and intensity, are major factors affecting yield and seed quality [17,19]. In addition, one of the major effects of global warming is the change of patterns and reliability of rainfall, which will likely increase drought severity, especially near the tropics [20]. An important approach to coping with drought constraints and maintain or increase productivity is through breeding and development of resilient cultivars [21,22]. Yield, and its components, are quantitative traits under polygenic control [23][24][25]. Pod Harvest Index (PHI), defined as the partitioning of biomass from pods to seeds, has been identified as a useful and stable trait for indirect yield selection in drought and non-drought stressed environments [19,26,27]. Previous efforts in QTL discovery for PHI have identified associated genomic regions in the presence or absence of drought stress [28][29][30] or only in well-watered conditions [31,32]. However, any differential expression of PHI QTLs between terminal and intermittent drought has not been investigated to our knowledge. Furthermore, variation in PHI is determined by variation in allocation of photosynthates to pod and seed, as well as the capacity of the remobilization from the pod wall to the seed; these three factors could be differentially controlled genetically. Although the role of epistasis in common bean has been quantified previously for yield components [23] and pod shape [33,34], there are no estimates of epistasis for PHI and its components. Therefore, the objectives of the research presented here were to: 1) estimate the differential effect of intermittent and terminal drought on yield, PHI, and its components; 2) identify QTLs for these traits; and 3) estimate the role of pleiotropy-tight linkage and epistasis in controlling PHI.

Phenotypic segregations
Comparison of parents across environments ICA Bunsi and SXB405 differ significantly for pod-wall weight (PWW), whole-seed weight (average seed dry weight per pod: WSW), whole-pod weight (average dry pod weight, including seeds: WPW), 100-seed weight (SW100), but not for days to flowering (DTF), number of seeds per pod (SPP), pod harvest index PHI), and yield (YLD; Table 1). SXB405 had higher phenotypic values for all traits, except for PHI. SXB405 had 61, 59 and 60% higher weight in pod wall, seed and whole pod weight, respectively. It also had a 100-seed weight that was 10 g heavier (59%). In all traits, there was transgressive segregation among the RILs, resulting in phenotypic values that went beyond parental values (Table 1; Fig. 1). For flowering time, some genotypes were six days earlier than ICA Bunsi, and 10 days later than SXB405. For PHI, the highest line showed a value of 0.76, compared with the high parent (0.71). For yield, the highest line had 3682 kg ha − 1 across environments, which is 24% higher than SXB405. For PWW, WSW, and WPW, the transgressive segregation was less marked.

Sources of variation
Among RILs, the effects of genotype, environment and their interaction were significant for all traits at the P < 0.001 level (Table 1). In general, the differences between the two years were larger than those between treatments. For days to flowering, the effect of year was significant. Within years, the drought treatments caused significantly earlier flowering, although the average difference was less than one day. For PWW, terminal drought led to heavier pod walls in both years compared to the control. For WSW and WPW, reductions were observed in both years compared to the control. For 100-seed weight (SW100), terminal drought led surprisingly to heavier seeds than intermittent drought in the first year, and to the control in the two years. SPP was lower under intermittent drought, compared to the control and terminal drought in the first year. The second year had significantly lower number of seeds per pod, with no statistical difference among treatments. PHI was significantly lower in the well-watered control and only in the intermittent drought in the first year. In the second year, PHI under terminal drought treatment was also lower than in the control treatment. In both comparisons, however, the reduction in PHI was only 0.01. For yield, there was a reduction of 27% under terminal drought in the first year and 23% in the second year, compared to the well-watered control (Table 1).

Correlations among traits
PHI (r = 0.5), followed by WSW (r = 0.25) and SPP (r = 0.22) (Fig. 1, Table 2) were the traits most positively correlated with yield. In contrast, DTF was negatively correlated with yield (r = − 0.22). PHI had a higher (negative) correlation with the weight of pods (r = − 0.47) than to the number of seeds (r = 0.032). There was a relatively high positive correlation of the number of seeds per pod and pod harvest index (r = 0.41). SPP was correlated with the whole-pod seed weight (r = 0.38) but not significantly correlated with whole-pod weight. Correlations between yield and yield-associated traits in each year (Table 2), were similarly significant between treatments (fully irrigated control and terminal drought) within each year. WSW, WPW, SW100, SPP, and PHI were significantly correlated with yield in all year/treatment combinations. In contrast, DTF and WPW weight were only significantly (and negatively) correlated with yield in 2014. For DTF, the correlation was higher in the drought treatment than in the full irrigation treatment in both years. The correlation between yield and PHI was similar between treatments but lower in 2013 than in 2014, with 0.28, 0.26, 0.55 and 0.53 for 2013I, 2013D, 2014I and 2014D, respectively.

Molecular linkage map and QTL analyses Molecular linkage map
The length of the linkage map was 951 cM using 378 markers across the 11 linkage groups of the common bean genome (Fig. 2). There were on average 34 markers per linkage group with an average and range of spacing of 3.0 and 1.7-5.3 markers/cM, respectively (Table 3). Table 1 Descriptive statistics and analysis of variance for 8 traits: days to flowering (DTF), pod wall weight (PWW), whole-seed weight (WSW), whole-pod weight (WPW), 100-seed weight (SW100), seeds per pod (SPP), pod harvest index (PHI), and yield (YLD), in three water regimes in two year By comparing the genetic map with the physical map using shared markers with the reference genome of G19833 (Version 1.0; [35]), the map covered a physical length of 479,410,006 base pairs or 82% of the total sequence length of the 587 Mb bean genome. The coverage of the genome among linkage groups ranged from 98 to 67% of physical length. The average maximum interval without markers across linkage groups was 19 cM long. The longest interval was on linkage group Pv03 with 29 cM between flanking markers. A relatively large part of the short arm of linkage group Pv08 (~20 Mbp from the beginning of the linkage group) was nearly devoid of polymorphic SNPs. The average physical distance (kbp) per unit of genetic distance (cM) was 516 kbp/cM and ranged between 362 kbp/cM (Linkage group Pv09) and 680 Kbp/cM (Linkage group Pv10).

Identification of additive QTLs and their distribution on the molecular linkage map
Quantitative trait loci analyses were performed for each treatment and year combination as well as for the overall mean across environments, as an indication of the degree of stability of each QTL (  There were QTLs in all linkage groups. Nevertheless, only one QTL was found in linkage groups 08 and 10. Across environments, the percentage of explained variation (PVE) ranged from 24% (SPP) to 73% (DTF) ( Table 4). PHI and YLD had a percentage of explained variation of 43 and 33%, respectively. Interestingly, PHI had higher PVE in 2013 than in 2014, while for yield PVE was lower in 2013 and higher in 2014.
For days to flowering, 10 QTLs were found in eight linkage groups. Five of these QTLs were significant across all Linkage groups and QTLs of the ICA Bunsi x SXB405 RI population for pod harvest index (PHI), pod wall weight (PWW), whole seed weight (WSW), whole pod weight (WPW), 100-seed weight (SW100), seeds per pod (SPP), days to flowering (DTF) and yield (YLD). Full bars represent QTLs that were detected in all years or in the overall mean across environments. Open bars represent QTLs that were not consistent across environments. The length of the bars represent the 99% confidence interval. The color of the bars represents the positive effect of the parent, ICA Bunsi in green and SXB405 in red. The positions in the linkage group are the physical positions of the markers in Mb and plotted with MapChart [36] environments. The remaining QTLs were year-specific, and not treatment-specific. In four of the 10 QTLs, ICA Bunsi contributed the negative allele (earliness). The QTL with the largest percentage of explained phenotypic variation was located on linkage group Pv01 with 26% of the additive variance across environments, an additive effect of 1.2 days, and ICA Bunsi contributing to earliness.

Photosynthate allocation and remobilization traits
For PWW, 11 QTLs were found distributed over nine linkage groups. The ICA Bunsi allele contributed positively for only three of them. Four QTLs were expressed in all environments. The QTL explaining the largest variation at 19 cM in Pv06, explaining 29%, with SXB405 allele having the positive effect. Three QTLs were expressed only in 2013, in the control and terminal drought treatments, but not in the intermittent drought treatment. There were 10 QTLs for WSW, distributed on seven linkage groups. SXB405 was the source of the positive allele in six QTLs.
Similar to PWW, the largest QTLs for WPW were located on Pv05 and Pv06, with 28 and 11% of the variance across environments explained, and SXB405 being the source of the positive allele. For WPW, there were 12 QTLs distributed on nine linkage groups. The positive allele was contributed by SXB405 for eight QTLs. Two relatively large QTLs were found on Pv05 and Pv06, with 24 and 17% of the additive variation. In both cases, SXB405 was the source of the positive allele. These two QTLs were the only ones being expressed across environments.
For PHI, there were nine QTLs on seven linkage groups. ICA Bunsi alleles contributed positively to five QTLs. Only the QTL at the 59 cM position on Pv07 was expressed across environments. This QTL also had the highest percentage of the variation explained (17%), with SXB405 contributing the positive allele. In general, the other QTLs were expressed similarly between treatments, but not between the two years. Five QTLs were significant in all treatments in 2014, but not in 2013, while one QTL was significant across treatments in 2013 but not in 2014. The QTL on Pv01 was significant only under terminal drought in 2014. The QTL in Pv04 was significant only in the control conditions in both years.
For SPP, four QTLs were found. All QTLs were expressed in 2013, but only one in 2014. The QTL explaining the largest variation was located on Pv05 (8%), the allele of SXB405 had the positive effect. While the other three QTLs positive effect was from the ICA Bunsi. There were 10 QTLs for SW100 distributed in   eight linkage groups. SXB405, the parent with the heaviest seed, was the source of the allele with positive effect in seven QTLs. The QTL on Pv05, had the largest explained variation with 30% of the additive variation, and limited environmental influence. The SXB405 allele increased the seed weight by 1.27 g per hundred seeds over the mean.

Yield
For YLD, eight QTLs were found across seven linkage groups. SXB405 was the source of the positive allele in six of these QTLs. None of the QTLs were significant in all environments. The differential expression was mostly due to the year than to the treatment effect; none of the QTLs were significant in both years. In 2014, all the QTLs expressed were significant in both treatments. However, in 2013, two of the three QTLs were significant only under terminal drought. Nevertheless, the analysis across environments detected three QTLs, of which the major QTL in Pv06 explained 11% of the additive variation. The positive allele was conferred by SXB405, increasing yield by 106 kg ha − 1 above the average.

Genetic interactions Epistatic interactions among QTLs for the same trait
Digenic epistatic interactions were analyzed for the average of the traits across environments (Table 5) and within each environment (Additional file 2: Table S2). Across environments, epistatic interactions were found for 13 pairwise interactions, including PWW, WSW, SW100, PHI, and YLD (Table 5). Six epistatic QTL combinations were detected for PHI, two for PWW, three for WSW, and one each for WPW, SW100, and YLD. Furthermore, we found 51 epistatic interactions when analyzing each environment separately (Additional file 2: Table S2). The interactions were more consistent across treatments than across years. When analyzing the overall mean, PHI had the largest number of epistatic interactions.

Pleiotropy and epistasis across traits
We further modeled pleiotropy and epistasis jointly with CAPE (Tyler et al. 2013(Tyler et al. , 2016. Analyzing the mean phenotypes across environments, there were genomic areas that had a pleiotropic effect for all traits, especially in Pv02, Pv05, Pv06 and Pv07 (Fig. 3). In addition, we found 47 loci, which were interacting in a network involving all linkage groups. In all traits, the main effects explained most of the variance (Fig. 3). However, for SW100, PHI and YLD, the interactions accounted for 18, 25, and 16% of the variation, respectively.

Discussion
This study had three major goals: 1) to estimate the differential effect of intermittent and terminal drought on yield, PHI, and its components; 2) to identify QTLs for these traits; and 3) to estimate the role of pleiotropy-tight linkage and epistasis in controlling PHI.
Differential effect of drought type on yield, PHI, and its components The full population was grown for two years (2013 and 2014); the first year, both a terminal drought and intermittent drought treatment were imposed, while in the second year, only a terminal drought was considered. The dry conditions of the field site in Davis, California, allowed us to control the drought level through irrigation as rainfall during the summer is very rare, localized, and limited. In general, for all the traits, we found a significant effect in the treatment effects within each year, indicating that the reduced number of irrigations had imposed a sufficient experimental drought stress. While the drought stress under terminal drought was encountered after flowering, the intermittent drought treatment overlapped during mid flowering. Nevertheless, the difference between years was even more notable. Heat waves are an occasional issue and can induce inflorescence and pod abortion, as well as reduced fertilization [37][38][39]. There was a heat wave during flowering in the second year causing flower drop and a lower number of seeds per pod. Nevertheless, the yields recovered, probably through re-growth and compensation that can occur in indeterminate cultivars [38,40]. We found a significant genotype-by-environment interaction in all traits; however, pod traits (PWW, WPW, WSW, and PHI), SW100, and DTF showed high broad-sense heritability, while, unsurprisingly, SPP and YLD showed lower heritability (Table 1).
When comparing the effects of intermittent and terminal drought to the control, PHI was similar in terminal drought and the control [29,40], but was lower in intermittent drought. These results suggest that photosynthate allocation towards pods is enhanced under TD but decreased under ID. The SPP was lower in ID than in control and TD, probably because of heat stress during the reproductive stage that has an effect on fertilization and seed abortion [38,41]. Therefore, photosynthate allocation could be hampered due to lower sink strength [14].
The study of correlated traits, which arise through pleiotropy and linkage disequilibrium [42], is important because they influence the progress from selection in breeding. Across environments, the correlations of YLD with the other traits were consistent. When comparing the mean of phenotypic traits across environments, PHI had the highest correlation with YLD (r = 0.5) among traits, suggesting the importance and usefulness of the trait for indirect selection for yield [6,27]. Among its component traits, PHI had a higher absolute correlation with PWW (− 0.47) than with WSW (0.32) and SW100 (0.17), suggesting both the importance of photosynthate allocation to pods and its remobilization to the seeds. The number of seeds per pod was also correlated with PHI (0.41), suggesting that sink strength due to the number of seeds is higher than sink strength based on seed size. Individual seed weight usually has lower plasticity and higher heritability [43,44]. The number of seeds per pod is highly affected by temperature [37,38,45] and drought stress [46]. The latter is consistent with our findings that the number of seeds per pod was significantly lower in the second year, when the plants suffered a high heat episode during flowering and the differences between treatments were only significant under intermittent drought, when drought stress was present at the beginning of flowering. Interestingly, there was no correlation between SPP and PWW or SW100, in addition to the high heritability of the latter two traits, and medium correlation between pod wall and seed weight. This observation suggests that pod walls have relatively low plasticity; hence, seed and pod morphologies may be inherited differently.

Identification of QTLs Map development
A genetic map was developed with the BARCBean6K_3 SNP chip platform with an average density of 3 markers per cM. The marker density of linkage maps of other populations genotyped with the same BARCBean6K_3 SNP chip had an average of 1.61 and a range from 0.64 to 2.4 markers per cM [25,29,30,47,48]. Overall, the sample size and marker density for our analyses were adequate for QTL mapping [49,50].
However, the distribution of SNPs was uneven over the eleven linkage groups. In most of the metacentric chromosomes ( [51]: Pv02, Pv03, Pv04, Pv05, Pv07, Pv08, Pv10, and Pv11, the SNPs were distributed towards the ends of the two arms. An exception was Pv08, a metacentric chromosome, with a large region of towards one end of the chromosome. This large monomorphic area may be due to common ancestry from ICA Bunsi, as this cultivar was not only one of the parents of this RI population but has also been used as a parent of at least two of the lines (A686, SEA15) used in the development of the other parent, SXB405 (S. Beebe, unpublished information). It may also contain genes that were under the same selection during the development of the two RI parents. Traits that have been mapped to these regions include principally growth habit and phenology traits. These included number of nods and pods on the main stem, harvest index, and days to flowering and maturity [52], lodging and seed weight [53], and shoot ant total weight, hypocotyl length, root weight, and SPAD index (chlorophyll content [54]. This overall focus on plant type is also reflected in the use of ICA Bunsi as a parent to obtain white mold resistance. Miklas et al. [55] mapped several QTLs for white mold (Sclerotinia sclerotiorum Lib. de Bary) avoidance to this same region on Pv08, including canopy porosity, height, and lodging, suggesting that selection for a more upright and porous plant type is controlled at least partially by this region on Pv08 and has been selected by multiple breeding programs, whether directly for growth habit or for avoidance of white mold. Table 5 Digenic-epistatic interactions across environments for pod wall weight (PWW), whole-seed weight (WSW), whole-pod weight (WPW), 100-seed weight (SW100), pod harvest index (PHI) and yield (YLD)

QTL mapping
Additive and epistatic marker-trait association analysis for PHI, its components, and YLD were performed for each irrigation treatment/year combination as well as for the overall mean across environments. In addition, we included a QTL by environment analysis to estimate the consistency of the marker QTLs across environments. Nine QTLs were detected for PHI with the positive alleles of the QTLs contributed by both parents. A major and stable QTL was detected on linkage group Pv07 explaining 17% of the additive variation across environments. It was detected in all environments and had a low (3%) additive-by-environmental percentage of explained variation. SXB405 contributed the positive allele. No QTL in Pv07 has been reported before in previous analyses [28][29][30]32].
The other eight QTL were not detected in all environments, but were detected mostly either in one year or the other.
We found eight QTLs for yield. None of the QTLs were detected in all environments; instead, they appeared in one year or the other. However, despite the inconsistency of yield QTLs, a frequent observation among QTL studies [32,56], three QTLs were detected when mapping the mean yield across environments. The co-location of some yield and PHI QTLs across populations (and environments) is encouraging. Furthermore, the sign of the parental effect of all the QTLs across environments was similar, that is, the allele that had a positive effect in one environment did not have a negative effect in another environment.

Clustering of QTLs and pleiotropy
Pleiotropy refers to the condition when one genetic element regulates more than one gene product and trait [57]. Pleiotropy can be assessed first by the clustering of genes controlling different traits to a specific sequence although co-location does not necessarily imply pleiotropy, as it could be due to tightly linked genes. The expectation of the underlying variation and genetic control of PHI would be that the photosynthate remobilization rate, as well the intrinsic allocation to pod walls and seeds co-localize. As expected, PHI QTLs co-localized with QTLs of some of its components. However, the components were co-localized differentially. Three QTLs co-localized solely with PWW, two with the seed and WPW, one with WSW and SPP, one with SPP and DTF, and one with pod wall, seed and whole pod weight and SW100. In addition, the major QTL at Pv07 did not co-localize with any of its components, but it co-localized with a yield QTL.
This confirms that the variation in PHI is not only a measure of photosynthate remobilization capacity, but also depends on the variation in allocation of photosynthates to seeds and pod-walls. For example, the seed weight in a pod is a function of the number of seeds and the individual seed weight. The former is a factor of the number of ovules and fertilization/abortion rates (Nakamura, 1986(Nakamura, , 1988. The latter is a factor of seed volume and density. The seed volume in turn is controlled by seed dimension variation (length, width, thickness, shape, and distance between ovules within the pod). The weight of pod wall per pod would also be a factor of the pod dimension and wall thickness and density. All these traits can have unique or pleotropic genetic control as well as being constrained by the availability of photoassimilates, water, and nutrients [14,58].
Four of the eight QTLs for yield co-located with other traits. Three of them co-located with PHI, underlying the importance of photosynthate allocation and remobilization for high productivity. No co-localization between QTLs for PHI with YLD has been observed before [28][29][30]32].

Evidence of epistasis
A substantial number of digenic epistatic interactions were observed across and within environments. However, the interactions were not usually consistent among environments, especially between the two years. In addition, from the 13 pairwise interactions found in the mean across environments, only 5 were detected within years, confirming the complexity and lack of consistency of epistasis. Nevertheless, we found significant epistatic networks in all traits and linkage groups. For some traits, like PHI and yield, epistasis accounted for a large part of the variation. Our results support previous findings of digenic interactions for yield and yield component traits in beans, as well the importance of interacting genomic regions that do not show main effects [23,34]. They also suggest the importance of higher-order epistasis [23]. Furthermore, epistasis, in addition to gene complementarity, would explain the significant transgressive segregation we found in all traits [59]. As epistatic effects are substantial, especially in complex traits like yield, genomic selection in beans would benefit from models that include epistasis, which have been proven useful in improving the accuracy in other species [60][61][62][63][64]. Accumulating the right epistatic combinations is complex because of the large number of factors and, thus, the low probability of recovery of progeny with the desired gene combinations. However, improved capabilities of high-throughput genotyping and phenotyping, allows us to relatively cheaply screen large numbers of recombinant individuals with a large number of markers.

Conclusions
Through the detailed phenotypic and genetic evaluations under different irrigation regimes, we have identified a substantial number of genomic regions associated with photosynthate allocation and remobilization under drought stress in common bean. We found a significant effect of the water regime on all the traits, at both the phenotypic and QTL levels. However, the year effect was stronger. With some exceptions, when the QTL was not consistent across environments, they tended to be detected across treatments within the same year. We found several QTLs for PHI, including a major QTL on Pv07. In addition, QTLs for PHI co-localized with its component traits, but pod wall or seed allocation and remobilization patterns varied among QTLs. PHI co-localized with three out of seven QTLs for YLD, underscoring its importance as a determinant of productivity. All the QTLs had the same sign of effect between environments, strongly suggesting there are no tradeoffs among water regimes. We also found substantial evidence for epistasis, especially for yield and PHI. Although epistasis was not consistent across environments, it explained a considerable part of the phenotypic variation. Our results support an approach of joint genomic and phenotypic selection of yield and its components. PHI is a valuable selection goal in both well-watered and water-stressed environments because of its high correlation with yield and its higher heritability than yield. Information about physiological traits can be useful for improving the efficiency of selection and abiotic stress resilience [11,65]. Ultimately, more robust beans, able to better withstand yield-reducing stresses, will enhance its nutritional role as a source of dietary protein, fibers, vitamins, and minerals.

Population development
To investigate the genetic basis of PHI [defined as seed dry weight per pod at harvest or whole-seed weight (WSW) over average dry pod weight, including seeds, at harvest or whole-pod weight (WPW)] and other yield-related traits, a bi-parental population was evaluated consisting of 226 F 9 Recombinant Inbred Lines (RILs) from the cross of ICA Bunsi and SXB405, developed at the Centro Internacional de Agricultura Tropical (CIAT), Colombia. This population was selected for a QTL analysis because a previous study using a subset of this population (78 lines, F 4:6 ) found large phenotypic variation and transgressive segregation for PHI, as well as relatively high heritability and correlation with yield [27]. This study provides a more extensive analysis on the entire population; it also combines phenotypic and genotypic data, leading to an as yet unperformed QTL analysis in this population.
ICA Bunsi is a navy type bean developed in 1968 by the Instituto Colombiano Agropecuario (ICA) from a cross between Magdalena 8 and Japón 3 [66]. It has been used extensively in Canadian breeding programs [67], due to its high productivity and resistance to white mold (Sclerotinia sclerotiorum) [68,69]. ICA Bunsi carries the I gene, which confers Bean Common Mosaic Virus resistance [70]. SXB405 is a cream-colored type experimental line developed at CIAT from a four-way cross (A 686/A 774//NXB 80/SEA 15) and selected for its high productivity under drought and common bacterial blight (Xanthomonas axonopodis pv. phaseoli) resistance [27,71]. Both genotypes belong to the Mesoamerican gene pool, are photoperiod-neutral, and have a type II (indeterminate, bush) growth habit [27].

Trial design
Field experiments were carried out at the Plant Sciences Farm at the University of California, Davis (38.53°N, 121.78°W) for a total of five environments in two years and three water regime treatments. The soil type of the site belongs to the Yolo series, a member of fine-silty loam, mixed, nonacid, thermic Mollic Xerofluvents, considered well-drained, with slow to medium runoff and moderate permeability (https://soilseries.sc.egov.usda.gov).
The seeding was carried out on the 5th of June in 2013 and 8th of June in 2014. The plants were harvested on the 10th of September in 2013 and 12th of September in 2014. All water was provided by irrigation as there were no rain events during the experiments. In 2013, there were three irrigation treatments: terminal drought (TD), intermittent drought (ID), and full irrigation (well-watered). The full irrigation treatment received four irrigations while the intermittent drought received the first and third irrigation, and the terminal drought received the first and second irrigation. Only terminal drought and full irrigation treatments were carried out in 2014, receiving 4 irrigations for the full irrigation treatment, and the first two irrigations for terminal drought. The second irrigation of the terminal drought was applied during early flowering. The timing of the irrigation was decided according to weather and evapotranspiration (Fig. 4). The agricultural management was according to standard practices [72]. The parental lines and 226 of the RIL lines were planted in a split plot design with 3 blocks and one repetition per block. The experimental unit was a plot of 60 plants grown in two 3 m-long rows and 76 cm between rows (density of 131,578 plants per hectare).

Phenotyping
From each plot, 15 to 20 pods were collected at random across the plot during harvest and dried for 5 days at 50°C. Seed and pod number and pod mass were recorded. The seeds and pod walls were weighed separately. Pod harvest index was calculated as the proportion of seed biomass to the overall pod biomass (sum of pod wall + seed biomass) [17]. Days to flowering (DTF) was assessed when at least 50% of the plants in the plot had an open flower. Agronomic yield was weighed after four days of drying at 40°C to standardize grain water content. Yield was not measured on the intermittent drought treatment in 2013 as the plots were mixed after cutting due to unexpectedly high winds.

Statistical analyses
The data were analyzed in a linear mixed model: genotype, irrigation treatment and their interaction were considered as fixed effects, while the block and the interaction of block and treatment were considered as random effects. Each irrigation treatment and year was combined into one treatment with five levels. This approach was chosen since one irrigation pattern (intermittent drought) was not present in the second year. In addition, stress levels can vary between years, as weather variables and the precise timing of irrigation (in relation to growth and development of the plant) might have differential effects.
Statistical analyses were performed in the R environment [73], with the "lme4" package [74]. The "lmerTest" was used to determine the effect significance and least squared mean calculation [75], with type III hypothesis testing using the Satterthwaite approximation for degrees of freedom. R squared was obtained with the "pie-cewiseSEM" package [76]. The coefficient of variation was calculated with "sjstats" [77]. Broad sense heritability was estimated with REML [78] as GT is the variance of the genotype by treatment interaction, σ 2 E is the variance of the experimental error, e is the number of environments and r the number of repetitions per environment. The correlation between traits was calculated and plotted with the "ggally"package [79].

Genotyping
A single F 7 plant of each of the 226 RILs and the two parents were used for extracting DNA with a modified ammonium acetate based protocol [80]. The samples were genotyped with 5398 SNP markers from the BARCBean6K_3 BeadChip SNP chip platform [81] at the USDA-ARS Soybean Genomics Improvement Laboratory, Beltsville, MD. After filtering, GenomeStudio Module v1.8.4 (Illumina Inc., San Diego, CA, USA) was used for automatic SNP calling. After subsequent manual adjustments, filtering for quality control, and a 0.15 Gencall score cutoff, 5186 SNPs remained.

Genetic map construction
The SNP segregation database was filtered for marker monomorphism, missing data less than 5%, co-located markers, genetic clones and segregation distortion. The construction of the genetic linkage map markers was carried out with the "qtl" and "asMAP" packages using the Kosambi function [82,83]. From 1039 polymorphic markers between the parents, three-hundred seventyeight markers remained after filtering for missing data and co-located markers.

QTL mapping
The identification of QTLs was performed using inclusive composite interval mapping (ICIM) with QTL IciMapping version 4.1 [84]. ICIM uses a two-step strategy using stepwise regression for marker selection followed by interval mapping using phenotypes adjusted from the markers selected in the first step except the flanking markers of the current marker interval [85,86]. All the traits/environment combinations, as well as the mean value across environments, were analyzed with the ICIM-ADD for testing additive interactions. A probability of 0.001 and a step of 1 cM was used for the stepwise regression and the significance LOD threshold was calculated by 1000 permutations in each trait at a significance level of 0.01. QTLby-environment interaction and stability was studied with the ICIM-QEI mapping function, which detects average and environment specific effects [87]. To study epistasis, two-dimensional scanning was carried out in each trait/ environment combination by the ICIM-EPI method [84], with a 1 cM step and a minimum p-value and LOD score of 0.01 and 3.5, respectively.
Because of evidence of both pleiotropy and epistasis, we further modeled them through the "CAPE" package [88,89]. CAPE infers epistatic networks of QTLs affecting one or more traits First, it integrates multiple related phenotypes as eigentraits through singular value decomposition to de-correlate traits [88,90]. We used first three eigentraits, capturing 85% of the variance across traits. A single-variant scan was performed with linear regression with 1000 permutations to obtain the significance threshold. After selection of markers as covariates, a pairwise scan, for interaction detection, was carried out with 500,000 permutations to obtain empirical p values and a maximum Pearson correlation between markers of 0.7 to avoid testing highly correlated markers. The regression coefficients of single markers and marker interaction on the eigentraits were combined and reparameterized to detect between marker-pair influence independently of the eigentraits. Then, the influences were translated to the original phenotypes. The final network was obtained with a minimum adjusted p value of 0.01, and a minimum threshold power of 1. We estimated the variance explained by the QTLs as main effects and the interacting markers using linear regression. The Stepwise regression with model selection using the Akaike information criterion (AIC) in the package "MASS" [91] was used first and only the significant markers and interactions were included in the final model.

Additional files
Additional file 1: Table S1. QTLs identified for various phenotypic traits in different years and under different water availability treatments. C: well-watered control; ID: intermittent drought; TD: terminal drought. PVE: proportion of variation explained by QTL. (XLSX 28 kb) Additional file 2: Table S2