Plant materials and DNA extraction

For this study, a cross was made between ‘Aragorn’ (PI 648006) and ‘Kiflica’ (PI 357292). ‘Aragorn’ is an agronomically desirable and widely grown variety with a low to medium concentration of mineral nutrients, while ‘Kiflica’ is a variety with a high concentration of mineral nutrients and less desirable agronomic characteristics [42]. Aragorn seed was provided by Plant Research (NZ) Ltd. Kiflica seed originally collected in Macedonia and donated to the USDA Western Regional Plant Introduction Station by Aladzajkov Lazar in 1970. Kiflica is freely available from the USDA (https://www.ars-grin.gov). The cross was made in Pullman, WA in 2010 and single seed descent was used to get a F₆ generation consisting of 158 recombinant inbred lines (RILs).

Fresh leaf tissue from each RIL was ground using a Geno/Grinder 2000 (SPEX SamplePrep, Metuchen, NJ) and total DNA were extracted using DNeasy 96 Plant Kit (QIAGEN, Valencia, CA). A NanoDrop ND-1000 spectrophotometer was used to quantify the DNA concentration of each extracted sample following the manufacturer’s instructions (Nano-Drop Technologies, Wilmington, DE).

SSR and allele specific marker analysis

A total of 114 simple sequence repeat (SSR) primer pairs from the work of Loridon et al. [43] and one eIF4E allele specific marker from the work of Smýkal et al. [44] were chosen to anchor this study’s linkage map to the SSR-based map of Loridon et al. [43]. PCR amplifications were performed with 4 ng genomic DNA, 1 × PCR buffer, 1.5 mM MgCl₂, 0.2 mM dNTPs, 0.05 μM forward primer, 0.25 μM reverse primer, 0.2 μM M13 primers with dyes of FAM, VIC, NED, PET, 0.6 U BIOLASE^TM DNA polymerase (Bioline), and 6.76 μl ddH₂O in a total volume of 12 μl. The cycling conditions included initial denaturization at 95 °C for 5 min, followed by 42 cycles, each of which consisted of 95 °C for 1 min, 56 °C for 1 min, and 72 °C for 1 min. The final extension was at 72 °C for 10 min. These PCR products were analyzed on an ABI 3730 DNA analyzer (Applied Biosystems) and data were scored using GeneMarker software version 2.2.0 (SoftGenetics).

SNP markers analysis

Two hundred fifty four gene-based SNP markers from the work of Deulvot et al. [13] were selected for use in this study to anchor this linkage map to the previously published gene-based map. The SNP genotyping was analyzed using the Sequenom MassARRAY iPLEX platform [45]. Eight iPLEX assays, each carrying 28–36 SNP markers, were developed with the software SpectroDesinger v3.0 [39].

The iPLEX GOLD reactions consisted of three parts: the iPLEX PCR reaction, the SAP reaction, and the iPLEX Extend reaction. The iPLEX PCR amplifications were performed with 25 ng DNA sample, 1 × PCR buffer, 2 mM MgCl₂, 0.5 mM dNTPs, 0.1 μM each PCR primer, 1 U Taq DNA polymerase (Bioline), and 1.8 μl ddH₂O in a total volume of 5 μl. The reaction was performed at 95 °C for 2 min, followed by 45 cycles, each of 95 °C for 30 s, 56 °C for 30 s, and 72 °C for 1 min. The final extension of 72 °C was for 5 min. After the iPLEX PCR, the SAP reaction was performed with 0.17 μl 10 × SAP buffer and 0.5 U SAP enzyme. Then, samples were incubated at 37 °C for 40 min, followed by 85 °C for 5 min. The iPLEX Extend reaction was performed with 1 × iPLEX buffer, iPLEX terminator, a primer mix containing extension primers with a final concentration between 0.625 and 1.5 μM, and 1.35 U iPLEX enzyme. The amplification conditions were performed as follows: 95 °C for 30 s; followed by 40 cycles, each of which consisted of 94 °C for 5 s followed by 5 cycles of 52 °C for 5 s and 80 °C for 5 s; and a final extension at 72 °C for 3 min. Then, 6 mg of resin was added in each well. The iPLEX extension products were dispensed on a SpectroCHIP through a RS1000 Nanodispenser (Sequenom). A matrix-assisted laser desorption/ionization time-of-flight mass spectrometry (MALDI-TOF) mass spectrometer (Sequenom) was then used for the SNP genotyping.

GBS library construction and SNP identification

The DNA of the 158 lines of the RIL population and the two parents were used to construct GBS libraries [7]. The concentration of genomic DNA was 100 ng/ul. ApeKI which recognizes GCWGC (where W = A or T) was used as the restriction enzyme. The libraries were sequenced on an Illumina HiSeq 2000 at the Cornell University Genomics Core Laboratory.

The raw data were analyzed with the universal network enabled analysis kit (UNEAK) pipeline, which was developed for non-reference GBS SNP calling [11]. In this pipeline, the following parameters were used: minimum number of tags was five, error tolerance rate was 0.03, minor allele frequency (MAF) was 0.4, and the sample calling rate was 0.5. The SNPs with unknown or heterozygous genotypes in one or two parents were also removed. Finally, using an in-house perl script, homozygotes for alleles of ‘Agarorn’ were recorded as “A”, homozygotes for alleles of ‘Kiflica’ were recorded as “B”, and heterozygotes were recorded as “H”.

Linkage map construction

Polymorphic markers from published maps [13, 43] and SNPs from this GBS study with less than 20% missing data per sample were used to construct the linkage map. The genetic linkage map was constructed using OneMap software [46] with LOD values arranged from 3 to 6 and a recombination frequency less than 0.3. The Kosambi mapping function was used to calculate centimorgan distances. The recombination counting and ordering (RECORD) algorithms were used for ordering the markers [47].

Comparative mapping

Comparative mapping was performed between pea and genetically close legumes. The mapped pea SNP sequences from this study were pairwise aligned using the BLAST algorithm (BLASTN, E value < 1e-10, percentage similarity > = 90%) with the Medicago truncatula genome v4.0 (http://www.jcvi.org/medicago/, [48]), the Cicer arietinum genome v1.0 (http://cicar.comparative-legumes.org/, [49]), the Phaseolus vulgaris genome v1.0 (http://phytozome.jgi.doe.gov, [50]) and the Lotus japonicus genome v2.5 (ftp://ftp.kazusa.or.jp/pub/lotus/lotus_r2.5/pseudomolecule/, [51]). The software Circos [52] was used to visualize the synteny between the closely related species. The cM distances on the pea linkage groups were multiplied by 250,000 to match the bp on the chromosomes of the four genomes.

Phenotyping

The field experiments were established in two locations, Whitlow (46°74’N 117°13’W) and Spillman (46°43’N 117°10’W) farms in Pullman WA in 2014. Plots in each environment were planted in a randomized complete block design with three replications. The seeds from each plot were harvested from pods, cleaned and dried at room temperature. One hundred seeds from each plot were weighed. To analyze mineral concentrations, 20 g of seeds were ground using stainless-steel coffee grinders to homogenize each sample. Then, 0.5 g sub-samples were digested to dryness using concentrated ultra-pure nitric acid and hydrogen peroxide as previously described [53]. Subsequently, the digestates were resuspended in 2% nitric acid. Each sample was analyzed for ten different elements, B, Ca, Fe, K, Mg, Mn, Mo, P, S and Zn, using inductively coupled plasma optical emission spectroscopy (ICP-OES; CIROS ICP Model FCE12; Spectro, Kleve, Germany). The instrument was calibrated with certified standards each day and blanks and certified tissue standards were run to verify the accuracy of the instrument. Mineral content per seed was calculated by multiplying the average sample elemental concentration by average seed weight.

Statistical and QTL analysis

All the trait data from tissue analysis under the two different environments were analyzed by analysis of variance (ANOVA) using the SAS program PROC MIXED (SAS Institute Inc., Cary, North Carolina, USA). Pearson’s correlations among the quantitative traits were calculated using the SAS program PROC CORR. The broad-sense heritability (H ² ) was calculated for each trait as H ²   = σ_G ²/[σ_G ² + (σ_GE ²/e) + σ_e ²/re], where σ_G ² = genotypic variance, σ_GE ² = variance due to interaction between genotype and environment, σ_e ² = error variance, e = number of environments, r = number of replicates. The variances were calculated by SAS program PROC MIXED with genotypes and environments considered as random effects. The QTL Cartographer 2.5 software [54] was used to identify and locate QTLs using the composite interval mapping (CIM) method by the permutation test (1000 times) at a P value of 0.05. The backward regression model was used to get cofactor and walk speed and window size were set to 1 cM and 10 cM respectively. The LOD score threshold for detecting QTLs was set between 3.0 and 3.2 for all the traits. Mapchart (version 2.2) [55] software was used to draw the genetic linkage map and the QTLs.

Polymorphism analysis

A total of 40 out of 114 SSR markers from the map of Loridon et al. [43] showed polymorphism between the two parents, ‘Aragorn’ and ‘Kiflica’. Of these, 27 of the SSR markers were successfully amplified and scored in the RIL population. The 254 Sequenom-designable SNP markers from the map of Deulvot et al. [13] were screened among the RILs, along with eight iPLEX assays, of which 50 were polymorphic. SNPs with more than 20% missing data were removed from further study. Finally, 47 informative SNPs were used for linkage map construction. The eIF4E allele specific marker was also polymorphic and used in linkage map construction. The sequences of the polymorphic markers are listed in Additional file 1.

SNP discovery

Linkage mapping

The 1609 SNPs identified from GBS were combined with the polymorphic SNP and SSR markers from previous linkage maps to give a set of 1684 markers for constructing the linkage map. Of these markers, only one SNP, TP56850, could not be assigned to a linkage group. 1683 markers were assigned to seven linkage groups with the identity and orientation of the linkage groups determined by the 75 previously mapped markers (Fig. 1). The markers were evenly distributed throughout the seven linkage groups with 99% of the intervals between the adjacent markers being smaller than 10 cM. The estimated map length was 1310.1 cM and the map had a density of 1.3 markers per cM (Table 2).

	LGI	LGII	LGIII	LGIV	LGV	LGVI	LGVII	Total
Number of markers	201	267	261	235	236	229	254	1683
Length (cM)	158	178.1	226.6	168.8	190.3	176.6	211.7	1310.1
Number of marker per cM	1.3	1.5	1.1	1.4	1.2	1.3	1.2	1.3
Number of gaps (>10 cM)	2	3	5	1	2	2	4	19

Comparative mapping

1608 corresponding DNA sequences from the mapped SNP loci were used for comparative genome analysis to evaluate syntenic relationships between pea and other closely related legumes. Comparison between pea and M. truncatula showed the closest genetic relationship (402 sequence matches) (Additional file 2). Pea linkage groups PsLG I and PsLG V were syntenic with M. truncatula chromosomes MtChr 5 and MtChr 7, respectively. PsLG II exhibited synteny with MtChr 1 with a large inversion. Some pea linkage groups were collinear with more than one M. truncatula chromosomes: PsLG III - MtChr 2 and MtChr 3; PsLG IV - MtChr 4 and MtChr 8; PsLG VI - MtChr 2 and MtChr 6; PsLG VII - MtChr 4 and MtChr 8 (Fig. 2a). In the case of pea and chickpea, 296 sequence matches were observed. Among pea linkage groups, PsLG II, PsLG IV, PsLG V and PsLG VII were collinear with C. arietinum chromosomes CaChr 4, CaChr 7, CaChr 3 and CaChr 6 respectively. Also, three pea linkage groups, PsLG I, PsLG III and PsLG VI, showed syntenic relationships with CaChr 2 and CaChr 8, CaChr 1 and CaChr 5, CaChr1 and CaChr 8 (Fig. 2b). Although there were 91 sequence matches (Fig. 2c) between pea and P. vulgaris, and 86 sequence matches between pea and L. japonicus (Fig. 2d), limited syntenic patterns were observed between these genomes.

Phenotypic analysis

Correlation coefficient analysis was performed between the seed traits among all the RILs in both locations (Additional file 5). Positive correlations were observed between all the seed mineral concentrations and between all the seed mineral contents. Negative correlations were shown between seed weight and all the mineral nutrient concentrations. The highest positive correlations between different mineral concentrations were observed between Ca and Mn (0.69), Mg and Mn (0.69) and the highest positive correlation between different mineral contents was seen between Fe and S (0.85), while the lowest positive correlations between mineral concentrations were between K and Mn (0.13), P and Mo (0.13), Ca and Zn (0.13) and the lowest positive correlations between mineral contents were between Mo and P (0.15), Mo and S (0.15), and Mo and Zn (0.15).

QTL analysis

Five QTLs were identified for B concentration, two of which were detected in both locations. It is worthwhile to note that the QTL, [B]-Ps5.1, explained 42% of the phenotypic variance and had ‘Kiflica’ as the contributing parental allele. Three B content QTLs were observed with explained variances of 15.1%, 10.3%, and 11.3% respectively, two of which had ‘Kiflica’ as the contributing parental allele. Several B concentration QTLs and B content QTLs colocalized with each other: [B]-Ps1.1 colocalized with B-Ps1.1, [B]-Ps6.1 colocalized with B-Ps6.1.

Five Ca concentration QTLs were detected, with the one on LG V, contributed by ‘Kiflica’, explaining 31% of the phenotypic variance. Four Ca content QTLs were identified with all of them observed in both locations. Three Ca concentration QTLs, [Ca]-Ps4.1, [Ca]-Ps5.1, [Ca]-Ps7.1, overlapped with the Ca content QTLs, Ca-Ps4.1, Ca-Ps5.1, Ca-Ps7.2, respectively.

Five Fe concentration QTLs were identified and three of them were observed in both locations. One QTL, [Fe]-Ps7.1, explained 19.4% of the phenotypic variance with ‘Kiflica’ as the contributing parent. There were five Fe content QTLs detected, four of which had ‘Aragorn’ as the contributing allele. Co-localizations were observed between [Fe]-Ps2.1 and Fe-Ps2.1, [Fe]-Ps5.1 and Fe-Ps5.1.

Six K concentration QTLs were identified, three of which were observed in two environments. It is noteworthy that the QTL [K]-Ps5.1 explained 43% of the phenotypic variance and had ‘Kiflica’ as the contributing allele. Four K content QTLs were observed but none of them was identified in both locations.

Four Mg concentration QTLs and four Mg content QTLs were identified. The concentration QTL [Mg]-Ps5.1 explained 43.3% of the phenotypic variance and had ‘Kiflica’ as the contributing parent.

Five Mn concentration QTLs were detected in both locations, one of which, [Mn]-Ps5.1, explained 29.9% of the phenotypic variance. ‘Kiflica’ was the contributing parent for the QTLs on LG II and LG V, while ‘Aragorn’ was the contributing parent for the QTLs on LG I, LG IV and LG VII. Three Mn content QTLs were identified and one with an explained variance of 19.8% was identified on LG VII. Co-localizations were observed between [Mn]-Ps7.1 and Mn-Ps7.1.

The Mo concentration QTL [Mo]-Ps5.1 on LG V was identified in both locations and it explained 34.2% of the phenotypic variance. The Mo content QTL, Mo-Ps5.1, was found with explained the variance of 24.6% and colocalized with the Mo concentration QTL.

Five P concentration QTLs and six P content QTLs were identified. One P content QTL explained 20.5% of the phenotypic variance. Four P concentration QTLs, [P]-Ps3.1, [P]-Ps5.1, [P]-Ps7.1, [P]-Ps7.2, overlapped with the P content QTLs, P-Ps3.1, P-Ps5.1, P-Ps7.1, P-Ps7.2, respectively.

Five QTLs were identified for S concentration, two of which were observed in both locations. All the QTLs for S concentration had ‘Kiflica’ as the contributing parental allele. Three QTLs were detected for S content, one of which explained 14.1% of the variance.

Five Zn concentration QTLs were identified. All the QTLs except the one on LG II had ‘Kiflica’ as the contributing parental allele. Four Zn content QTLs were detected but none of them were identified in both locations. Co-localizations were observed between [Zn]-Ps3.1 and Zn-Ps3.1, and [Zn]-Ps5.1 and Zn-Ps5.1.

For 100-seed weight, six QTLs were identified, one of which in LG V explained 23.2% of the phenotypic variance and had ‘Aragorn’ as the contributing parental allele. The co-localizations between seed weight and different nutrient traits were observed on LG III, LG V, LG VI and LG VII. The most common co-localization was found on LG V for most of the traits (B, Ca, Fe, K, Mg, Mn, Mo, P, S, Zn concentrations, Ca, Fe, Mo, P, Zn contents and seed weight).

SNP discovery using GBS

GBS has been shown to be an effective method that allows simultaneous discovery and genotyping of a large number of novel SNPs. The first SNP discovery study using GBS in pea was conducted by Boutet et al. [20] with a total of 64,263 SNPs mapped in a pea genetic map, indicating the GBS approach could significantly facilitate genetic studies and improvement of genomic resources. In our study, a set of 1609 high-quality SNPs were discovered between the two parental lines and well distributed throughout the seven linkage groups of pea. The complexity reduction involving the methylation-sensitive restriction enzyme ApeKI was adapted in this study in order to avoid and lower repetitive regions of the genome. The variation of reads per sample observed in some other reported GBS studies [7, 57] was also found in this study. These likely resulted from differences in the DNA quality and quantity among different samples. In this study we provide evidence that ApeKI is indeed well-suited for GBS in pea. Alternate restriction enzymes [8] and modified GBS library preparation [58] will provide a further set of tools to customize GBS for SNP discovery in pea. The entire process from DNA extraction, library preparation to next-generation sequencing and SNP calling through the UNEAK pipeline was relatively simple and fast. In addition, the overall cost of GBS was economically efficient. Given the benefits of using GBS as shown here and in other studies [8–11], it may be expected that utilization of this technology will become widely adopted in pea research for marker discovery and application in marker-assisted breeding.

Genetic linkage mapping

So far, significant efforts have been put into pea genetic linkage map construction and many genetic linkage maps have been successively developed through different types of molecular marker technologies. In this study, the ‘Aragorn’ × ‘Kiflica’ linkage map exhibited a length of 1310.1 cM and an average marker density of 1.3 per cM. It included 1608 SNPs and the 75 anchor markers which were used to determine the orientation of the seven linkage groups. The anchor markers provided a bridge that allowed us to combine our linkage map with the recently published pea consensus map [19]. Comparison of the anchor markers in the ‘Aragorn’ × ‘Kiflica’ map with the previous map [14, 19] revealed high consistency of marker order with only minor exceptions. Such minor inconsistencies might be caused by missing data or chromosome rearrangements [59].

Comparative mapping

With the availability of sequenced legume genomes, synteny studies can be used to help identify candidate genes underlying specific traits in the less characterized legume species such as pea. Using the mapped pea SNP sequences identified in this study, comparative maps were generated using the whole genome sequences of four legumes from the closely related M. truncatula to the distant L. japonicus. Common bean and L. japonicus belong to the Phaseoleae and Loteae tribes, which are more distantly related to pea, while M. truncatula and chickpea belong to the Trifolieae and Cicereae tribes, which are more closely related to pea according to phylogenetic affinities [60]. Consistent with previous comparative mapping studies in pea [15, 19], this study also found pea (~4.3 Gbp [2]) to be genetically closest to M. truncatula (~470 Mbp [61]) and chickpea (~740 Mbp [62]). We found that a high level of conservation was observed with comparisons to MtChr 1, MtChr 5, and MtChr 7. In contrast, a low level of conservation was shown in MtChr 2 and MtChr 6, also consistent with previous studies [15, 60]. For common bean and L. japonicus, there were fewer large syntenic blocks as compared with M. truncatula and chickpea.

Genetic variation and correlation between traits

Having significant genetic variation in a RIL population is critical in QTL identification studies and cultivar improvement through traditional breeding practices. In this study, while the parents of the RIL population did not exhibit significant differences in some of the mineral nutrients (Ca, Mg, Mn, Mo, P and Zn concentrations; all mineral nutrient contents) large variation was observed for all the mineral nutrients in the RILs, ranging from 1.6-fold to 21-fold. Furthermore, transgressive segregation was observed for most of the traits except Fe, Mo and S concentration and all the traits revealed continuous distributions, which indicated the traits may be controlled by multiple genes, also indicating the potential for mineral nutrient improvement through breeding. Positive and negative correlations are able to be observed in all the traits in a population because of factors such as linkage or environmental differences [63]. It is interesting that positive correlations were observed between all the nutrient concentrations and contents, while the correlations between seed weight and nutrient concentrations were all negative. Similar results were found in the study of L. japonicus [36] and common bean [64], where negative correlations were observed between seed weight and Ca, Mg and Mn concentrations. The high positive correlations found between Ca and Mn concentration (0.69), Mg and Mn concentration (0.69), and Fe and S content (0.85), suggests that improving one of these mineral nutrients could simultaneously improve the other.

Mineral nutrients and seed weight QTLs

In this study, we conducted a comprehensive QTL analysis in a RIL population for pea mineral nutrients (B, Ca, Fe, K, Mg, Mn, Mo, P, S and Zn) to identify genetic loci associated with these traits. A total of 89 QTLs were identified for twenty-one different traits: 46 QTLs for seed mineral concentrations, 37 QTLs for seed mineral contents, and 6 QTLs (SW-Ps3.1, SW-Ps5.1, SW-Ps5.2, SW-Ps6.1, SW-Ps6.2, SW-Ps7.1,) for 100-seed weight. The previous studies on other legumes found 46 QTLs for eight seed mineral concentrations (Ca, Cu, Fe, K, Mg, Mn, P and Zn), 26 QTLs for seven seed mineral contents (Ca, Cu, K, Mg, Mn, P and Zn) and 3 QTLs for seed weight for M. truncatula [30]; 34 QTLs for ten seed mineral concentrations (Ca, Cu, Fe, K, Mg, Mn, Ni, P, S and Zn) and 48 QTLs for nine seed mineral contents (Ca, Cu, Fe, K, Mg, Ni, P, S and Zn) for L. japonicus [36]; 79 QTLs for ten seed mineral concentrations (B, Ca, Cu, K, Mg, Mn, Na, P, S and Zn) and 9 QTLs for iron and zinc contents for common bean [26, 31–35]. Of the three previously reported association studies for mineral nutrients in pea [37–39] only one [38], which looked at Fe, Zn and Se mineral nutrients, had common markers that could be used to compare co-localization of the corresponding mineral nutrient QTLs identified in this study. On the basis of the common markers between the linkage map from Diapari et al. [38] and the one from the current study, the corresponding positions of markers on these two linkage maps were able to be estimated. Based on common marker Agps1, [Zn]-Ps3.1 appeared to co-localize with PsC7872p386, a marker associated with Zn concentration on LG III. Three Fe concentration QTLs ([Fe]-Ps5.1, [Fe]-Ps5.2 and [Fe]-Ps7.1) appeared to be at the same end of the LG V and LG VII with PsC5316p234 and PsC12961p224. However, the positions of the Fe concentration QTLs in the other two studies [37, 39] could not be precisely estimated because of a lack of common markers. For seed weight, five QTL studies [65–69] have been reported in pea and some QTLs were located in the same linkage groups (LG III, LG V, LG VI, and LG VII) with the QTLs found in the current study. However, it was difficult to estimate if these QTLs co-localized with each other, because of insufficient common markers. Adding common markers between previous and present linkage maps will be needed to facilitate further comparisons of the mineral nutrients and seed weight QTLs among studies.

QTL co-localization and candidate gene prediction

Co-localizations of the mineral concentration and content QTLs were observed throughout the seven linkage groups in this QTL analysis. Most interesting was the co-localization of B, Ca, Fe, K, Mg, Mn, Mo, P, S, Zn concentration QTLs with Ca, Fe, Mo, P, Zn content QTLs and the seed weight QTL in LG V (29.6–52.4 cM). This finding was not unexpected as previous mineral nutrient QTL studies in other legumes have shown similar co-localizations on L. japonicus chromosome 1 [36], P. vulgaris chromosomes 1 and 6 [31], and M. truncatula chromosomes 1 and 7 [30]. During seed development, seed mineral accumulation relies on continued uptake, translocation, and remobilization of mineral nutrients from different vegetative and reproductive tissues to sink tissues (seeds), with most of these minerals ultimately entering seeds via the phloem pathway [70]. The LG V co-localization region discovered in this study might point to a locus that governs the whole-plant mobilization of several of these minerals (e.g., some global aspect of phloem translocation). Alternatively, for some of the metals, this locus may condition for a metal chelator [71] or a common transporter (e.g., the ZIP transporter protein family members can transport Fe, Mn, and Zn [72] or the CAX transporter family can be bi-functional for Ca and Mn [73]).

To estimate whether co-localizations found in other legumes were syntenic with the co-localization discovered in this study, the comparative maps between those species were developed by mapping the pea SNP sequences to the genome sequences of these legumes. The co-localization of Ca, Cu, Fe, K, Mg, Mn, P, Zn concentration QTLs observed on M. truncatula chromosome 7 by Sankaran et al. [30] was found to be syntenic with the LG V mineral co-localization region found in the current study. Putative mineral-related genes located in the vicinity of the co-localization on M. truncatula chromosome 7 have been identified, which included Ca²⁺ cation antiporter (CAX7), heavy metal transporting P-type ATPases (HMA), and the cation/H⁺ exchanger (CHX). Also, the co-localization of Fe, Zn concentration QTLs and seed weight QTLs observed on P. vulgaris chromosome 1 by Cichy et al. [31] was collinear with the LG V mineral co-localization region in pea. Unfortunately, no putative mineral-related genes were reported in the study of Cichy et al. [31].

In addition, a QTL associated with seed surface (wrinkled/round) with an explained variance of 78.1% was co-localized with the LG V mineral and seed weight co-localization region (unpublished data). This observation was also validated through the comparative mapping between pea and Medicago truncatula: a gene coding starch branching enzyme I which controls wrinkled/round seed trait [74] was detected in the LG V mineral co-localization region. Furthermore, a SNP marker Tip_SNP2, located in the LG V mineral co-localization region has significant marker-trait association for the seed surface trait [39]. According to the documented characteristic of pea core collection through the Germplasm Resources Information Network, there is significant difference in terms of seed mineral concentration between round-seed germplasm and wrinkled-seed germplasm: mineral concentrations of wrinkled-seed germplasm are higher than round-seed germplasm. Also, it was reported that there was significant difference in terms of seed weight between round and wrinkled pea [75, 76]: pea seeds with round surface have higher seed weight than pea seeds with wrinkled surface. Although the co-localization was found here in terms of seed surface, seed weight and seed mineral nutrients, further study is needed to confirm a mechanistic association between these traits.

Choice of QTLs for MAS to increase mineral nutrients

This comprehensive mineral nutrient study provides a foundation of QTLs for use in discovery of markers to improve mineral concentration and content in pea cultivar development. The QTLs, [B]-Ps5.1, [Ca]-Ps5.1, [K]-Ps5.1, [Mg]-Ps5.1, and [Mo]-Ps5.1, identified in this study explained 42%, 31%, 43%, 43.3%, and 34.2% of phenotypic variations respectively, with ‘Kiflica’ as the parent contributing the high-concentration alleles. Additionally, these QTLs were stable over two environments and were located in a similar region of the genome. This locus is possibly a major-effect region that can be used for improvement of multi-mineral levels together. However, further effort will be necessary before use in breeding. Saturating the linkage map through the addition of more DNA markers will increase the coverage of the whole pea genome and reduce confidence intervals of the QTLs. Furthermore, since QTLs have been reported to be lost or not effective in different population backgrounds [77, 78], it would be important to validate the targeted QTLs in different mapping populations. The future release of a pea reference genome sequence [79] will enable a more complete analysis of these trait loci.