Characterization of Genetic Diversity and Genomic Prediction of Secondary Metabolites in Pea Genetic Resources

1. Introduction

Field pea (Pisum sativum L.) is a globally important cool season grain legume that is valued for its high yielding ability [1,2], a moderately high seed protein content reportedly in the range of 23-30% [3,4], and various agronomic benefits, such as nitrogen fixation, that enhance the sustainability of cropping systems [5,6]. An increasing interest for pea production in the European union has been driven by the demand for high-protein feedstuff to replace imported plant protein [7], and by the surge in demand for plant-based foods, driven by health, environmental, and ethical concerns [8,9,10].

Domesticated approximately 10,000 years ago in the Near East, pea gradually expanded its cultivation across the Mediterranean basin, China, Northern and Eastern Africa [11]. The genetic diversity of modern cultivars is very narrow [11]. Landraces and old cultivars, although largely available in germplasm collections, remain largely underutilized in plant breeding because of some unfavorable agronomic characteristics (e.g., the leafy trait) and the possibly high content of seed antinutritional factors (ANFs) [12,13]. Information on the diversity for seed quality traits of large pea germplasm collections focused largely on primary metabolites, especially protein content e.g. [14,15], whereas limited information is available for micronutrients and secondary metabolites. Among the various ANFs found in pea seeds, saponins, raffinose-family oligosaccharides (RFOs), and phenolic compounds are some of the most problematic, with negative effects on digestibility in monogastrics. However, some of these molecules have also been reviewed for their positive health-related effects for human diets [16].

Saponins are a diverse group of secondary metabolites that are widely spread in plant species [17]. Pea seeds contain triterpene saponins such as soyasaponin of the B class (of which Soyasaponin I is found in pea) and DDMP saponin (soyasaponin βg), which contribute to bitterness and astringency [18,19,20,21]. Saponins in pea and other grain legumes were also reported for their damaging effect on cell membranes [22]. Saponin content in different pea genotypes was found to be related to biotic stress response [23], confirming the role of these compounds in plant defense mechanisms. Variety differences in saponin content have been reported in pea, yet their variation across broader germplasm remains understudied. Several steps of DDMP saponin biosynthesis have been studied in different grain legumes, including pea [24]. [25] identified a major gene, BAS1, as being responsible for most of the DDMP saponin content in pea. In a recent work, [26] successfully developed a CRISPR/Cas9 protocol for the development of mutant lines with reduced DDMP saponin content in the seed.

The raffinose family of oligosaccharides (RFOs) includes α-D-galactosides of sucrose that are widely accumulated in pea seeds. The digestive tract of humans and monogastric animals lacks the α-D galactosidase, an enzyme that hydrolyzes the α-D (1 -> 6)-glycosidic bonds between molecules of galactose moieties present in the oligosaccharides. Therefore, RFOs are digested with the involvement of the bacterial microflora homing further sections of the intestine. Excessive amounts of carbon dioxide and hydrogen are produced as a result of sugar decomposition and fermentation of the released monosaccharides, causing flatulence and discomfort [27]. Therefore, RFOs are widely regarded as antinutritional factors causing flatulence in humans and monogastric animals [28]. In extreme cases, the increased content of RFOs in feed ratios for monogastric livestock can provoke diarrhea and reduce energy use.

Pea seeds are rich in polyphenols, a group of secondary metabolites primarily localized in the seed coat that has potential health benefits for humans. Polyphenols are known for their antioxidant activity and are involved in plant defense mechanisms [29]. Despite the positive health effects associated with polyphenols consumption, these secondary metabolites are also known as ANFs. In particular, condensed tannins hinder the digestibility of pea in monogastric animals [12] and decrease the iron absorption through chelation [30]. The core of the polyphenol biosynthetic pathway is well conserved across plant species [31]. The phenolic profiling of pea and other pulses was recently performed through liquid-chromatography and mass spectroscopy (LC-MS) [30,32]. The variation in pea seed coat polyphenols is associated with flower color, with purple color genotypes having specific classes of polyphenols, known as tannins, that are lacking in white flower genotypes, thus resulting in higher total polyphenols content in the former group [33]. The genetic regulation of flower color has been widely studied, notably as a character used by Mendel in the study of inheritance in pea, with two major genes being identified in the pea genome as responsible for the white color mutation [34].

The rapid development of next-generation sequencing technologies has revolutionized the study of genetic diversity, with genotyping-by-sequencing (GBS) emerging as a cost-effective and efficient approach for large-scale SNP discovery and genome-wide association studies (GWAS). Quantitative trait loci (QTL) and candidate genes were identified for pea seed quality traits relative to concentrations of crude proteins, amino acids, minerals and fibers, as well as protein digestibility [15,35,36,37,38,39,40,41]. Genomic prediction models, which integrate phenotypic and genotypic data to estimate the breeding values of untested genotypes [42], have also been envisaged for selection of traits, such as protein content, that exhibited a polygenic control [38,39,43], given the higher efficiency of genomic selection relative to marker-assisted selection for such traits [44].

While genome-wide association analysis (GWAS) and genomics have been applied extensively for the study of agronomic traits in pea, secondary metabolite traits remain underexplored within those frameworks. To address these gaps, this study was designed to (i) characterize for seed saponin, RFOs, and total phenolic compounds (TPC) and antioxidant activity (AA) of 156 pea accessions belonging to a collection previously described by [14], (ii) conduct GWAS to unveil the genetic architecture of these traits and identify associated genomic loci and potential candidate genes, and (iii) assess the scope and predictive ability of genomic prediction models for these traits.

2. Results

2.1. Phenotypic Trait Variation

The 156 pea accessions were grouped into 20 germplasm pools that included 19 regional landrace/old cultivar pools and one modern cultivar pool (Supplementary Table S1). The mean and range values of the germplasm pools for the seed content of eight secondary metabolites is shown in Table 1. The analysis of variance (ANOVA) showed significant variation among germplasm pools for the content of all traits (Table 1). The results of Tukey HSD mean comparisons are provided in Supplementary Table S2. On average, the highest and lowest contents of TPC were observed in accessions from West Asia and East Balkans, respectively. High TPC values also occurred in accessions from Greece and Afghanistan, whereas low values were also a feature of modern cultivars and accessions from the UK and Nepal. Highest and lowest values in vitro antioxidant activity (AA) featured the accessions from West Asia and Germany, respectively. On average, the soyasaponin βg saponin (hereafter referred to as Ssβg) displayed about 20-fold higher content than the soyasaponin I saponin (hereafter referred to as Ss1). The accessions from Germany and the UK had particularly high values of the former saponin, while those from France exhibited the lowest values. Materials from the UK had by far the highest content of the Ss1 saponin, while those from Georgia and France had the lowest values. Within RFOs, verbascose showed the highest mean value in the data set (8.52 mg/g), followed by stachyose (6.74 mg/g) and raffinose (2.47 mg/g). The mean value for sucrose was 5.90 mg/g. On average, UK accessions had the highest verbascose and stachyose contents, while those from Georgia and Ethiopia had the lowest content of these compounds. The accessions from China and Ukraine showed the highest and lowest content, respectively, of sucrose. Finally, the accessions from West Asia and North Africa had the highest and lowest content, respectively, of raffinose. Besides showing low TPC values, the germplasm pool including the modern cultivars was characterized by fairly low AA, intermediate content of saponin compounds, and fairly low levels of all RFOs (Table 1).

Several landrace/old cultivar germplasm pools (Turkey, Spain, France, UK, Italy, Greece) exhibited large within-pool variation across traits as indicated by the mean value of genetic coefficient of variation (CVg) values; in contrast, generally low within-pool variation was observed for a few landrace/old cultivar germplasm pools (Germany, Russia), as well as for the modern cultivar pool (Supplementary Table S3). On average, the within-pool genetic variation was particularly high for AA and the two saponins (Supplementary Table S3).

The association among traits was investigated through correlation analysis and principal components analysis (PCA). The highest correlations were found between metabolites belonging to the same category. In particular, high correlations (≥ 0.63; P < 0.01) were observed between the two saponin traits, between TPC and AA, and between the RFO compounds stachyose and verbascose or stachyose and raffinose (Table 2).

The first 2 axes of the PCA accounted for about 73% of the overall trait variation (Figure 1). The first axis tended to represent a contrast between germplasm pools with high saponin and verbascose contents (featuring with high score values) versus pools with high phenolic content and high antioxidant activity (featuring low score values) (Figure 1). The second PCA axis tended to represent a contrast between material with high values of the three soluble sugars (sucrose, raffinose, and stachyose) and, to a lesser extent, high phenolic content and antioxidant activity, featuring high score value and represented especially by the germplasm pool from West Asia, versus material with the opposite characteristics (featuring low score values; Figure 1). The PCA plot indicated the similarity of a few landrace/old cultivar germplasm pools with close geographical origin, such as those from Central Asia and India, or from Russia and Ukraine, or from Italy and Spain. On the whole, however, the ordination of the germplasm pools in the space of the first two PCA axes was not closely related to their region of origin. The accessions from UK and, to a lesser extent, those from Germany featured high score values on the first axis, whereas those from Afghanistan, Georgia, Turkey and France had low score values on this axis (Figure 1). The accessions from West Asia were strongly distinct along the second PCA axis.

To further explore the similarity of the germplasm pools based on their metabolic profile, we conducted a hierarchical cluster analysis that resulted in the optimal grouping of the pools into seven groups, partly reflecting the geographic similarity of the pools. The heatmap reported in Figure 2 simultaneously visualizes the seven groups and their association with their variation in metabolite contents. One group included German and UK accessions, featuring high saponin and RFO contents and low TPC and AA. A second, contrasting group included only the material from West Asia, characterized by high TPC and AA and relatively low saponin content. Two other groups (indicated by orange and light blue colours in Figure 2) tended to include germplasm pools of various origins having high TPC and AA and low RFO and saponins contents. A fifth group included only the germplasm from China, featuring very high sucrose content. A sixth group included the accessions from Spain and Italy, of which the main characteristic was a higher-than-average RFO content. One last group included the modern germplasm along with landrace/old cultivar groups from various regions of Eastern Europe (Ukraine, East Balkans), Russia, Central Asia, India, and Nepal (Figure 2). Most of this material featured average or below-average content of all secondary metabolites.

2.2. Analysis of Population Structure

Genomic data were available for only 151 of the 156 accessions characterized with phenotypic analyses. The structure analysis, which was based on 10,249 SNPs of the 151 pea accessions, indicated an optimal number of 11 clusters (Supplementary Figure S1; Figure 3). On the whole, the classification of the individual accessions partly reflected the origin and the type (old cultivar/landrace or modern cultivar) of the accessions. The cluster with the highest number of accessions (light green in Figure 3) comprised most of the accessions from Ethiopia and several from India. The second most numerous cluster (dark green in Figure 3) included several accessions from France and East Balkans and other single accessions from other pools. Most of the accessions from Italy, several from Spain and two from West Asia belonged to the same cluster (pink in Figure 3). Most of the accessions from Ukraine, two from Russia, and two from East Balkans belonged to the same cluster (gray in Figure 3). Also, most of the accessions from China and Afghanistan, as well as two from Central Asia, were classified in the same cluster (light blue in Figure 3). Three accessions from India and two from Nepal were assigned to the same cluster (brown in Figure 3). Finally, most of the modern cultivars were grouped in the same cluster (yellow in Figure 3), which also included an accession from Spain, one from the UK, and one from Germany. A total of 62 accessions were classified as admixed.

2.3. Linkage Disequilibrium Decay and Genome-Wide Association Study

The linkage disequilibrium (LD) decay was very fast in all chromosomes (Supplementary Figure S2). The distance at which the LD was likely due to physical linkage averaged 4669 bp with a range from 3390 to 5846 bp.

A GWAS was carried out for each traits according to two models, namely, BLINK and FarmCPU, selecting the significant SNPs according to the Bonferroni threshold at P = 0.05. Figure 4 reports the Manhattan plot results for BLINK, while those for FarmCPU are reported in Supplementary Figure S3. At least one significant SNP was identified for any trait. Three SNPs on chromosome 7, one on chromosome 6, and one of chromosome 5, were identified for TPC by the BLINK model. One SNP on chromosome 5 was found for antioxidant activity by both BLINK and FarmCPU. For this trait, BLINK identified two additional SNPs on chromosome 6, whereas FarmCPU identified four additional SNPs, of which two on chromosome 4, one on chromosome 3, and one in chromosome 7. For Ssβg, BLINK identified one SNP on chromosome 7 and a second one on chromosome 5 just below the threshold of significance, whereas FarmCPU identified one SNP in each of the chromosomes 1, 3, 5, 6 and 7. For Ss1, BLINK and FarmCPU identified the same two SNPs on chromosomes 2 and 5, while FarmCPU identified two additional SNPs in chromosomes 1 and 2. For verbascose, FarmCPU identified six SNPs (one in chromosomes 1, and three and two in chromosomes 5 and 7, respectively), while BLINK identified one SNP on chromosome 3. Interestingly, the same SNP was also significant for stachyose according to both BLINK and FARMCPU. For the same trait, FarmCPU identified one additional SNP on chromosome 3 (found just below the threshold of significance from BLINK), one in chromosome 1 and two in chromosome 7. For raffinose, BLINK identified two SNPs on chromosomes 4 and 5. Finally, for sucrose, one SNP was identified by the BLINK model on chromosome 4. Out of a total of 37 different significant SNPs identified by the GWAS, 35 were located within annotated genes, one was within the threshold of LD decay performed in our analysis, and one was not in the vicinity of any annotated gene in the pea genome. The list of the significant SNPs and their putative associated candidate genes are provided in Supplementary Table S4.

2.4. Genome-Enabled Prediction

The predictive ability of two statistical models (rrBLUP; Bayesian Lasso) envisaged for genome-enabled prediction is reported in Figure 5 for each of the eight traits. Bayesian Lasso exhibited greater predictive than rrBLUP for five traits out of eight, but the difference between models was always modest. Considering the best-predicting model, we observed high predictive ability (> 0.6) for verbascose content, and moderately high predictive ability (> 0.4) for all other traits except raffinose amount.

3. Discussion

Our study provided an assessment of the extent of pea genetic variation for eight important quality related traits that is unprecedented with respect to the number and the diversity of origin of the germplasm under study. Its results indicated a significant and considerably large variation for all traits, which was partially associated with the geographic origin of the accessions.

The range of TPC was between 0.35 and 1.07 mg GAE/g (Table 1), aligned with the values reported by [45] for 100 pea accessions (0.49-1.28 mg GAE/g) but lower than the range reported by [46] for ten pea varieties from China (0.66-2.66 mg GAE/g). The AA results, ranging from 0.16 to 2.40 µmol TE/g, were lower than those reported by [47], who found 0.19–3.97 µmol TE/g in 81 pea genotypes. Furthermore, the results were significantly inferior to the range of 6.06–12.49 µmol TE/g found by [46] in 10 pea varieties. Interestingly AA was the trait with the highest average CVg across germplasm pools, suggesting wide margins for its improvement through breeding.

Phenolic compounds, which include phenolic acids and flavonoids, stand out as extensively studied functional components of pea seeds [48]. These compounds and particularly the condensed tannins were traditionally considered ANFs because of their negative effects on protein digestibility in monogastric animals [49], but they have lately been reconsidered as valuable antioxidants able to effectively inhibit free radicals and to prevent oxidative reactions at cellular level. In addition, phenolic compounds play a pivotal role in modulating the gut microbiota by promoting the growth of beneficial species and inhibiting harmful bacteria [50], and exhibit a range of biological activities, including anti-diabetic, anticarcinogenic, cardioprotective and anti-neurodegenerative effects [51,52]. Hence, the consumption of phenolics-rich pea goes beyond a merely nutritional function, offering various potential health benefits [48,50,51]. A highly positive correlation was observed between TPC and AA, as known [53]. However, a previous study including pea and other pulses [30] highlighted a different contribution to the antioxidant activity of different types of polyphenols, suggesting to target specific types rather than general TPC by breeding work. Indeed, the increase of these compounds may result in negative effects such as iron chelation and protein precipitation. The complexity of the phenolic compounds biosynthetic pathway in grain legumes was recently highlighted by a study in faba bean [54]. This study indicated how different genotypes with the same genetic background at the low tannin gene locus (zt) had very different phenolic profiles, and how this was partially affected by the environment, suggesting that several genes are regulating the secondary branches of this pathway. A qualitative characterization of phenolic compounds of the accessions would represent a useful follow-up of this study.

The current overall variation for RFOs agree substantially with previous studies by [55] and [27], which analyzed the content of sucrose and RFO in a set of Spanish breeding lines and a pea collection from the Polish gene bank, respectively. Stachyose and verbascose consistently accounted for the highest total amount of RFO, although the maximum absolute values that we found (up to 16.20 for verbascose and 10.30 mg/g for stachyose) were lower than those reported in those studies. We identified low variation in raffinose content across germplasm pools as well as among individual accessions. This result may partly depend on the generally low raffinose content, since FT-IR based models are less effective in detecting trait variation for metabolites with low concentration. The observed high correlation between stachyose and verbascose is consistent with earlier findings and the fact that these compounds are synthesized through the same biosynthetic pathway, although different enzymes are active in the biosynthesis of verbascose in different pea genotypes [28].

Using the information on seed morphology and seed color reported for our panel by [39], we verified through t tests that the variation in RFO content was significantly associated with major seed traits. The accessions with a transparent seed coat (recessive allele) showed higher RFO levels than those with a colored testa; those with a wrinkled seed (recessive allele) had higher RFO content compared with smooth-seeded lines; and those with a green cotyledon (recessive allele) accumulated more RFO than those with yellow or dark cotyledons. These findings are consistent with previous reports indicating that seed morphology and seed color are linked to variation in oligosaccharide composition in pea, with substantially lower RFO levels observed in lines carrying the dominant alleles controlling seed shape, seed coat color, and cotyledon color [27]. Similarly, [55] reported an association of seed coat with the concentration of soluble sugars.

The observed variation in saponin content was wide, both for Ssβg (often reported as DDMP saponin in pea; range of 0.04-1.0 mg/g), which represented the main compound, and Ss1. However, a wider range of total saponin concentration, namely, 0.8-2.5 mg/g, was previously reported [19]. Advances in gene editing such as CRISPR/Cas9 might help to quickly introduce the low-saponin mutation in elite pea germplasm [26], but regulation issues might be a limiting factor for practical applications of such technologies. Therefore, the identification of accessions, such as those from the French and Georgian pools, with considerably lower saponin content than that of modern cultivars is remarkable from a breeding perspective, encouraging the selection for reduced saponins through crossing of these genetic resources with elite lines.

The analysis of population structure based on SNP data indicated a partial association between genetic diversity and geographical origin of the material. Our results agree with those in [39] based largely on the same plant material. Modern cultivars tended to be genetically distinct from landrace germplasm pools in the structure analysis. Interestingly, the metabolic profile of the set of modern cultivars was characterized by an intermediate content of saponins and RFOs and a fairly low content of TPC and DPPH, indicating scope for an improvement of the nutritional profile of modern germplasm by introgression of useful variation from old varieties and landraces.

The GWAS identified various SNPs associated with each of the focus traits, a result that supported their polygenic control. The investigation of gene ontology (https://plants.ensembl.org/) suggested that most of the candidate genes containing or in LD with the significant SNPs seem to play a role in cell metabolism, general stress response or DNA transcription regulation. Among the genes that were characterized by more specific functions, Psat5g137080, associated with Ss1 saponin and identified by both the BLINK and FarmCPU models, encodes a Acetyl-CoA carboxylase protein that plays a role in the fatty acid biosynthetic process [56]. The Ss1 saponin is lipid-derived molecule to which a sugar moiety is added in a later stage of the biosynthetic pathway to give it both hydrophobic and hydrophilic properties. Psat7g192120, here associated with TPC, has an orthologue in Arabidopsis thaliana (NCRK) that is known to code for a cysteine-rich receptor-like kinase (RLK) that is involved in plant signaling pathways, including those responding to stress and pathogen attacks [57]. Phenolics compounds are well-known to be active in plant stress response and their accumulation in plant organs is often enhanced in response to both biotic and abiotic stresses. The Psat1g013760 gene was associated with verbascose and codes for a Histidyl-tRNA synthetase protein. The Arabidopsis thaliana orthologue of this gene (HRS1) encodes a transcription factor involved in integrating nitrate and phosphate signaling to regulate root development and seed germination [58]. RFOs were shown to contribute to the control of seed germination in Arabidopsis,[59] and were correlated with enhanced seed vigor in Zea mays [60]. Our inability to identify genes with a clear disruptive function in the biosynthetic pathways of the secondary metabolites, and the identification of several genes with broad metabolic roles, indirectly supported the conclusion that these traits are under polygenic control. Such a result would support the usefulness of genome-enabled prediction models for breeding line selection and/or the identification of useful genetic resources.

Although Bayesian Lasso is theoretically more suitable than rrBLUP for genomic prediction of traits that are not controlled by many genes [61], its advantage occurred only for five of the eight secondary metabolites and, anyway, in the presence of similar predictive ability of the two statistical models for all traits. A similar predictive ability of the two models was reported earlier for other pea traits which, although polygenic, might be less genetically complex than crop yield, such as protein content, onset of flowering, seed weight, frost resistance, and tolerance to rust [39,62,63,64]. The moderately high predictive ability (> 0.4) observed for all traits except raffinose content has high perspective interest for selection of breeding lines or identification of genetic resources with a desired level of secondary metabolites. On the other hand, raffinose could probably be seen as the least important of the eight secondary metabolites, as it was the RFO with lowest content in the seed. Our results, combined with earlier results for the same germplasm collection showing moderately predictive ability (0.55) also for protein content, would allow to envisage a multi-trait genome-enabled selection for a combination of primary and secondary traits. The predictive ability of the current models for reference populations other than the current global germplasm collection (and similar material) is pending verification. The prediction model for protein content constructed from data of the current germplasm collection displayed a moderate fall of predictive ability (from 0.55 to 0.28) when applied for prediction of a completely different, genetically narrow reference population of breeding lines evaluated under climate conditions quite different from those used for the germplasm collection [39].

In conclusion, the information provided by this study on the extent and geographical pattern of genetic variation, the genetic control, and the genome-enabled prediction of the eight focus traits can contribute to more targeted and efficient breeding efforts aimed to select breeding lines and identify genetic resources with a desired level of secondary metabolites, filling a research gap and meeting the growing demand for pea production with suitable traits.

4. Materials and Methods

4.1. Plant Material

The study included 156 ecotypes or old cultivars of P. sativum subsp. sativum subdivided into 19 regional pools, each represented by 4–14 entries, and 7 modern cultivars bred in France (Attika, Dove, Isard, Spirale), Spain (Cigarron, Viriato) or Germany (Santana) (Supplementary Table S1). This collection was previously described and characterized for agronomic traits by [14]. Plants were grown in 7.5L pots filled with peat soil in greenhouses at the Norwegian University of Life Sciences (NMBU) in 2023. Each accession was represented by two pots containing 4 seeds each. The pots were randomized in the greenhouse, setting the growing conditions to 16 h photoperiod, 300 μmol/ m^-2 s^-1 photosynthetic photon flux density (PPFD), and 21-16 °C day/night temperature. The plants were fertilized once during the experiment with 3 g/pot of NPK (19-4-12) fertilizer and watered regularly to maintain optimal growth. They were progressively harvested when they reached optimal maturation, pooling the seeds belonging to the same accession. Seeds were oven-dried at 60 °C for 24 h before sample preparation for the chemical analyses.

4.2. Chemical Analyses

A dried seed sample of 100 g from each accession was milled using Retsch Twister (Retsch GmbH, Germany). Flour from the same sample was used for all chemical analyses.

4.2.1. Total Phenolic Compounds and Antioxidant Activity

The extraction of TPC was carried out according to [65]. In detail, 1 g of sample was extracted with 5 mL of an 80:20 methanol/water solution (v/v). The suspensions were submitted to ultrasound (CEIA international S.A., 115/230 Vac 1 - 50/60 Hz–400 Watt, Viciomaggio, Italy) for 15 min at room temperature, then shaken for 30 min and centrifuged (Thermo Fisher Scientific, Osterode am Harz, Germany) for 10 min at 12000 × g at 4 °C. The supernatants were filtered through 0.45 μm nylon filters. The TPC were then quantified as described in [66] with some modifications. Specifically, 200 µL of the filtered extract was added to 800 µL of deionized water and 100 µL of Folin-Ciocalteu reagent and kept in the dark for 3 min, then 800 µL of 7.5% Na₂CO₃ was added, followed by incubation for 60 min, again in the dark. The spectrophotometric quantification was carried out at 720 nm, using a Cary 60 UV–Vis spectrophotometer (Agilent Technologies, Santa Clara, CA, USA) and TPC were expressed as mg gallic acid equivalents (GAE)/g of sample. The analysis was carried out in triplicate.

To assess the in vitro AA, the phenolic extracts were submitted to radical scavenging assay using 1,1-diphenyl-2-picrylhydrazyl (DPPH) radical, according to [66]. A 0.08 mM solution of DPPH in ethanol was freshly prepared. For the analysis, 50 µL of the extract was added to 950 µL of DPPH solution. After 30 min of incubation in the dark, the absorbance was read at 517 nm. The results were expressed in µmol of 6-Hydroxy-2,5,7,8-tetramethylchroman-2-carboxylic acid (Trolox) equivalents (TE)/g of sample. The analysis was carried out in triplicate.

4.2.2. Saponins

Two replicates, each consisting of 50 ± 0.2 mg of dried flour in a 2 mL vial, were prepared for each accession. In addition, an extraction blank, i.e. a vial without pea flour, was prepared for each extraction day. Blanks and sample replicates were randomized, then each was extracted in 500 µL 80/20 methanol/water (v/v, 80% MeOH containing 1 µg/mL Ginsenoside Rb1 as ISTD). The extracts were shaken for 30 min (1400 rpm, 20 °C), then centrifuged for 5 min (12 000 rpm, 20 °C, Heraeus Fresco21, Thermo Fisher Scientific). After removing 400 µL of the supernatant to a second vial, the pellet was resuspended and re-extracted as before in 500 µL 80 % MeOH (without ISTD). Then, 500 µL of the supernatant were combined with the previous 400 µL, vortexed, filtered through 0.2 µm regenerated cellulose (RC) syringe filters (Phenomenex Inc), and finally diluted 1:10 in 80 % MeOH. A QC sample was prepared by pooling aliquots from each extract. Extracts were analyzed in newly randomized order and after injection of extraction blanks on a Vanquish Horizon UHPLC coupled to an Orbitrap IQ-X iontrap orbitrap mass spectrometer (both Thermo Fisher Scientific). The chromatographic system consisted of an autosampler, pump, an Ascentis Express C18 column (150 mm × 2.1 mm, 2.7 µm, Merck) with a 5 mm guard cartridge maintained at a stable 20 °C in a column compartment, and a photodiode array detector (PDA) recording at 205 nm and 227 nm. The mobile phase was (A) water (0.2 % FA) and (B) ACN:MeOH 3:1 (v/v, 0.2% FA) with the following gradient: 0 min 15 % B; 10-14 min 90 % B; 14.1-16 min 100 % B; 16.1-24 min 15 % B. The flow was 0.25 mL/min and 1 µL of sample was injected. The heated ESI spray chamber settings were as follows: spray –3500 V, sheath gas 40 au, aux gas 15 au, sweep gas 2 au, ion transfer tube 350 °C, vaporizer 325 °C. The MS recorded in full scan mode in the time range 1.5-16 min and mass range m/z = 500-1700 Da. Saponins were first identified by their molecular ion as well as their tandem mass spectra (MS2) acquired representatively on the QC sample. Isolation and targeted fragmentation (tMS2) were performed on 5 candidate molecular ions (Supplementary Table S5) identified from PDA and MS chromatograms. Fragmentation used collision induced dissociation (CID) with 10 ms activation time at normalized collision energies of 30 %, 45 %, and 60 % and 5 s cycle time. The isolation window was m/z = 3 Da and isolation offset m/z = +1 Da. The resolution was set to R = 60 000 FWHM (at m/z = 200) for full scan and R = 30 000 FWHM for MS2 scans. Three scans were averaged to improve MS2 spectrum quality. MS2 spectra were subsequently compared to published data (Pollier et al., 2011). Individual saponins were then quantified based on their parental extracted ion chromatogram, corrected for ISTD, as soyasaponin-I equivalents based on an 11-point standard dilution series (0.005-5 µg/mL). The QC sample was analyzed after every 10th sample injection. Quantification was done in TraceFinder v.5.1 (Thermo Fisher Scientific) and corrected for sequential extraction loss.

4.2.3. HPLC Quantification of Sucrose and RFO for FT-IR Calibration

For each genotype, one sample consisting of 50 ± 2 mg of finely ground pea flour was weighed into 15 mL polypropylene centrifuge tubes. Samples were extracted in 10 mL of 50 % (v/v) ethanol containing 50 µg/mL melibiose as internal standard (prepared from a 1 mg/mL stock solution in Milli-Q water with 0.02 % sodium azide). Tubes were vortexed and subsequently shaken horizontally at 250 rpm and 50 °C for 1 h in an incubator. After extraction, samples were allowed to stand vertically at room temperature for 15 min and then centrifuged at 4000 × g for 10 min. The resulting supernatant was diluted 1:10 (v/v) with Milli-Q water, filtered through 0.22 µm PVDF syringe filters (Merck KGaA) into HPLC vials, and used for chromatographic analysis.

Soluble sugars were quantified by High-Performance Anion-Exchange Chromatography coupled with Pulsed Amperometric Detection (HPAEC–PAD) using a Dionex ICS 5000+ system (Thermo Fisher Scientific, USA) equipped with an AS-AP autosampler, ICS 5000+ SP pump, and ICS 5000+ DC column oven and detector compartment. The detection system consisted of a pulsed amperometric detector with a gold working electrode and an Ag/AgCl reference electrode, operated using the following waveform: 400 ms at 0.11 V; ramp 10 ms from 0.10 V to −2.00 V; 10 ms at −2.00 V; ramp 10 ms from −2.00 V to 0.60 V; ramp 10 ms from 0.60 V to −0.10 V; and 60 ms at −0.10 V. Separation was achieved on a CarboPac PA-1 analytical column (Thermo Fisher Scientific) with a CarboPac PA-1 guard column, maintained at 25 °C. The mobile phase consisted of A) MilliQ-water, B) 200 mM NaOH and C) 100 mM NaOH, 500 mM Na-acetate using the following A:B:C (v/v) gradient program: 0-5 min 49:49:2, 28-31 min 48:48:4, 31-34 min 0:0:100, followed by re-equilibration at start conditions. The injection volume was 25 µL, and the flow rate was 1 mL min⁻¹. Quantification was performed using external calibration with mixed sugar standards (sucrose, raffinose, stachyose, and verbascose) prepared in the same solvent system, with melibiose as internal standard.

4.2.4. FT-IR-Based Modeling for Sucrose and RFO

The content of soluble sugars, namely sucrose, raffinose, verbascose and stachyose content were estimated by Fourier Transform Infrared Spectroscopy (FT-IR) using Partial Least Squares (PLS) calibration models using the HPLC results for individual sugars as reference.

For the FTIR measurements a Bruker Invenio spectrometer (Bruker Optics, Billerica, MA, USA) equipped with a Pike MIRacle diamond crystal ATR accessory with a Deuterated Triglycine Sulfate (DTGS) detector was used for data collection. The spectra were recorded in the region between 4000 and 600 cm-1 with a spectral resolution of 4 cm-1 based on 32 scans. Before each sample a background spectrum of the empty ATR crystal was recorded to compensate for water vapour and CO₂. Three replicates were measured for each sample. All FTIR spectra were imported to the Aspen Unscrambler version 14 (AspenTech Inc., USA) and spectra were processed by applying a Savitzky-Golay 2nd derivative with a 2.degree polynomial and a window size of 9. PLS models were calculated from average spectra for the spectral range 1800-700 cm-1 using the individual sugars as reference. This resulted in the following models: for sucrose root mean square error of cross validation RMSECV=0.9 and R2=0.83 (7-factor model); for verbascose RMSECV=1.27 and R2=0.77 (6-factor model), for raffinose RMSECV=0.26 and R2=0.77 (7-factor model) and for stachyose RMSECV=1.2 and R2=0.46 (3-factor model). For further statistical analyses sugar content from three individual replicate spectra was predicted from these models

4.3. Statistical Analyses

For each trait, a linear model with the random effects accession and replicate was fitted using the lmer() function from the lme4 package [67] of the R software with the following formula:

$$Y_{ij}=\mu +u_{i }+v_j+{ϵ}_{ij }$$

where $Y_{ij}$ is the phenotype; $\mathrm{\mu }$ is the overall mean; $u_{i }$represents the random effect of the i-th genotype; $u_{i }\sim N(0,{\sigma }_u^2)$; $v_j$represents the random effect of the j-th replicate; $v_j\sim N(0,{\sigma }_v^2)$; and ${ϵ}_{ij }\sim N(0,{\sigma }_{ϵ}^2)$ is the residual error term. Best linear unbiased predictors (BLUPs) adjusted means for each accession were calculated by adding the random effect deviations to the overall intercept and used as input data for the subsequent correlation analyses, multivariate analyses, and genome-wide association studies (GWAS). For TPC, DPPH and soluble sugars, the models were fitted using data from three analytical replicates per accession, while two replicates per accession were used for saponins.

The response variables were evaluated for normality using both visual inspection (histogram and Q–Q plot of model residuals) and the Shapiro–Wilk test. Although for some traits the Shapiro–Wilk test indicated departures from perfect normality (p < 0.05), the visual assessment suggested approximately normal distributions without major deviations, supporting the suitability of the data for mixed model analysis.

An analysis of variance (ANOVA) that included the fixed factor germplasm pool and the random factors accession within germplasm pool and replicate assessed the variation between the 20 germplasm pools (hence, holding accession within germplasm pool as the error term for germplasm pool comparison).

To assess the magnitude of genetic variation within each germplasm pool, the variance between accessions of each germplasm pool was estimated and the genotypic coefficient of variation (CVg) was calculated as the square root of the genotypic variance divided by the trait mean value for the relevant pool.

Pearson’s correlation coefficients were estimated from BLUPs of individual accessions to test the association among traits. Multivariate patterns of variation among germplasm pool accessions were investigated by principal component analysis (PCA) and hierarchical clustering analysis (HCA). These analyses were performed on the BLUPs matrix of germplasm pool by trait after trait standardization to zero mean and unit variance. PCA and HCA were performed with PCA() and hclust() function of the factormineR [68] and the stats package function of R, respectively. For HCA, the distance matrix was calculated using Euclidean distance with the dist() function of the stats package according to Ward clustering algorithm. The optimal number of clusters was chosen both through visual assessment of the dendrograms and the silhouette method after generating the relative plot with the package cluster() [69]. Finally, a heatmap was created with the Heatmap() function from the ComlexHeatmap package [70], to simultaneously visualize HCA results for the metabolites and accessions. All other statistical analyses were performed by the R package stats (R Core Team, 2025).

4.4. GBS, SNP Calling, Marker Filtering and Imputation

Information on DNA isolation and GBS-based genotyping can be found in [71]. SNP calling was performed using Legpipe2 pipeline [72] with default settings for diploid species, with a preliminary filtering for mapping quality (MQ < 40). For alignment, we used the Pisum sativum L. (2n = 14) reference genome version 1a [73]. The whole set was filtered for minor allele frequency (MAF) > 5%, missing rate < 50% and imputed with the k-Nearest Neighbors (KNN) method. This process retained 10,249 SNP markers.

4.5. Linkage Disequilibrium, Population Genetic Structure and Genome Wide Association Study (GWAS)

Genomic data were available for only 151 of the 156 accessions characterized with phenotypic analyses; therefore, only 151 accessions were considered for the following genomic analyses. The linkage disequilibrium (LD) for each single chromosome was calculated using only the SNP markers with known chromosomal position relative to the reference genome [73]. LD was estimated as the squared allele frequency correlations (r^2) for each pairwise combination of markers distanced within 100 Kbp with the pcor.shrink () function of the package corpcor in R studio. LD was then plotted against the genomic distance between marker pairs and the LD decay was visualized with the LOESS regression model. The LD decay was estimated as the point where the fitted curve reached half of its maximum value [74].

The population structure was analyzed with the snmf() and Q() functions of the R package LEA [75] which output the results in a similar way as the STRUCTURE software although it employs different algorithms and was indicated as being more accurate for self-pollinating species [39]. Filtered and imputed SNP data, as described in the section above, were used as input in the analyses. The optimal number of genetic clusters was selected after visual assessment of the cross-entropy parameter plot estimated through cross-validation.

A GWAS was conducted for the eight phenotypic traits for 151 pea accessions using two models, namely FarmCPU and Blink, using the R package “GAPIT” [76]. The population structure was accounted for by including the first 2 principal components of a PCA conducted on the SNPs data (Supplementary figure 2). The appropriate compensation of the population structure was assessed through the visual examination of quantile-quantile (QQ) plots. For trait which showed an over-compensation of population structure (Ssβg and Ss1) GWAS models with reduced number of principal components (1 or 0) were performed. Significant SNPs were selected according to the Bonferroni threshold at 5%. Candidate genes associated with the significant SNPs and their putative function were searched with the pea genome browser (https://www.pulsedb.org/), investigating chromosomal regions corresponding to the chromosome-specific LD extent flanking the significant SNP.

4.6. Genomic Prediction

The SNP markers retained after SNP calling and filtering were employed to build GS models for the eight phenotypic traits. We compared two regression models, namely ridge-regression BLUP (rrBLUP) and Bayesian Lasso. The former model [77] assumes a linear mixed additive model in which each marker is assigned an effect as a solution of the following equation:

$$y=1\mu +Wq+\mathcal{ }\mathcal{E}$$

where y is the vector of observed phenotypes, µ is the mean of y, W is the genotype matrix (e.g., {0,1,2} for biallelic SNPs), q ∼ N (0, ${I\sigma }_q^2$) is the vector of marker effects, and ε ∼ N (0, ${I\sigma }_e^2$) is the vector of residuals. This model, which is solved in a restricted maximum likelihood (REML) context, assumes that the effects of all loci have a common variance, making it suitable for traits influenced by a large number of minor genes. According to a Bayesian context, Bayesian Lasso assigns different prior densities to marker effects, with a strong shrinkage for regression coefficients of marker effects with small values [78]. Predictive ability values were computed as Pearson’s correlation between the observed phenotypic values and the breeding values predicted by the regression model. We implemented a 10-fold cross-validation scheme repeated 50 times for numerical stability. Reported values are the resulting averages. All regressions and cross validations were implemented using the GROAN R package [79].