Genome-wide SNP identification, linkage map construction and QTL mapping for seed mineral concentrations and contents in pea (Pisum sativum L.)

Background Marker-assisted breeding is now routinely used in major crops to facilitate more efficient cultivar improvement. This has been significantly enabled by the use of next-generation sequencing technology to identify loci and markers associated with traits of interest. While rich in a range of nutritional components, such as protein, mineral nutrients, carbohydrates and several vitamins, pea (Pisum sativum L.), one of the oldest domesticated crops in the world, remains behind many other crops in the availability of genomic and genetic resources. To further improve mineral nutrient levels in pea seeds requires the development of genome-wide tools. The objectives of this research were to develop these tools by: identifying genome-wide single nucleotide polymorphisms (SNPs) using genotyping by sequencing (GBS); constructing a high-density linkage map and comparative maps with other legumes, and identifying quantitative trait loci (QTL) for levels of boron, calcium, iron, potassium, magnesium, manganese, molybdenum, phosphorous, sulfur, and zinc in the seed, as well as for seed weight. Results In this study, 1609 high quality SNPs were found to be polymorphic between ‘Kiflica’ and ‘Aragorn’, two parents of an F6-derived recombinant inbred line (RIL) population. Mapping 1683 markers including 75 previously published markers and 1608 SNPs developed from the present study generated a linkage map of size 1310.1 cM. Comparative mapping with other legumes demonstrated that the highest level of synteny was observed between pea and the genome of Medicago truncatula. QTL analysis of the RIL population across two locations revealed at least one QTL for each of the mineral nutrient traits. In total, 46 seed mineral concentration QTLs, 37 seed mineral content QTLs, and 6 seed weight QTLs were discovered. The QTLs explained from 2.4% to 43.3% of the phenotypic variance. Conclusion The genome-wide SNPs and the genetic linkage map developed in this study permitted QTL identification for pea seed mineral nutrients that will serve as important resources to enable marker-assisted selection (MAS) for nutritional quality traits in pea breeding programs. Electronic supplementary material The online version of this article (doi:10.1186/s12870-016-0956-4) contains supplementary material, which is available to authorized users.


Background
Pea (Pisum sativum L.), an important pulse crop, is widely grown for human and animal consumption. It is the plant used by Gregor Mendel to illustrate the principle of genetics [1] and has long been considered a good source for protein, carbohydrates, minerals and vitamins. Associated with high nitrogen fixation, pea plays a vital role in the crop rotation system. In recent years, pea yield production worldwide has exceeded ten million tons. In 2014, the major producers of dry peas were Canada (30.4%), China (13.9%), Russia (13.3%), the United States (6.9%) and India (5.3%) (FAOSTAT, 2014).
With an estimated genome size of~4.3 Gbp [2] and a high repeat component estimated to be between 50% and 60% [3,4], the improvement of genetic and genomic resources for pea is required for marker-assisted breeding (MAB). MAB, which involves the use of DNA markers to predict trait performance is widely used in crop breeding [5]. Identification of genome-wide markers has undergone a revolutionary transition over the last few years with the advent of low-cost and high-throughput genotyping by sequencing technology [6]. In comparison to traditional marker discovery, GBS can be combined with marker genotyping, allowing marker discovery and genotyping to be completed at the same time. This assay was developed by Elshire et al. [7] and has been used as a tool in linkage mapping, QTL discovery, genomics-assisted breeding, and genomic diversity analysis in a large range of crops, including barley and wheat [8], rice [9], sorghum [10] and switchgrass [11]. While more than fifty-two genetic linkage maps are available for pea [12], eight are highdensity SNP-based [13][14][15][16][17][18][19][20] with only one [20] developed using GBS.
Mineral nutrients are inorganic elements essential for plant and animal growth and development [21]. Based on their quantitative requirements, plant mineral nutrients are classified into two groups, macroelements and microelements. Macroelements, generally found in plant tissues in the mg/g dry weight range, include nitrogen, phosphorus, potassium, calcium, magnesium and sulfur. Microelements include boron, copper, iron, chloride, manganese, molybdenum, and zinc, and are found in plants at the μg/g dry weight or lower range. For humans, plant foods are an important source of essential minerals, but unfortunately, mineral deficiencies are a major concern in global health [22] with over two-thirds of the world's population estimated to experience inadequate intake of one or more mineral nutrients, with more than half considered iron deficient and over 30% zinc deficient [23]. Nutritional deficiencies are especially prevalent in developing countries where people do not have the resources to adequately diversify their diets with vegetables, fruits and animal products. These mineral nutritional deficiencies can lead to stunted growth and development in children, lower resistance to disease, and increased mortality rates [24]. Improving the levels of minerals in foods, through the process of biofortification, has been proposed as a strategy to help combat these dietary deficiencies. Biofortification through traditional plant breeding or biotechnology can be a powerful and sustainable approach to significantly increase nutrient concentrations in crops [25]. Food legumes provide essential nutrients and usually contain higher concentrations of mineral nutrients than do cereals and root crops [26]. Pea is one of the crops targeted for biofortification and has long been recognized as a valuable, nutritious food for the human diet. According to a study conducted with six different cultivars across seven locations by Amarakoon et al. [27], a single serving of cooked pea seeds (100 g fresh weight) can supply 58-68% of the recommended daily allowance (RDA) of iron for male aged from 18 to 50 years, and 26-30% of the RDA of iron for female in the same age group; 36-58% of the RDA of zinc for male, 48-78% of the RDA of zinc for female. The mineral variation within pea germplasm provides the potential to create new pea cultivars with greater mineral density.
To begin to improve levels of mineral nutrients in pea seeds, an understanding of the genetic basis of these traits is required. The accumulation of mineral nutrients in seeds is determined by a series of complex processes that begin with uptake from the rhizosphere, membrane transport in the roots, translocation and redistribution within the plants through the xylem and phloem systems, and import and deposition in the seeds [28]. To date, genes associated with translocation of several elements have been identified in Arabidopsis thaliana, but only limited research has been done in pea [23,29]. Identification of QTLs provides a valuable platform to help identify the genetic basis underlying phenotypic traits. Previous studies on QTLs for mineral nutrients in legumes have been reported on the model legume Medicago truncatula [30], common bean [26,[31][32][33][34][35] and Lotus japonicus [36]. However, so far, there are only three QTL studies dealing with mineral nutrients in pea, all of which used association studies in diverse populations [37][38][39]. Kwon et al. [37] discovered ten DNA markers for seven mineral nutrients (Ca, Cu, Fe, K, Mo, Ni and P), while Diapari et al. [38] discovered nine SNP markers associated with iron and two related with zinc in seeds. In addition, Cheng et al. [39] found five SNP markers associated with calcium and magnesium.
Comparative genetics is used to identify syntenic regions controlling traits of interest among closely related species [40]. Within the legumes, the sequenced genomes of Medicago truncatula, Cicer arietinum, Phaseolus vulgaris, and Lotus japonicus can be used to transfer knowledge such as trait loci and underlying genes to less studied crops like pea.
The focus of this study was to develop a series of genomic tools to enable mineral improvement in pea through marker-assisted cultivar development. Increasing seed mineral concentration can be influenced by several factors, including seed weight, slow plant growth and low seed yield [36,41]. Additionally, a previous QTL study of mineral nutrients in Lotus japonicus identified several seed mineral concentration QTLs co-localized with QTLs associated with average seed mass. This suggested that higher seed mineral concentrations might be inversely correlated with seed weight [36]. Therefore, to avoid the utilization of loci associated with high seed mineral concentration but low seed weight, this study also assessed QTL for 100seed weight and seed mineral content. The objectives of this study were to (1) develop genome-wide SNPs using a GBS approach, (2) construct a high-density genetic map using a RIL population, (3) establish comparative maps between pea and the closely related legumes, and (4) identify QTLs associated with seed weight and mineral concentration and content.

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 6 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 SSRbased map of Loridon et al. [43]. PCR amplifications were performed with 4 ng genomic DNA, 1 × PCR buffer, 1.5 mM MgCl 2 , 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 2 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 Spectro-Desinger 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 2 , 0.5 mM dNTPs, 0.1 μM each PCR primer, 1 U Taq DNA polymerase (Bioline), and 1.8 μl ddH 2 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].

Phenotyping
The field experiments were established in two locations, Whitlow (46°74'N 117°13'W) and Spillman (46°43'N 117°1 0'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 2 ) was calculated for each trait as where σ G 2 = genotypic variance, σ GE 2 = variance due to interaction between genotype and environment, σ e 2 = 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
Two GBS libraries were constructed using 96 barcodes and the ApeKI restriction enzyme to generate SNP data on the two parents and the RILs. A total of 384.7 million reads were obtained from the high-throughput sequencing and 349.6 million reads (91% of total raw reads) met the UNEAK pipeline's quality control. The number of reads per sample ranged from 0.8 million to 5 million reads with an average of two million reads. In the analysis, identical reads were defined as a tag. A total of 45 million tags were identified in the entire set of reads with the number of tags per sample ranging from 130 K to 481 K, for an average of 255 K (Table 1). 1.2 million tags corresponding to 336.2 million reads (87% of total raw reads) met the minimum standard of 5 reads per tag and were used for SNP calling. Following pairwise alignment, 104 K tag pairs were identified. A total of 3,095 SNPs with a MAF > 0.4 and a sample calling rate > 0.5 were called by the UNEAK pipeline. In order to ensure the SNPs were high quality, only SNPs with homozygous genotypes in both parents and with less than 20% missing data per sample were kept for further analysis. High-quality SNPs (1609) were identified and used for linkage map construction. The sequence reads have been uploaded to the NCBI SRA database with accession number SRP092012.

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

Phenotypic analysis
The mean values of mineral nutrient concentration, mineral nutrient content and 100-seed weight for the two parents and the RILs across the two locations are listed in Table 3. Also, the table shows the coefficient of variation and ranges of nutrient concentration, nutrient content and seed weight for the RILs. 'Kiflica' had higher nutrient concentration and content than ' Aragorn' , while 'Kiflica' had lower seed weight than ' Aragorn'. Seed mineral concentration and content ranged from 1.6-fold to 21-fold across the RILs and seed weight varied 2-fold. All the seed traits showed high degrees of correlation between the RILs grown in both locations ( Table 3). The P concentration showed the lowest value of correlation (0.27), while Ca concentration showed the highest correlation (0.91). From the frequency distribution histograms shown in Additional file 3, all the traits revealed continuous distribution in two locations and transgressive segregation except for the Fe, Mo and S concentrations. The ANOVA table shown in the Additional file 4 indicates that all the genotypes had significant differences in all the traits. In terms of the environmental effects, there were significant differences (P < 0.05) in all the traits with the exceptions of Mn concentration and content, B content and 100-seed weight. Genotype by environment interactions had no significant effect (P < 0.05) in most of the mineral concentration traits

QTL analysis
QTL analysis was performed for all the seed traits in the ' Aragorn' x 'Kilfica' RILs across the two locations. A total of 46 QTLs were identified for seed mineral concentrations, 37 QTLs for seed mineral contents, and 6 QTLs for 100-seed weight (Tables 4 and 5). The QTLs were named following the convention of Hamon et al. [56]. The QTLs explained from 2.4 to 43.3% of the phenotypic variance. Co-localizations of QTLs on the seed traits were detected in both locations.
Five QTLs were identified for B concentration, two of which were detected in both locations. It is worthwhile to note that the QTL, Five Fe concentration QTLs were identified and three of them were observed in both locations. One 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 colocalization 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 methylationsensitive 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][9][10][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  [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  [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][32][33][34][35]. Of the three previously reported association studies for mineral nutrients in pea [37][38][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.  [37,39] could not be precisely estimated because of a lack of common markers. For seed weight, five QTL studies [65][66][67][68][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 colocalizations 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 bifunctional 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 colocalization region found in the current study. Putative mineral-related genes located in the vicinity of the colocalization on M. truncatula chromosome 7 have been identified, which included Ca 2+ 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 colocalization 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.  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.

Conclusion
In this study we developed 1609 high-quality SNPs though GBS, constructed a linkage map and identified 89 QTLs for mineral concentration, content and seed weight in pea. Comparative mapping between pea and four sequenced legumes showed regions of synteny that can be utilized for identifying candidate genes and beneficial alleles to improve mineral nutrient levels in pea cultivar development. The QTLs discovered in this study have potential use in MAS through trait predictions in breeding germplasm. As seed nutrient traits are physiologically complex, further