Introduction

Legumes contribute an average of 2.5 and 7.5% of total protein intake in developed and developing countries, respectively. In 28 countries belonging to the latter category, this figure is greater than 10% (Akibode 2011). The type of legumes consumed, however, varies from region to region. In sub-Saharan Africa, Latin America, and the Caribbean, dry beans (Phaseolus vulgaris L.) are the major legumes produced and consumed. In 2006–2008, worldwide dry bean production reached 15 million ha with a harvest of 10.65 million t (Akibode and Meridia 2011). In Canada, dry bean is an important rotation and cash crop, planted on 119,000 ha in 2016. Canada is also the world’s fifth-largest exporter of dry bean with an annual production of 229,000 t in 2014 (Statistics Canada 2016). In the near future, with greater emphasis on sustainable agriculture and dietary diversification, the importance of dry bean cultivation is likely to increase. Moreover, the effect of climate change on crop productivity could favor dry bean production in the northern hemisphere, including Canada (Ramirez-Cabral et al. 2016).

Cooking quality is an important factor for bean consumers worldwide. Beans are consumed after traditional cooking or in canned form. Therefore, both cooking and canning quality attributes are important for the success of bean varieties (Castellanos et al. 1997; Kelly and Cichy 2013). Traditional cooking and canning quality attributes of dry bean seed include cooking time, absence of stone seeds, hydration capacity, texture, and appearance etc. Like most legumes, bean seeds are prone to seed hardness. Seed hardness in legumes refers to the phenomenon requiring extended cooking time to allow softening to a desired texture. Hard seeds that do not imbibe any water during hydration or cooking are also known as “stone seeds”. Seed hardness affects bean seed cooking time and the hydration process in canning (Aguilera and Stanley 1985; McWatters et al. 1987). In addition, seed hardness trait has other negative impacts such as increased cost of consumption, loss of nutritional quality, canning quality, and uneven germination when seed is planted for field production (Stanley 1992).

Seed hardness is heritable trait while it is also affected by environmental factors during production and even under seed storage conditions (Argel and Paton 1999). The genetic factors affecting seed hardness are not well understood and could vary from a simply inherited trait such as ASPER (Asp) gene to an unknown number of major or minor genes. The Asp gene is associated with seed coat luster without affecting the color. Genotypes which lack Asp gene display dull/matte/opaque seed coats (Bassett 1996). Seed coat shininess has been associated with a low rate of water uptake in genotypes (Bushey et al. 2000). The impact of seed coat luster on water uptake was shown through the use of isogenic lines that differed only at the Asp locus (Konzen and Tsai 2014). However, in the same study, it was also shown that Asp was not the sole determinant of water uptake rates in black bean lines of different genetic origins. Water uptake of a variety is one of the most important criteria for consumers and the misconception that shiny seeded varieties are always poor for water imbibition has affected consumer preference. Breeding for shiny seeded varieties is beneficial on other accounts. Shiny seeds probably handle the environmental and storage stresses better than dull seeds (Diamant et al. 1989). Shiny-seeded varieties display a thicker palisade cell layer under microscope (Konzen and Tsai 2014). Shiny-seeded varieties are also known to contain more anthocyanins in their seed coats; however, this trait does not seem to have an impact on the anthocyanin content of the canned beans (Cichy et al. 2014). The Asp locus has been genetically mapped to chr 7 in several studies (Pérez-Vega et al. 2010; Cichy et al. 2014). However, the underlying gene and its mechanism have not been identified.

The environmental factors impact seed hardness through the phenomenon known as the hard-to-cook (HTC) defect. This defect occurs when legume seeds are maintained under adverse storage conditions, such as high temperature and high humidity. The mechanism of the HTC induction is still not fully understood (Liu et al. 1992). Many theories have been postulated to explain the origins of the HTC defect in legume seeds, the most documented among these is pectin-cation-phytate-phytase theory (Galiotou-Panayoutu et al. 2007; Kinyanjui et al. 2015). This theory postulates that the activity of phytase in seeds, under adverse conditions, leads to degradation of phytic acid causing the release of metal cations. These metal cations, chiefly Ca++, migrate to intercellular spaces to bind pectins and thereby rendering them as insoluble pectates. Although environment induced, the HTC phenomenon itself is not independent of genetic influences, since some varieties are more prone to HTC defect than others (Shiga et al. 2004).

It is possible that similar genetic factors are involved in seed hardness under normal conditions and hard-to-cook defect induced by adverse environment. Therefore, to develop superior varieties, it is vital to understand the underlying genetic factors that render varieties susceptible to seed hardness and/or the cause of the HTC defect. Testing for the seed hardness trait is not practical during the early breeding stages as it requires a large seed samples for effective determination. Moreover, the trait is highly influenced by environment, and therefore, multiple biological replicates are required to make an accurate assessment of the trait. Again, this is generally not feasible in early stages of breeding process. This difficulty makes seed hardness an ideal trait for molecular marker tagging and assisted breeding selection. Unfortunately, apart from Asp gene, there are no known molecular markers associated with the seed hardness trait in dry beans. However, QTLs have been reported for the calcium and magnesium content of the seed coat and water absorption (Pérez-Vega et al. 2012; Casañas et al. 2013). It is possible that these may be correlated to the incidence of stone seeds as this also agrees with the mechanism of HTC in pectin-cation-phytate-phytase theory. It is likely that more unexploited variation for this trait is present in germplasm collections. Germplasm of diverse origins should be used in QTL studies to identify genes that impact this trait. Therefore, a black bean recombinant inbred population was generated by crossing parents different for seed hardness trait.

The main objective of this study was to map the QTL for seed hardness traits and develop molecular markers for use in breeding. Multi-environmental field data was used to map three significant QTL for seed hardness and a QTL for the appearance of cooked beans. Coincidentally, several QTLs for seed weight and seed yield were also mapped.

Material and methods

Plant material and field experiment conditions

A total of 114 F2:7 RILs were derived through a single-seed decent method from a cross between black bean lines BK004-001 and H68-4. BK04-001 and H68-4 were breeding lines developed at the Morden Research and Development Centre (MRDC), Morden, Manitoba, Canada. BK04-001 showed a low incidence of stone seeds and H68-4 was identified with the highest stone seed count among the breeding lines screened in the year of selection (2010). The initial reciprocal F1 cross-pollination was made in 2011 in the greenhouse at MRDC. Seed was scarified before planting to encourage the uniform germination during each generation of the RIL advance. Randomly selected 85 RILs and two parents were grown in the field at Morden (49.1923°N, 98.0977°W, elevation 297.50 m) and Carman (49.5086°N, 98.0017°W, elevation 268 m) sites for 3 years (2014–2016) using randomized complete block design (RCBD) with three replications. Average growing season minimum and maximum temperatures in Morden were 12.4, 13.5, 13.4, 23.9, 25.7, and 24.6 °C, respectively, during the 3 years of the study. Average growing season minimum and maximum temperatures in Carman were 10.8, 11.7, 11.9, 23.6, 25.7, and 24.3 °C, respectively, during the 3 years of the study. Total precipitation for 2014, 2015, and 2016, respectively, during the growing seasons was 243.9, 171.5, and 364.6 mm in Morden and 291.7, 253, and 252.6 mm in Carman in the 3 years, respectively. However, data from the Carman trial in 2015 was excluded from analysis due to extensive flooding damage to the trial in early July.

Phenotyping

Seeds for phenotypic analyses were harvested at maturity. Seeds were harvested with a combine in 2014, but manually harvested in 2015 and 2016. Manual harvest was used to avoid any potential damage to the seed coats, which was suspected to have a major impact on the traits under study. In 2016 at the Carman site, seeds were harvested twice at the interval of 2 weeks to study the effect of timing of harvest after full maturity. Seed hardness was measured as two negatively correlated traits, namely stone seed percentage (SSP) and hydration capacity (HC) following AACC method 56-35.01 (AACC International 2007). HC was defined as the ratio of hydrated seed weight to dry seed weight. Stone seeds are the seeds, which stayed completely un-hydrated after 16 h soak in water (at 22 ± 2 °C). These traits were measured as follows: field harvested seeds were air-dried to approximately 10% moisture level before phenotyping. A random sample of 100 intact-looking seeds from each plot were weighed and soaked in water for 16 h at room temperature (22 °C). After 16 h, seeds were drained of excess water, strained, and weighed to calculate HC (ratio of hydrated seed wt. to dry wt.). Numbers of stone seeds were also counted to calculate the stone seed percentage. Three components of color, L*, a*, and b*, were measured using CM-5 Spectrophotometer (Konica Minolta INC, Japan). Parameter L* is the lightness component and can range in values between 0 (black) and 100 (white). Parameter a* varies between green to red and parameter b* from blue to yellow, with values ranging from − 120 to + 120. If two seed sample shares the same values of L*, a*, and b*, that indicates the samples will be perceived as having exactly same color by the human eye. Color loss as a result of cooking on the beans also was evaluated. Due to the lack of a canning facility, we used an alternative method for assessing cooking quality analysis. Samples from 2016 harvest were cooked for 20 min in a water bath at 96 °C. Cooking was done after a 16-h soak at room temperature. Due to the unevenness of color loss in cooked samples and presence of stone seeds, it was not feasible to get an accurate reading of seed color using spectrophotometer. Therefore, the cooked seeds were visually scored for color loss and visual appeal. A score between 0 and 5 was assigned to each sample, in which low score indicated good visual appeal and better color retention (including uniformity of color) and a high score indicates a poor visual appeal and low color retention after cooking. Two agronomic traits were also measured; total seed yield from the field plots was converted to kilograms per hectare, and weight of 100 randomly selected dry seeds from the field samples was measured to obtain 100-seed weight (SW) (g). In addition, the flowering dates were recorded as the days from planting to 50% plants with at least one flower; the maturity date was recorded as the days from planting to 90% plants matured for harvest; growth habit was recorded following the description of van Schoonhoven and Pastor-Corrales (1987).

Scanning electron microscope

Seed-coat microstructure of the parental lines was studied using the Quanta 650 FEG Environmental Scanning Electron Microscope (FEI, Hillsboro, Oregon) at the Engineering Department, University of Manitoba (Winnipeg, Canada). The seeds were sliced with a glass knife on a microtome. The sliced/split seeds and cross-sections were mounted on scanning electron microscope stubs (aluminum) with the help of a double-sided Carbon Tape. Later, the samples were coated with a thin (10–15 nm) layer of Au-Pd using a Denton Desk II Sputter Coater. The coated samples were viewed in high vacuum mode with ETD detector (Everhart Thornley Detector).

Genotyping by sequencing and linkage map construction

A genotype-by-sequencing (GBS) approach was used to generate SNP markers for the RIL population. The 114 RILs and parents were grown in the greenhouse at MRDC. Young leaves were used to isolate genomic DNA using DNeasy Plant Mini kit (QIAGEN, Valencia, CA). Isolated DNA samples were checked for quality using agarose gel electrophoresis and then sent to Genome Diversity Facility (Cornell University, Ithaca, NY) for sequencing. The GBS technique is described in detail by Elshire et al. (2011). Briefly, genomic DNA was fragmented using a type II methylation sensitive restriction endonuclease enzyme, ApeKI that recognizes the cut site CWGC. The fragmented DNA was ligated with barcoded adapters and was amplified using appropriate primers in a polymerase chain reaction (PCR). The amplified fragments were then sequenced on an Illumina Hiseq 2500 instrument (Illumina Inc., San Diego, CA). The raw sequence data were then filtered for quality and aligned with the P. vulgaris reference genome (Phytozome v11, Schmutz et al. 2014) using the Burrows-Wheeler Alignment tool (BWA V0.78-r455) integrated in the GBS analysis pipeline as described in Glaubitz et al. (2013).

A total of 80,398 single-nucleotide polymorphisms (SNPs) were identified from the 114 RILs and two parents with a minor allele frequency (MAF) > 0.01 and missing data rate per site < 90%. A total of 3115 SNPs were retained after removing those with a missing data rate per site > 10% for downstream analyses. The redundant SNPs which had strong linkage disequilibrium (LD) were further removed and only one SNP marker was retained in the same LD block. This resulted in the final 619 SNP markers for the linkage map construction. Due to significant segregation distortion from the expected 1:1 (p > 0.05), eight SNPs were also excluded from the genetic mapping. Consequently, a genetic linkage map was constructed from the data of 611 SNPs in 114 RILs (including the 85 RILs used in phenotypic analysis) using IciMapping V 4.1 (Li et al. 2007). To construct the linkage map, SNP markers were grouped based on LOD score of 4.0 and ordered using REcombination Counting and ORDering (RECORD) and COUNT algorithm (Van Os et al. 2005). Genetic map distance was estimated in centi-Morgan (cM) based on Kosambi mapping function. Linkage groups were oriented and assigned to the chromosomes using anchoring markers from the P. vulgaris consensus map (Galeano et al. 2011), and matching SNP coordinates from GBS data with the genome sequence information (DOEJGI, www.phytozome.net). A few SNP markers closely linked to seed hardness traits were also converted to dCAPS markers (Neff et al. 1998) for marker-assisted selection purposes. A graphical representation of genetic map was constructed using MapChart (Voorrips 2002). For the QTL effect and interaction plots, the phenotypic values were plotted against the genotypes at one and two linked markers, respectively. The QTL effect and interaction plots were generated using “plotPXG” and “effectplot” function of the “qtl” package in R (Broman et al. 2003).

Field experiment data analyses

To reduce the heterogeneity of variance and meet the normality assumption, stone seed percentage data were converted to proportion and transformed using log transformation. The data for other traits were used in their original form as they met the assumptions of normality (W-test stat > 0.90). Analysis of variance for each year and site was performed separately. For calculation of heritability parameter, all effects were considered random. The basic model used was as follows:

$$ {Y}_{ij}=\upmu +{\beta}_i+{\tau}_j+{\varepsilon}_{ij} $$

where μ is the mean, β i is the block effect of ith block, and τ j is the treatment effect of the jth treatment. ε ij is the error, following a normal distribution N(0, \( {\sigma}_e^2 \)). \( {\sigma}_e^2 \) is the error variance. Broad-sense heritability was calculated on the entry mean basis as follows:

$$ H=\frac{\sigma_g^2}{\sigma_g^2+\frac{\sigma_e^2}{r}} $$

where σ2 g and σ2 e are the genetic and error variances, respectively, and r is the number of replications. Statistical calculations were performed using the META-R (Multi-Environment Trial Analysis with R for windows) (Alvarado et al. 2015). Best linear unbiased estimates (BLUEs) were calculated for various traits using META-R. BLUEs were calculated using ordinary mean squares and considering all effects as fixed. Combined data from Morden and Carman sites in 2016 were used to calculate Pearson’s correlation coefficients between traits using “rcorr” function in “Hmisc” package in R (R Development core team 2008).

QTL analysis

Quantitative trait loci (QTL) analysis was performed on BLUEs for all traits. For the stone seed trait, BLUEs from both untransformed and transformed data were used for QTL mapping. The results were similar; therefore, only the results from the original untransformed data were adopted. For QTL mapping, inclusive composite interval mapping (ICIM), a mapping function that considers all markers simultaneously to compute stepwise regression of the markers was used (Li et al. 2007). To test the statistical significance of QTL candidates, the logarithm of odds (LOD) score was estimated through 1000 permutation tests (Churchill and Doerge 1994) and a type 1 error rate of α ≤ 0.05. Only QTL with LOD score above the threshold was retained. The genotypic variation explained by each QTL was calculated from the ratio: R2/H where R2 is the percent phenotypic variation explained by the QTL and H is the broad-sense heritability estimate for the trait.

Results

The parental lines, BK04-001 and H68-4 of the RIL, had similar growth habits (type I) and required the same number of days to flowering (48–51) and to reach maturity (95–100). However, they displayed significant differences in seed hardness traits (Table 1). Two parents differed in their seed coat lustres, with shiny for H68-4 and opaque/matte for BK04-001. Their seed hardness traits were also significantly (P < 0.05) different in all environments while they did not produce significant differences in agronomic or color traits in all the cases (Table 1). Among the RILs, significant (P < 0.05) differences were observed for all traits. Transgressive segregation was also observed for all the traits, indicating that the parents were genetically diverse for these traits (Table 1). Seed hardness was measured in terms of SSP and HC. Although there was an inverse relationship between these two traits, they were also complementary to each other as HC measurement took into account all seeds, including fully and partially hydrated while SSP was only based on completely un-hydrated seeds (Supplementary Fig. 1). In 2014, when trials were harvested with a combine harvester, the SSP was lower than that in years 2015–2016. As indicated from the SSP means of the RILs, SSP values were slightly higher in the Carman (CA) trials as compared to Morden (MO) traits (Table 1). In terms of harvesting time in 2016, the seeds from the second harvest had fewer stone seeds than those from the first, although their SSP range and heritability estimates were similar (Table 1).

Table 1 Means and ranges of seed traits for the parents and RILs, and broad-sense heritability estimates in RILs based on three replications grown at two locations during 2014–2016 in Manitoba, Canada
Full size table

As expected, a similar trend was evident for HC, given the inverse and proportional relation of this trait with SSP. Heritability estimates for SSP and HC were greater than 90% in all environments, indicating strong genetic basis for the observed variability. Color parameters were only measured in 2016 in seed samples from the Morden trial and the first harvest of Carman trial (CA16). Among the color parameters, L* had the lowest heritability, and a* had the highest (Table 1). There were also significant differences among sites for heritability of color parameters. Among agronomic traits, the heritability of 100 seed weight (SW) was higher than that of seed yield (SY).

Pearson’ correlation coefficients indicate that SSP was significantly correlated with all traits except L* (Table 2). HC was also correlated to all traits except L* and visual appearance score (VSC). L* was not correlated with a*; however, both L* and a* were significantly correlated with b*. As expected, SY was significantly correlated with SW. VSC was correlated with SSP, a*, and SW. These correlations, positive or negative, suggest linkages of corresponding trait QTL on the genetic map.

Table 2 Pearson correlation coefficients between various traits in the RIL population
Full size table

The genetic map was constructed based on the mapping data generated from 114 RILs and 611 SNP markers (Fig. 1). One phenotypic marker, Asp (asper) based on seed coat luster (Lamprecht 1940), was also scored. In addition, some InDel markers (Moghaddam et al. 2014) and derived cleaved amplified polymorphic sequence (dCAPS) markers were also developed to fill in the gaps between two SNP markers to reduce the size of intervals surrounding certain QTL. These SNP and dCAPS markers were mapped to a total of 114 bins. Total map length was 1024 cM with markers present on all 11 chromosomes (Fig. 1 shows only chromosomes in which QTL were mapped). A total of 27 unique QTL for eight traits were mapped in this study (Fig. 1). Table 3 provides the summary of QTL analysis results.

Fig. 1
figure 1

Linkage map of identified QTL for all traits. QTL are depicted left of the chromosomes with solid bars indicating 1-LOD interval and outer whiskers indicating 2-LOD interval. The QTL labels are derived by joining trait name abbreviation and site with an underscore

Full size image
Table 3 Chromosome locations and effects of all significant QTL discovered in this study using the RIL population derived from BK04-001/H68-4 and phentoyped at two sites during 2014–2016 in Manitoba, Canada
Full size table

QTL affecting seed hardness traits, SSP and HC, were mapped onto chromosomes (chrs) 1, 2, 5, 7, 8, and 10. Among these, QTL on chrs 1, 2, and 7 were detected in all environments except in Morden 2014. In contrast, QTL on chrs 5, 8, and 10 were only detected in single environments. The QTL on chr 7 was mapped in the vicinity of Asp gene, as reported by Cichy et al. (2014), with an interval of 326.6 kb. Phenotypic variations explained (PVE %) by the major QTL on chr 7 for SSP ranged from 28.5 to 43.9% in different trials (Table 3). The Asp phenotypic marker and left-flanking marker NDSUInd07c1410 also shared the same bin, suggesting a tight linkage between seed coat luster and Asp for the seed hardness. The QTL for SSP and HC on chr 1 resided at the distal lower arm, with an interval (smallest) of 13.6 kb. The PVE explained by this QTL ranged from 9.6 to 23.6%. The location of QTL on chr 2 was located at the upper arm with an interval of 2 Mbp, and the PVE accounted for by the QTL ranged from 10.2 to 19.3%, respectively. QTLs on chromosomes 1 and 7 were contributed by the hard-seeded H68-4 parent and the QTL on chr 2 was contributed by the soft-seeded parent BK04-001. In addition, QTL for SW were detected on chrs 1, 3, 8, and 10. QTL for SW also were contributed by both parents, as indicated by transgression in the range of values (Table 1). QTL increasing SW on chr 3 was contributed by H68-4, and QTL on chrs 8 and 10 were contributed by BK04-001. In contrast, all three SY QTL for increasing yield were contributed by H68-4. Interestingly, a major VSC QTL that accounted for 28.3% of PVE was detected on chr 4 in the 2016 Carman trial. Similar VSC QTL peaks were also detected in Morden 2016 trial, but their LOD scores were below the detectable threshold. QTL for color parameters, b*, were also detected on chrs 4 and 7. Figure 2 shows graphically the QTL effects and interactions for the 2016 Carman test for three QTL affecting SSP trait and the major QTL on chr 4 affecting VSC. The effects and patterns were similar for the Morden 2016 test. Supplementary Fig. 2 reports the QTL effects and their interactions for data collected from the second harvest of the 2016 Carman test. As seen in the figure, the QTL effects were not affected by the harvest time.

Fig. 2
figure 2

QTL effects (a–d) and interactions (e–h) estimated from the first-harvest seeds at Carman site in 2016. Error bars indicate SE

Full size image

In order to determine the role of seed coat on the production of stone seeds, scanning electron microscope (SEM) analysis of BK04-001 and H68-4 seeds was conducted (Supplementary Fig. 3). The shiny seeds of H68-4 displayed smoother seed coat surface as compared to the BK04-001 seeds which had a rougher surface. No prominent wax layer was observed on the seed coat of either parent. Thickness of the palisade layers was measured and there were no significant differences among the parents (BK04-001 = 29.69 μm ± 0.31 SE, H68-1 = 28.66 μm ± 0.28 SE; P = 0.1498).

Discussion

Stone seeds have been identified as a major production issue under short season growing conditions in the bean breeding program at MRDC. Stone seeds were a serious concern in all bean market classes, especially black beans. Furthermore, the trait was also shown to be highly affected by environment, hampering selection efficiency in breeding. The aim of this multi-year study was to find stable genetic factors mainly affecting seed hardness, and therefore cooking quality in black beans. Using a RIL population, several QTL were mapped for seed quality traits as potential targets for selection in breeding programs. The mean number of stone seeds was affected by factors such as year, location and harvesting time (Table 1). Mean SSP values were always higher at the Carman site, as was SW. Although effect of harvesting time was only studied in 2016 at one location, those findings match the general experience over the years at MRDC. This indicates a strong role of the local environment at maturity in generating stone seeds. However, the heritability estimates and QTL analyses indicate that regardless of sites, the genetic control of SSP is still stable and strong (Tables 1 and 3).

There was a small but significant positive correlation of SSP with SY and SW in 2016. This was also evident in the QTL mapping results, where SY and SW QTL are co-localized with SSP QTL on chr 7 (Fig. 1). This indicates a positive association between these traits due to close linkage or pleiotropy, which requires further investigation. The correlation of SSP and HC with VSC was 0.13 and − 0.087, respectively. Only one out of three VSC QTL was co-localized with SS and HC QTLs (Fig. 1). The VSC QTL on chr 7 linked with Asp explained 13.7% of PVE indicating that seed coat luster had only a small effect on the cooked bean appearance as previously reported by Cichy et al. (2014).

QTL analysis indicated that a gene at/near the Asp locus plays a major role in stone seed production and hydration capacity. A QTL at this location was also identified for HC by Cichy et al. (2014). However, no QTL at this location was identified by Pérez-Vega et al. (2010) even though the RIL parents were segregating for seed coat luster. This suggests that a closely linked gene to Asp may be responsible for the effect on HC. In addition, two novel and stable QTL were identified in the RIL population for SSP and HC. For SSP, the QTL on chrs 1 and 2 together explained more of the variation than the major QTL on chr 7 (Table 3 and Fig. 2a–c). The effect of the chr 2 QTL was completely additive to the QTL on chr 7 (Fig. 2f). However, the effect of chr 1 QTL was dependent on presence of the chr 7 QTL (Fig. 2e). As seen in the Supplementary Fig. 2, the QTL effects were stable; however, there were some shifts in the interaction patterns between the QTL. The VSC QTL, however, remained independent of Asp.

A direct role of Asp for the seed hardness trait cannot be ruled out. It was reported that homozygous lines for the recessive asp had a rougher cell surface (Beninger and Hosfield 2000; Konzen and Tsai 2014). Theoretically, a rougher seed coat would increase the contact surface and therefore increase water uptake rate during the soaking treatment. The SEM analysis in this study also confirmed differences in surface of seed coats between shiny and matte seeds as previously reported (Beninger and Hosfield 2000; Konzen and Tsai 2014). Those studies also reported differences in thickness of palisade cell layer between shiny and matte lines. However, this was not observed in the current study likely due to the use of non-isogenic lines for SEM analysis. The SEM analysis of BK04-001 and H68-4 parents also did not indicate presence of epicuticular wax layer in either line (Supplementary Fig. 3). Based on the available genomic sequence, both Asp and chr 7 QTL were mapped over a 326.6 kbp region. This region contains 41 genes (P. vulgaris genome V1.0, www.legnumeinfo.org). In future, the genetic and phenotypic characterization of recombinants in the Asp region can be used to determine if the gene underlying Asp phenotype and the major QTL are the same. Breeding programs use mostly an elite and narrow genetic base. This usually results in low genetic diversity between the parents resulting in uneven marker density during genetic mapping. The QTL region on chr 2 spans multiple Mbs due to a lack of SNP markers in that region. The QTL region of chr 1 is small enough to suggest some candidate genes. The minimum QTL interval (mapped in 2016 Carman test) contains only the two genes, Phvul.001G264300 and Phvul.001G264400. In Arabidopsis, Phvul.001G264300 is homologous to AT2G32000.1 (DNA topoisomerase, type IA, core) and Phvul.001G264400 is homologous to AT2G31980.1 (PHYTOCYSTATIN 2). Among those two genes, Phvul.001G264400 shows exclusive and high level of expression in seed tissues (http://plantgrn.noble.org/PvGEA/SearchVisual.jsp, O’Rourke et al. 2014). Gene sequencing from both parents revealed a mutation in the 3′ end of the gene (Supplementary Fig. 4). This gene is therefore a good candidate for the chr 1 QTL.

In conclusion, the seed hardness trait in black beans is an oligogenic trait that is significantly influenced by the environment. Screening diverse germplasm can help identify the factors affecting seed quality traits. While the environmental factors that impact stone seed production are still unknown, the incidence of the stone seed trait can be improved by using marker-assisted selection. In addition, the cooked appearance of black beans was largely inherited independently from the seed coat luster (Asp) and the QTL for hydration capacity. This suggests only a limited trade-off between various seed quality traits.