An interspecific F₂ population containing 121 individual plants was generated from WY7 and TCX7 parent materials. The female parent WY7 is a germplasm resource introduced from the UK with a medium seed size and dark-purple seed coat color. The male parent TCX7 has a large seed size, with a white seed coat, is of good quality, and is cultivated by the Jiangsu Yanjiang Institute of Agricultural Sciences, China. The F₂ individual plants and their parents were planted in Xueyao, Jiangsu Province, China from 2020–2021, and the F_2:3 plant lines and their parents were planted in Xueyao and Jiuhua respectively, Jiangsu Province, China, from 2021–2022. Each faba bean line was planted in one row of 2.4 m in length, with a row distance of 0.8 m, and plant spacing of 0.2 m. Field management was consistent with local production practices throughout the whole growth period. Ten seed phenotypic traits and five nutritional quality traits of the parents, F₂ individual plants, and F_2:3 families in two environments (Xueyao and Jiuhua) were investigated. Ten plants in each F_2:3 line and their parents were harvested. The seed shape traits assessed were SA, SP, SL, SW, SLWR, and ST. The average indicators of HSW, SA, SP, SL, SW, and SLWR used an automatic seed testing system (SC-A1, Hangzhou Wanshen Detection Technology Co., Ltd., Hangzhou, China). The average values of the 10 thickest parts of the seeds were regarded as ST. The seed coat color traits including SC-R, SC-G, and SC-B were measured by spectrophotometer (YS3020, 3NH, China). Nutritional quality traits PC, StC, FC, LC, and TC were determined using a DA7250 NIR analyzer (Perten Instruments, Hägersten, Sweden) with three replicates.

Statistical analysis of the data, such as frequency distribution, coefficient of variation, standard deviation, skewness and kurtosis analysis, was performed using the ANOVA function of IciMapping 4.2.53. The phenotypic correlation between these traits was obtained by Pearon’s correlation analyses using SPSS software. Ver. 26 (IBM SPSS Statistics, Chicago, IL, USA) and R software (version 3.2.2, http://www.r-project.org).

The total genomic DNA of the F₂ individuals and their parental lines was extracted from fresh leaves using the CTAB method (Doyle and Doyle, 1987). A NanoDrop spectrophotometer (Thermo Fisher Scientific, USA) was used to determine the optical density ratios of OD260/280 (>1.8) and OD260/230 >1.5). A Qubit was used for precise quantification, and gel electrophoresis was used to monitor and assess the quality and contamination of all DNA samples. The Illumina sequencing library was constructed by binding biotin-labeled RNA probes to spliced DNA fragments using restriction enzymes and was sequenced using the China Golden Maker (Beijing) Biotech Co. Clean data were derived from the raw sequencing data after quality control (filter parameters: trimmomatic-0.36.jar PE -phred33 ILLUMINACLIP: fa: 2:30:10:8: true LEADING:3 TRAILING:3 SLIDINGWINDOW: 4:15 MINLEN:100) and then matched to the faba bean transcriptome (Wang et al., 2021) by using BWA software (version 0.7.17) with parameters: MEM -T 4 -K 32 -M -R). Based on the results of the sequence alignment, SNPs from the populations genomic data were detected with GATK (version 4.1.2.0) and filtered with VCFtools (version 0.1.13). The detailed criteria and analysis methods were in accordance with Wang et al. (2021).

The harvested genotypes of the samples were firstly filtered before genetic map construction. Based on the filtered genotypes, for each loci, the individuals were coded as “A” (if same with parent TCX7), “B” (same with the parent WY7), “H” (heterozygous containing 2 alleles from each of the parents) or “missing” (all other scenarios). The discarded loci include 1) the loci which were heterozygous in either parent, and 2) the loci with the missing rate above 80% in the population. This was done by using a python script from Li et al. (2021). The genetic map construction used a similar procedure as Li et al. (2021). Briefly, the coded “ABH” genotype matrix was firstly filtered to discarding distortion loci with the threshold P value = 0.01, and then was fed to Lep-Map3 (Rastas, 2017). The default parameters and a logarithm of odds (LOD) score of 12 were used in Lep-Map3. Linkage groups (LGs) with markers less than 100 were removed, and Kosambi function was applied to covert the recombinant rate into LG length (cM, centi-Morgan).

Based on the genetic map constructed above and the phenotypes from multiple environments, we conducted QTL mapping in 2 programs, i.e. QTL Cartographer 2.5 (Wang S et al., 2012) and IciMapping 4.2.53 (Meng et al., 2015). In QTL Cartographer, CIM (Composite interval mapping) method was used, and the parameters were set up as: control markers = 5, window size = 10.0 cM, walk speed = 1.0 cM, and the LOD threshold was determined by 500 times permutation tests. For the Icimapping program, ICIM (Inclusive composite interval mapping) method was selected, and the flowing parameters were used: “missing phenotype = Deletion”, “mapping step = 1 cM” and “LOD threshold = 1000 times permutation at type I error 0.05”. In the mapping result, VG/VP value reflects the explanation rate of phenotypic variance, and the confidence interval of a QTL was determined by the outermost 2 markers above threshold. The QTLs were named as follows: q + trait abbreviation + chromosome number + QTL number.

The splice junction sequences in the Faba_bean_ 130 K SNP TNGS genotyping platform were searched within the QTL intervals and then mapped to 243,120 unigenes (Wang et al., 2021), which were referred to in order to obtain the candidate genes and their gene annotations.

Sequence of the genes in genetic map alignment with reference genome (https://projects.au.dk/fabagenome/genomics-data) and the candidate genes were visual mapped to the reference genome using TB tools software.

The two parent materials showed significant differences in HSW, SA, SP, SL, SW, SC-R, SC-G, SC-B, FC, TC, StC and LC ( Table 1 ). The statistical results of the phenotypic variations in the seed-related traits among the parents, F₂ populations, and F_2:3 individuals ( Supplementary Table S1 ) suggested that HSW, SA, SP, SL, SW, SLWR, PC, StC, FC, LC, and TC showed continuous variation. The absolute values of skewness and kurtosis were almost less than 1, approximately conforming to the normal distribution, meeting the requirements of QTL analysis ( Figure 1 ; Table 1 ).

Significant Pearson’s correlations (p < 0.01) for the same trait showed a significant positive relationship between the F₂ and F_2:3 populations in Xueyao and Jiuhua ( Supplementary Table S2 ). Phenotypic correlations (p < 0.01) among the different traits are shown in Figure 2 . Seed shape traits, including HSW, SA, SP, SL, and SW, were positively correlated with PC, and negatively correlated with StC. Seed coat color traits, including SC-R, SC-G, and SC-B, were positively correlated with FC and StC and negatively correlated with LC. There was no significant correlation between the seed shape traits and seed coat color traits in this study.

A total of 121 F₂ plants and their parents were genotyped using 130,514 SNPs in the Faba_bean_ 130 K SNP TNGS genotyping platform, showing excellent results, quality, and matching scores ( Supplementary Tables S3 ; S4 ). There were 12,023 SNP-tagged gene microarrays with polymorphism between parents ( Supplementary Table S5 ), and they were successfully genotyped into “A,” “B,” and “H” types in the population. All co-isolated markers were defined as one bin, and 1106 bin markers were used to construct a genetic map containing 6 LGs. The overall length of the genetic map was 1,182.65 cM with an average marker spacing of 0.098 cM. Each LG range was from 157.08–296.82 cM, and the average distance between markers was from 0.079–0.114 cM. LG1 had the largest number of markers with 3,325 SNPs. The smallest gap identified in the map was 0.826 cM, the total number of gaps > 5 cM was 9, and the largest gap was 11.78 cM LG6. Additionally, the ratio of marker intervals < 5 cM for all LGs was > 97% ( Figure 3 ; Table 2 ).

QTL mapping was performed using QTL lciMapping and QTL-Cart CIM, and 65 ( Supplementary Table S6 ) and 50 ( Supplementary Table S7 ) QTLs were identified for all 15 seed-related traits detected in the F₂ and F_2:3 populations, respectively. Together, these two mapping strategies identified 28 overlapping QTLs ( Supplementary Table S8 ). Of these, the QTL intervals observed using the CIM method were usually wider, whereas the intervals from the ICIM method were narrower. Consequently, the results obtained using the lCIM method were used in this study. The genetic effect (the explanation rate of phenotyte variance or VG/VP) of the QTLs detected using ICIM for 15 seed-related traits ranged from 4.90–73.99%, with peak LOD values ranging from 4.48–35.25 ( Supplementary Table S6 ). Among the 65 loci, there were 11 QTLs that were detected for more than two traits ( Supplementary Table S9 ). There were 39 QTLs identified that individually accounted for > 10% of the phenotypic variation ( Table 3 ) and 1 QTL explained < 5% of the phenotypic variation ( Supplementary Table S6 ). A total of 41 and 21 QTLs were found to have positive and negative additive effects, respectively.

The genes in the QTL intervals were screened using the Faba_bean_ 130 K SNP TNGS genotyping platform ( Table 4 ). The results showed that 333 genes and 610 SNPs were detected at 65 QTL intervals. Among the 333 genes, HSW, seed shape, seed coat color, and nutritional quality traits contained 8, 117, 100, and 109 genes, respectively, and 173 genes were functionally annotated by database comparison. A total of 67 candidate genes within the environmentally stable QTL intervals were detected, including 2, 20, 39, 3 genes related to HSW, seed shape, seed coat color, and nutritional quality traits, respectively. The results showed that 213 genes in 41 QTLs explained > 10% of the observed phenotypic variance, and they were further assessed ( Supplementary Table S11 ). There were 6 genes related to HSW within these QTL intervals, and 5 were annotated, including the CCCH-type zinc finger protein and calcium-binding protein. There were 53 seed shape-related genes and 30 genes were annotated, including serine/threonine phosphatase, bHLH transcription factor, calcium-binding protein Ca²⁺/H⁺-exchanging protein, and other functional genes. Seed color-related genes included 39 and 19 genes that were annotated, including ubiquitin-like protein, the WD40 family, and transcription factors. There were 79 genes associated with nutritional quality traits, and 41 genes were annotated, including numerous genes encoding enzymes, functional genes, and some transcription factors.

Sequences of gene in our genetic map were well alignment with the recent published reference genome of faba bean ( Supplementary Table S12 ). It was found that about 60% of the genes in each LG were mapped to the corresponding chromosome. Specifically, LG1–LG6 were assigned to chromosome 1, chromosome 3, chromosome 2, chromosome 5, chromosome 4 and chromosome 6, respectively.

Candidate genes were mapped to the reference genome and most annotated genes were located on other five chromosomes except the chromosome 3 ( Supplementary Figure S1 ; Supplementary Table 11 ). Twenty-five genes were located on chromosome 1L (the long arm of chromosome 1) and 17 genes were located on chromosome 1S (the short arm of chromosome 1). Seven, eleven, one and seven of these annotated genes were located on chromosome 2, chromosome 4, chromosome 5, chromosome 6, respectively. Furthermore, there were also seven genes located on free chromosomes (the unassigned scaffolds that cannot be placed on any known chromosome).

Owing to the rapid development of high-throughput sequencing technologies, sufficient molecular markers can now be obtained to facilitate the mapping of high-density genetic maps and research on map-based gene cloning (Yang et al., 2012; Zhang et al., 2016; Zhou et al., 2018; Gaur et al., 2020; Gu et al., 2020; Sa et al., 2021). The molecular genetic analysis of faba bean is currently lagging in comparison to that of many other crops due to its large genome size (Adhikari et al., 2021). Establishing a reliable linkage map between genetic markers and traits is one of the key approaches to improve molecular breeding without a reference genome (Chapman et al., 2022). In this study, a genetic map of faba beans was constructed using high-throughput genotyping platforms. To date, genetic map construction using microarray chips has been successfully reported in several crops, such as pea (Tayeh et al., 2015), wheat (Liu et al., 2018; Ren et al., 2021), cotton (Gu et al., 2020) and pepper (Cheng et al., 2016). In addition, the 130 K liquid-phase gene chip used in this study was developed using transcriptome data, which contains large-scale information. Furthermore, all marker sequences provided valuable gene information, indicating that this liquid-phase gene chip is an effective and feasible tool to utilize for genetic map construction.

There have been more than 20 genetic maps reported for faba beans. Of these, the genetic map constructed by Carrillo-perdomo et al. (2020) containing 1,728 markers, with a total length of 1,547.71 cM and an average genetic distance of 0.89 cM. To date, one of the two SNP genetic map constructed by Li et al. (2023) had the highest density, containing 5,103 markers, with a total length of 1,333.31 cM and an average genetic distance of 0.26 cM. In the present study, an ultra-dense genetic map was constructed, encompassing 12,023 markers in 6 LGs, with an average distance of 0.098 cM. The number, density, and distribution quality of the new molecular markers was thus significantly higher when compared with previous genetic maps. The presented genetic map only has 9 gaps > 5 cM, and thus, it can be effectively utilized for faba bean gene mapping and MAS breeding.

QTL mapping and the analysis of candidate genes within QTL intervals is an effective strategy to investigate numerous crop traits (Bornowski et al., 2020; Chen et al., 2021), and can contribute to the development of molecular marker-assisted breeding (Torres et al., 2010). The 100-seed weight, seed shape, and nutritional quality of faba beans are all quantitative traits susceptible to environmental influence. To improve the accuracy of QTL mapping for seed traits, QTL analysis was performed in the F₂ and F_2:3 populations in two locations. There was a total of 65 seed trait-related QTLs detected ( Supplementary Table S6 ), of which, 11 were repeatedly detected in different environments ( Supplementary Table S9 ).

Patto et al. (1999) used a genetic map constructed using the F₂ population and found that most of the QTLs related to seed weight were located on chromosome 6 for faba bean. Using the recombinant in bred line (RIL6) population constructed using Vf6 and Vf27, Ávila et al. (2017) identified 5 QTLs for HSW, which were located on 4 different chromosomes. Tian et al. (2018) identified two QTLs for seed weight using an F₂ population derived from Yun122/TF42, which were located on two different LGs. In this study, we identified three QTLs linked to HSW, one at LG4, and two at LG5. These results indicate that faba bean seed weight is controlled by multiple main-effect QTLs. qHSW5.1, one of the three QTLs related to HSW, was also associated with SA, and qHSW5.2 was associated with SA, SP, SL, and SW, which indicated that these two QTLs are also involved in controlling seed shape ( Table 3 ).

Seed shape traits are among the most important factors used to determine seed size. The localization and cloning of seed shape genes are of great importance when aiming to increase crop yield and improve appearance quality (Austin and Lee, 1996; Song et al., 2007; Verma et al., 2015; Cheng et al., 2017; Murube et al., 2020). According to the Gramene website (http://archive.gramene.org/qtl/), more than 400 rice grain shape-related genes/QTLs have been identified through genetic mapping and correlation analysis. However, few studies have reported QTL mapping for the seed shape traits of faba bean, a seed length-related QTL and a seed width-related QTL were identified by Tian et al. (2018), 8 QTLs related to seed length, 9 QTLs related to seed width and 8 QTLs related to seed thickness were identified by Li et al. (2023). In this investigation, 28 QTLs for 6 seed-shape traits were identified using linkage analysis, and most were located on LG5( Table 3 ; Supplementary Table S6 ). Compared to these QTLs reported, those identified as controlling seed shape in this study were new, and could thus be applied to the subsequent fine mapping of seed shape traits and the investigation of related genes in faba bean. qSA5.1, qSLWR6.1 and qST5.1 were stable QTLs explained > 10% of phenotypic variation, while qSA5.1 was also associated with SP, SL, and SW. which indicated that these QTLs can be used for further fine mapping and superior gene discovery of seed shape traits.

Seed coat color is a key factor affecting seed quality (Yoshimura et al., 2012; García-Fernández et al., 2021). Different seed coat colors may have different functions (Debeaujon et al., 2003), and the different seed coat colors of faba beans may also be associated with different nutritional qualities. The results of the correlation analysis among seed traits showed that seed coat color was positively correlated with FC and StC, and negatively correlated with LC. Mendel first proposed that the seed coat color of peas is controlled by a pair of genes and considered a qualitative trait (Myers, 2004), while the seed coat color of soybeans is controlled by multiple genetic loci (Choung et al., 2001), and more than 30 molecular marker loci on different chromosomes that control seed coat color in soybean have been detected (Yuan et al., 2022). However, few studies on the QTLs for seed coat color in faba bean have been reported. WY7 and TCX7, the parents used in this study, have purple and white coats, respectively. A total of 12 QTLs, mainly located on LG1, were detected by quantitative measurement of the SC-R, SC-G, and SC-B. qSC-R1.2 was also located with SC-G and SC-B ( Table 3 ), which could explain the > 50% phenotypic variation. qSC-R1.1 is located with SC-G, and qSC-R1.3 is located with SC-B ( Supplementary Table S6 ). These three QTLs are key objects for further study of grain coat color traits.

The main nutrients in faba bean seeds are protein and starch, with low lipid and fiber content levels, as well as tannin (Zanotto et al., 2020), pyrimidine glucoside, and other bioactive substances (Björnsdotter et al., 2021). QTL mapping for quality traits can help to improve the utilization and value of faba beans. At present, there are relatively few studies on the QTL mapping of quality traits in broad beans. Only five genes that control grain proteins have been identified (Macas et al. 1993b). In this study, 22 QTLs linked to quality traits were detected using SNP markers for the first time, including 7 QTLs for FC, 7 for StC, 4 for LC, 2 for PC, and 2 for TC ( Supplementary Table S6 ). In particular, qFC1.1, qTC1.1, qLC1.1, and qLC3.1 could explain > 20% of the phenotypic variation, and qStC1.1 was also associated with PC ( Table 3 ). These QTLs could thus be used to identify the candidate genes for faba bean quality traits.

To identify candidate genes for seed-related traits in faba bean, we focused on 213 genes within 41 QTL intervals that explained > 10% of the phenotypic variation. According to the results of the functional annotation, 57.28% of these genes had been annotated. Signaling pathways that regulate seed size in plants include the ubiquitin-protease pathway, mitogen-activated protein kinase signaling pathway, transcriptional regulation, G-protein signaling pathway, IKU pathway, and plant hormones (Gnan et al., 2014; Li and Li, 2016; Li et al., 2019). Jayakodi et al. (2023) identified 15 marker–seed size associations, and most prominent signal was located on chromosome 4 within the Vfaba.Hedin2.R1.4g051440 gene. In this investigation, there were 30 genes annotations among the 59 genes linked to HSW and seed shape ( Supplementary Table S11 ). Thirteen of these annotation genes located on chromosome 4 by whole genome sequence alignment. dou_TRINITY_DN52935_c3_g4 and hua_TRINITY_DN119282_c0_g1 encode serine/threonine phosphatase and the transcription factor bHLH, respectively, which are reportedly involved in regulating seed size (Savadi, 2018). dou_TRINITY_DN38848_c0_g1 encodes a CYP gene and CYP is involved in protein folding, signal transduction, and RNA processing (Krücken et al., 2009). There are also two calcium signaling pathway genes, including a calcium-binding protein gene ye_TRINITY_DN120969_c0_g1 and a Ca²⁺/H⁺-exchanging protein gene hua_TRINITY_DN154 119_c0_g2, which may be involved in the Ca signaling pathway to regulate seed development. Other unannotated candidate genes could also potentially regulate seed size.

The seed coat color of plants is affected by numerous factors, but flavonoids are the decisive pigments (Lepiniec et al., 2006). In this study, there were 34 candidate genes for seed coat color, 19 of which were annotated ( Supplementary Table S11 ). Among these genes, the translated product of ye_TRINITY_DN150431_c0_g1 is a ubiquitin-like protein that plays an important role in pigment accumulation (Tang et al., 2015). ye_TRINITY_DN150347_c0_g1 and ye_TRINITY_DN139828_c0_g1 are WD40 family genes, which have been suggested to regulate the formation of proanthocyanidins in seed coats (Shirley et al., 1995; Walker et al., 1999). Furthermore, the other 16 annotated genes and 15 unannotated genes may also be required for the pigment composition of different seed coat colors, but this requires further verification.

In this study, 79 candidate genes were associated with five nutritional quality traits, of which, 41 were annotated ( Supplementary Table S11 ). There were 7 genes for LC, 4 of which were annotated, but no functions related to lipid synthesis and accumulation were reported. hua_TRINITY_DN145176_c0_g1, a crude fiber candidate gene, is a triose-phosphate transporter gene that reportedly affects starch and glucose transport in transgenic tobacco (Häusler et al., 1998). dou_TRINITY_DN53089_c3_g1 and ye_T RINITY_DN155843_c1_g1 are GDSL esterases that may also be involved in fiber metabolism. Condensed tannins, also known as proanthocyanidins, exhibit antioxidant, antibacterial, anticancer, and anti-mutation activities (Gutierrez et al., 2020). The two genes zt-1 and zt-2 are the most studied for controlling tannin content in faba bean (Gutierrez et al., 2006; Gutierrez et al., 2007; Gutierrez et al., 2008). Of the candidate genes related to tannins, dou_TRINITY_DN58315_c1_g1 encodes a bHLH transcription factor gene, which is reportedly involved in the mechanisms of tannin biosynthesis in faba bean (Gutierrez et al., 2020). Other tannin-annotated genes obtained in the target intervals have not been reported in faba bean, and thus may be candidate genes affecting tannin content. Further studies are required to confirm the functions of these genes.

Compared to chromosomes and gene locations of the reference genome, the number of linkage groups in our genetic map was consistent with their respective chromosomes, but there were variations in the order of genes on the chromosome, and about 25% of them were not found in the genome ( Supplementary Table S12 ). Eighty-five candidate genes within the QTL interval were mapped to the reference genome, seven of which were located on the contigs ( Supplementary Figure S1 ). Therefore, a part of contigs on the reference genome can be assembled to the genome of faba bean based on the map constructed in this study, which is conducive to the further improvement of the physical map of faba bean.