A high-density consensus linkage map of white lupin highlights synteny with narrow-leafed lupin and provides markers tagging key agronomic traits

White lupin (Lupinus albus L.) is a valuable source of seed protein, carbohydrates and oil, but requires genetic improvement to attain its agronomic potential. This study aimed to (i) develop a new high-density consensus linkage map based on new, transcriptome-anchored markers; (ii) map four important agronomic traits, namely, vernalization requirement, seed alkaloid content, and resistance to anthracnose and Phomopsis stem blight; and, (iii) define regions of synteny between the L. albus and narrow-leafed lupin (L. angustifolius L.) genomes. Mapping of white lupin quantitative trait loci (QTLs) revealed polygenic control of vernalization responsiveness and anthracnose resistance, as well as a single locus regulating seed alkaloid content. We found high sequence collinearity between white and narrow-leafed lupin genomes. Interestingly, the white lupin QTLs did not correspond to previously mapped narrow-leafed lupin loci conferring vernalization independence, anthracnose resistance, low alkaloids and Phomopsis stem blight resistance, highlighting different genetic control of these traits. Our suite of allele-sequenced and PCR validated markers tagging these QTLs is immediately applicable for marker-assisted selection in white lupin breeding. The consensus map constitutes a platform for synteny-based gene cloning approaches and can support the forthcoming white lupin genome sequencing efforts.

White lupin (Lupinus albus L., WL) is a legume crop cultivated in ancient Greece and Egypt more than three thousand years ago 1 . Cultivars being developed currently are appreciated as a valuable source of protein (38-42% in seeds) 2 , while having moderate seed content of oil (10-13%) with excellent food quality 3 . These attributes make this species valuable for human food and animal feed 4 . Lupins have positive impact on soil fertility through symbiotic nitrogen fixation and efficient mobilization of soil phosphorus, and are extremely beneficial in crop rotations 5 . Elite germplasm resources with strengthened cold tolerance, determinate growth habit and dwarf architecture have been developed [6][7][8] . Moreover, untapped landraces with higher yielding ability than current cultivars await thorough exploitation in breeding programmes 9 . However, despite such a long domestication history and considerable opportunities for crop improvement, WL breeders have to struggle with various challenges hampering worldwide cultivation of the species that include, inter alia, the exploitation of optimal vernalization requirement for floral initiation and the identification of genetic resistance to devastating diseases, such as further analysis, while 3597 markers were retained. The average observed heterozygote frequency for RIL lines for the set of 3597 markers was 0.38%. Marker sequences are provided in Supplementary Table S1. These markers increased 8-fold the total number available for WL, confirming the capability of GBS to expand considerably the marker suite applicable for genetic studies. Likewise, GBS provided thousands of SNP markers in another cool-season grain legume species without a reference genome sequence, pea 50 . Genotyping data for 465 published WL markers 4,11,31,33 were verified by preliminary linkage mapping. Marker phase was revealed to be inverted for 76 markers (75 "Pms" markers and Lup270). Segregation data, Chi-square values and allele frequencies for all markers are provided in Supplementary Table S2.
Consensus WL linkage map carries 3669 sequence-defined markers. Some 4062 markers (3597 developed in this study, and 465 from two previous maps) were subjected to linkage mapping. 196 RILs were assayed, including 189 RILs genotyped in this study and 7 RILs missing new sequencing data but having genotype and phenotype scores from previous surveys. A total of 166 complete multipoint mapping runs were performed (from 5 to 10 per linkage group) to establish the final version of the map. Some 99.7% of the analyzed markers were mapped (4050), using 97.9% as "skeleton" (i.e., contributing to calculation of recombination frequencies), and 2.1% as "attached" (with position assigned approximately in relation to adjacent skeleton markers). The standard deviation (SD) of skeleton relative marker position (in Kosambi cM), calculated for the last five Joinmap outputs, ranged from 0 to 7.27, with an average value of 1.08 (Supplementary Table S3). Markers from two previously published maps 11,31 had higher average SD of relative position than those developed in this study (1.34 vs 1.05). Such a difference may have resulted from the different number of lines genotyped (94 vs 189), as well as from differences in the reliability of genotyping method applied. Twelve previously published markers failed to match any linkage groups. STS markers constituted 90.59% of the consensus linkage map; 98.04% of which were generated in this study ( Table 1). As marker grouping was performed using high independence LOD threshold (11.0) and was followed by several rounds of linkage group reconstructions with the removal of markers showing the nearest neighbour fit/stress values higher than 3.0 cM, only well-fitting markers were incorporated to linkage map. An average distance to markers with distorted segregation was 0.80 (Joinmap 4.1 51 ) or 1.61 cM (DistortedMap 1.0 52 ), whereas an average distance to markers with normal segregation was 0.73 (Joinmap 4.1) or 1.54 cM (DistortedMap). The influence of distorted markers on the linkage map length was revealed to be negligible, therefore we chose to retain such markers in the consensus map.
In total, 25 linkage groups were constructed (ALB01 -ALB25), consistent with the 25 chromosome pairs reported for this species 53 . All linkage groups from two previous maps (hereafter, '2007' and '2013') were matched (Fig. 1). No linkage group split was observed for the 2013 map 31 but three splits were observed relative to the  (Fig. 2). Some NLL chromosomes revealed fragmented macrosyntenic anchors to three WL linkage groups, such as NLL-04, NLL-08 and NLL-11. However, no WL linkage group showed large collinearity links to more than two NLL chromosomes (Fig. 2, Supplementary Table S8). The synteny pattern may indicate that WL is more similar to the ancestral form than NLL. Given that there appeared to be a whole genome triplication event around 24.6 mya 41 , it seems likely that NLL has undergone more chromosome fusions than WL. This would be analogous to A. thaliana whose 5 chromosome pairs appears to have fused from 8 ancestral chromosome pairs, hence the value of the more un-rearranged genome configuration of A. lyrata in comparative genomics 54,55 . This would increase the value of WL as a model, especially once its genome becomes available.
Vernalization responsiveness in WL is under polygenic control. QTL mapping was performed using three methods: interval mapping, IM (MapQTL 6 56 ), composite interval mapping with background markers, CIM (Windows QTL Cartographer 2.5) and genome-wide composite interval mapping, GCIM (mrMLM package 57 ). Four QTLs were identified for vernalization responsiveness in the Polish field experiment by IM, three by CIM and three by GCIM (Table 2). Four QTLs were found in the Polish greenhouse experiment. All methods provided similar localizations and additive effect estimates for three major QTLs. Four to five QTLs were found in the Australian field experiment by these methods. Loci of three QTLs were consistent across all experiments (nonv05, nonv16 and nonv15 traits on linkage groups ALB02, ALB13 and ALB16; Fig. 3), suggesting the involvement of the same set regulatory factors (genes or heritable epigenetic marks). The explained phenotypic variance for each QTL ranged from 6.5% to 41.9%. The positions of two major QTLs (on linkage groups ALB02 and ALB13) correspond to coordinates of those previously published 11,31 . Phenotype scores used for QTL mapping are provided in Supplementary Table S9 and QTL mapping results are provided in Supplementary Table S10. Kiev Mutant has been reported to carry a single recessive early flowering gene, brevis 1 . However, instead of a 1:1 segregation ratio expected of a single gene in a RIL population, continuous variation for flowering time was observed in RILs here and in other crosses involving Kiev Mutant 10 . In the absence of vernalization, we found the most influential QTL located on ALB13 and controlling up to 41.9% of phenotypic variation (Fig. 3). Taking into consideration the IM/CIM LOD distribution across the linkage group ALB13 and the resolution of interval mapping methods, the presence of a pair of quite closely linked QTLs might be hypothesized. It is possible that the previously announced gene brevis may relate to this gene (or genes). However, it is clear from the occurrence of multiple QTLs for flowering time that early flowering in Kiev Mutant is under polygenic control, which supports previous observations of low frequencies of early flowering plants in F 2 and downstream progenies of Kiev Mutant × late lines crosses (~4%) 10,11 . Thus, while the brevis allele may be influential, it comprises a relatively modest proportion of the early flowering trait in Kiev Mutant. In addition to brevis, other recessive genes conferring early flowering phenotype (floridus, festinus, and contractus) were described in WL germplasm; however, their relation to Kiev Mutant or P27174 genes remains unknown 58 . Recent investigation of a world WL collection identified several accessions with considerably reduced time to flowering, including Start (cultivar) and P28283 (French breeding line) 29,59 . Crossing of Start and P28283 with late flowering lines revealed that early flowering in these lines is controlled by two complementary dominant genes Ef1 and Ef2, different than those of Kiev Mutant and Ultra 59 . Besides QTLs overlapping in both Polish and Australian phenotyping experiments, one additional minor effect QTL in the Australian experiment was identified by all methods (nonv05_5). Such an observation may result from differences between environmental conditions between these two locations. Przebędowo experienced average minimum and maximum temperature about 3 °C higher than Perth during phenotyping surveys, accompanied total precipitation lower by 32-42% (Supplementary Table S11). White lupins are considered as facultative long day plants, accelerating flowering under long photoperiod both in controlled environment and field conditions 60 . Theoretical day length based on geographic latitude was approximately 5-6 hours longer in Przebędowo than in Perth. Thus, Polish conditions were more favorable for plants requiring long day to initiate flowering than those recorded in Perth. Days-to-flowering trait may be responsive to plant density. It was reported for A. thaliana The Polish experiment based on vernalized plants revealed the existence of only one weak QTL confirmed by IM and CIM methods and located in different region than the QTLs found without the vernalization treatment. Such an observation strongly suggests that the major genes influencing induction of flowering in non-vernalized experiments were those involved in the vernalization pathway. As vernalization is the main factor accelerating flowering of wild lupins in general 12 , the presence of non-overlapping QTLs between vernalized and non-vernalized plants is an expected outcome. For all QTLs, tightly linked, sequence-defined markers were identified in the study. NLL genome regions syntenic to those carrying WL vernalization responsiveness QTLs encode homologs of known flowering regulatory genes. All QTLs identified for vernalization responsiveness of WL have been located in regions showing shared collinearity to the NLL genome. Moreover, corresponding NLL regions revealed also preserved collinearity to reference legume genomes: Arachis duranensis, Glycine max and Medicago truncatula (Fig. 4). Analysis of syntenic regions revealed that none of the WL early flowering QTLs match the region of chromosome NLL-10 recently shown to contain the vernalization independence Ku locus 39,42 ( Table 2, Supplementary Table S8). The NLL Ku locus encodes a FLOWERING LOCUS T homolog, LanFTc1. NLL vernalization independence is thought to be conferred by 1.4 kb insertion/deletion in the LanFTc1 promoter region. Based on synteny analysis, WL vernalization independence appears not to be controlled by the same FT homologue as in NLL 42 .
It has been hypothesized that vernalization pathways evolved in the temperate Cenozoic era when global cooling occurred and differed between plant families which were already separated at that time [66][67][68] . Whole-genome triplication, which could have provided novel copies of regulatory genes, putatively occurred in the genistoid lineage ~24.6 million years ago (Mya) 41 . Considerable decreases in global temperature are estimated to have happened near the Eocene/Oligocene boundary (~36 Mya), in the Miocene (~15 Mya) and near the end of the Pliocene (~2.5 Mya) 69 . Therefore, vernalization pathways in WL and NLL might have evolved independently to each other, but the divergence time between these two species has not yet been estimated.
Anthracnose resistance of WL is controlled by two major QTLs. Two major QTLs (confirmed by all three methods) were revealed for the 2004 anthracnose resistance experiment and four for the 2005 experiment (Table 3; Fig. 3). Tightly linked flanking markers were identified for all these QTLs. Two major QTLs explaining jointly more than 40% of phenotypic variance were localized at highly overlapping positions in both experiments, indicating the involvement of two heritable factors located on linkage groups ALB02 and ALB04. These two QTLs correspond to two loci identified previously in the 2007 map on linkage groups LG4 and LG17, respectively, which previously were found to explain 31% and 26% of trait variance 11 . Due to high fragmentation of two preceding maps, no linkage between anthracnose resistance and vernalization responsiveness was observed 11,31 . However, improvement of linkage map resulted in merging several linkage groups, including those carrying major QTLs for early flowering (ALB02, ~4.7 cM and ~104.7 cM) and anthracnose resistance (ALB02, 0.4 cM). Combining these traits by conventional breeding proved difficult and crop improvement was not achieved at the desired level with the use of Kiev Mutant as early flowering donor and Ethiopian landraces as the source of anthracnose resistance 10  Ethiopian accessions represent a unique genepool 70 . Replacement of Kiev Mutant by French germplasm in WL breeding resulted in milestone achievement in WL crop improvement, as no linkage between these two traits was observed 29 . The present study provides a genetic explanation for this phenomenon, as well as delivers molecular markers to track both Kiev Mutant early flowering alleles and P27174 anthracnose resistance genes in breeding progenies.
Polygenic control of anthracnose resistance is uncommon. Resistance to Colletotrichum spp. in legumes is usually conferred by single dominant genes, recognizing narrow range of pathogen races and following the gene-for-gene model 71 . RCT1 in Medicago truncatula is a single dominant gene acting directly against C. trifolii  LanrBo in line Bo7212 44,46,74 . The localization of WL anthracnose resistance QTLs in blocks of conserved sequence collinearity enabled us to address the positions of these NLL loci (Fig. 4). None of WL anthracnose resistance QTLs matched NLL genome regions carrying Lanr1 or Anman genes, providing high and moderate levels of resistance to C. lupini in NLL 41,46,49,75 . QTL antr05_3 may hypothetically correspond to multi-strain C. lupini resistance LanrBo locus in NLL, located in the middle of linkage group NLL-11 74 . However, intervals between anchor markers flanking LanrBo in Bo7212-derived linkage map are too large to draw any firm conclusions. Analysis of the NLL genome region showing collinearity to WL QTLs provided two candidate genes for antr04_1/antr05_1 locus, encoding LRR protein kinase (Lup026475 and Lup026461). No candidate for antr04_2/ antr05_2 was found.

Major low alkaloid genes in WL and NLL are different. Alkaloid content was observed as a binary
trait (bitter or sweet) and segregated in a 1:1 ratio (Chi-square P = 0.64) consistent with Mendelian expectations for a single locus. It was mapped by linkage mapping in the same way as the molecular markers. Alkaloid locus was localized with high confidence LOD value of 37 in the linkage group ALB18 (Table 4, Supplementary  Table S3). We also observed no unexpected double crossovers in relation to adjacent molecular markers, which demonstrates the unambiguous position of this Mendelian trait locus. One co-segregating and two tightly linked flanking markers were identified for this locus, which refers to pauper gene, widely introduced to modern cultivars 4 . White lupin reveals quite complex pattern of minor alkaloids but a clear influence of domestication on the major alkaloid component (lupanine) was observed 76 Besides pauper, several other natural recessive low alkaloid mutations in WL germplasm were described, including exiguus, mitis, nutricius and reductus genes. Exiguus and nutricius were also bred into certain cultivars 77,78 . Three natural low alkaloid genes were early identified in NLL, namely depressus, esculentus and iucundus 79 , however only iucundus was exploited for NLL breeding 1 . NLL genome region encoding iucundus gene (NLL-07) revealed conserved synteny to WL linkage group ALB03, whereas pauper gene was localized on ALB18. This finding is in line with recent quantitative and qualitative seed alkaloid content assays, which revealed the partly different pattern of alkaloid compound variation in these species and, occasionally, a different major component influencing total alkaloid content of these species, i.e., lupanine in WL, and lupanine or 13-hydroxylupanine in NLL 14,76,80 . Moreover, it is hypothesized that pauper and iucundus genes differ by function because these species differ by lysine profiles among wild and sweet accessions 81 .
The WL pauper locus matches the L. angustifolius genome region encoding sequences annotated as HXXXD-type acyl-transferases (Lup021583 and Lup021586). Lup021586 gene reveals 100% nucleotide identity to LaAT gene (AB581532.1). LaAT belongs to BAHD acyl-CoA dependent acyltransferase superfamily and was demonstrated to be highly expressed in quinolizidine alkaloid-producing plants but undetectable in the sweet ones 82 . As the quinolizidine alkaloid biosynthesis pathway is only partially elucidated, no candidate gene for pauper has been identified so far 81 . Verification of hypothetical involvement of WL LaAT homolog in low alkaloid profile requires extensive studies involving differential transcriptome and metabolome profiling of mapping population parents, representative subset of RILs as well as selected sweet and bitter accessions.

WL Phomopsis stem blight resistance QTLs do not match NLL Phr1 and PhtjR genes.
Localization of QTLs carrying hypothetical genes conferring increased level of WL D. toxica resistance was tentatively estimated by incorporating data on published flanking markers 34 . Numerous skeleton markers were mapped within the ranges of Phomopsis stem blight resistance QTLs, but the resolution of the previous map 31 was too low to provide narrow flanking DNA landmarks directly targeting QTL loci. Therefore, the presented positions of D. toxica resistance loci (Table 4) should be considered as preliminary, tentative evaluations requiring further assessment. Several donors of Phomopsis stem blight resistance exist in WL collection, including wild accessions from Turkey (P27664), Syria (P27840, P28096) and Ethiopia (P28507, P27174) as well as cultivars from Australia (Andromeda, WALAB2008), France (Lublanc), Portugal (Madeira), Ukraine (Borki), and Germany   84,85 . The approximate positions of Phr1 and PhtjR genes in NLL genome assembly were recently determined 41,48,49 . Interestingly, no collinearity link between any of WL Phomopsis stem blight resistance QTLs and NLL regions carrying Phr1 and PhtjR genes was identified in this study ( Table 4).

Validation of SNP genotypes for PCR-based genotyping.
Attempts were undertaken to implement early flowering and anthracnose resistance QTL-targeting markers in plant breeding using PCR amplification and sequencing (Table 5). Moreover, to supplement this analysis Kiev Mutant and P27174 transcriptome assemblies were aligned to confirm particular SNPs. Amplicons were obtained for all 17 markers analyzed. Expected Kiev Mutant and P27174 sequences were revealed for all PCR product except one (TP253114). Positive verification was obtained for markers matching all major QTLs, namely antr04_1/antr05_1 and antr04_2/antr05_2 for anthracnose resistance, as well as nonv05_1/nonv15_1/nonv16_1, nonv05_2/nonv15_2/nonv16_2, nonv05_3/ nonv15_3/nonv16_3 and nonv05_4/nonv15_4/nonv16_4 for early flowering. Markers for some minor QTLs were also investigated and positively verified, including antr04_3, antr05_3 and antr05_5 for anthracnose resistance as well as nonv05_5 and vern15_1 for early flowering. Marker sequences were deposited in GenBank under accession numbers: MF630883-MF630916.

Conclusions
GBS was confirmed as a breakthrough technique for inexpensive genotyping of WL, a species without a reference genome sequence, making available nearly 3,600 markers with normal segregation. The WL consensus map that we generated constitutes a platform that can enable synteny-based gene cloning approaches and can support future WL genome sequencing efforts. A second important finding of our study is the occurrence of high collinearity between WL and NLL genomes. Using this new understanding of genome collinearity, we showed that four key traits are controlled by different regions in WL and NLL. A third, key finding of this study is the different control of early flowering and anthracnose resistance in L. albus compared with reference L. angustifolius accessions on the grounds of different number of genes involved (several QTLs vs single genes), different inheritance patterns (recessive vs dominant) and putatively different genes or homologs involved (visualized by lack of matching of corresponding trait positions between these genomes despite well-conserved synteny). This has considerable implications for our ability to exploit the genetic information made available by a close, genome-sequenced species such as NLL for the purpose of WL improvement, and reinforces the need for genome sequencing of WL rather than relying on NLL as a proxy genome. Finally, this study served the WL breeding community by making available a batch of allele-sequenced markers tightly linked to QTLs for vernalization independence, anthracnose resistance, low alkaloid profile and Phomopsis stem blight resistance, whose application can be devised immediately for marker-assisted selection of the crop.

Methods
Plant material. Genetic Table S11). Each line was grown in 2 m rows (5 cm between plants in row and 30 cm between rows). Flowering was recorded when 50% plants in the row had one or more open flower on the main stem inflorescence 11,31 . The scale from 1 (the earliest) to 6 (the last) was used. Due to high level of similarity (average variation between these two sets was 0.29), observations from Australian surveys were averaged (arithmetic mean of all 2004 and 2005 repeats) and used as one input, "nonv05" (190 RILs). Different number of RILs with phenotyping data results from two factors: seed availability and plant survival during experiments.
Anthracnose resistance. Anthracnose resistance profiling was performed in two experiments (Perth, years 2004 and 2005), involving three independent replicates. RILs and parental lines were grown in 2 meter rows. Each line had three replicates. Inoculum seedlings (Kiev Mutant) were maintained in a glasshouse at 20 °C. Fourteen days after sowing, plants were sprayed with spore suspension (10 5 /mL) of C. lupini isolate 96A4 (IMI375715) and incubated for 24 h. Ten days after inoculation, infected seedlings were transplanted into disease nurseries. Disease assessment was done when plants had begun to senesce. The percentage of plants with anthracnose symptoms on the main stem or lateral branches were recorded three times as continuous variables employing a scale from 1 (highly resistant) to 5 (extremely susceptible) (Thomas and Sweetingham 2004, Yang et al. 2010). Only partial results were used before due to linkage map constraints 11,33 . Here, the whole dataset was exploited. Resistance scores from three repetitions were averaged for each year and implemented as series "antr04" (152 RILs) and "antr05" (191 RILs).
Alkaloid content. Alkaloid content in 189 RILs and two parental lines was scored using two techniques (Dragendorff paper test and UV fluorescence) and encoded as a binary trait (bitter or sweet). Methods and partial results were presented elsewhere 4 . Here, we incorporated the whole dataset developed in that study.
Transcriptome profiling of Kiev Mutant and P27174 lines. Seeds of parental lines were vernalized as described above. RNA was extracted with TRIzol reagents (Invitrogen) using the method optimized for lupins 40 from young leaves, floral buds and developing pods and then pooled in equal measure. The quantity and quality of total RNA were measured using Qubit (Invitrogen) and BioAnalyzer (Agilent) assays. HiSeq. 2000 sequencing of 100PE cDNA libraries was performed by Macrogen Korea Inc. Transcriptome assemblies were generated by Macrogen using the Trinity assembler 86 . . Evaluation of quality scores was done using fastqc (http://www.bioinformatics.babraham.ac.uk/projects/ fastqc/). Demultiplexing was performed by fastq-multx 87 . SNP calling was executed in Tassel UNEAK 88,89 (parameters in Supplementary Table S12). Genotype scores were obtained in ad-hoc likelihood based SNP caller using the estimated sequencing error rate, 0.01; the minimum number of reads per locus, 2; and the minimum ratio of called genotype likelihood over all candidate likelihoods (a measure of the SNP calls reliability), 0.796 (the minimum value to accept a 2-reads locus). A missing value per marker threshold of 0.6 was applied. As attempting to impute missing data points in this condition would be very prone to errors 90 , no imputation was forced. Genotype codes were assigned as follows: a, P27174; b, Kiev Mutant; h, heterozygote; -, no data. A Chi-squared test of segregation distortion from expected for RIL 8 a/h/b ratio (49.61/0.78/49.61) was applied. Markers with P value below 1E-07 were discarded from further analysis.
Construction of the consensus linkage map. Segregation data from two previous maps 11, 31 were imported to Map Manager QTXb20 91 and manually inspected. To match all three datasets and phenotyping surveys, a consensus list of 196 RILs was established. Corrected segregation files, together with those developed in this study, were imported to JoinMap 4.1 51 . After grouping under independence LOD of 11.0 multipoint mapping was performed (parameters in Supplementary Table S12). Markers with nearest neighbour fit and nearest neighbour stress values higher than 3.0 cM were removed and mapping was repeated as many times as required to maintain these values below 3.0. Four additional mapping runs were performed using the established set of markers to calculate standard deviation of relative marker position. Markers removed during recombination frequency optimization procedure, were localized in the linkage map using Map Manager QTXb20 (map function Kosambi, search linkage criterion P = 1e-3) and classified as "attached". Their approximate positions were calculated by linear interpolation of adjacent skeleton markers. Additionally, DistortedMap 1.0 52 was used to analyze genetic distances between distorted markers.
Comparative mapping to NLL genome. Marker sequences were aligned to WL transcriptome assemblies, both the reference (LAGI01) 30 and those developed in this study (Kiev Mutant and P27174), allowing one nucleotide mismatch per marker and one lacking nucleotide per alignment. One best hit per marker was allowed. A Fasta file with marker sequences replaced by assigned transcripts (if applicable) was used for comparative mapping by BLAST 92 to the NLL genome assembly 41 (Supplementary Table S12). Based on the results of comparative mapping to NLL genome, minor alterations in marker order within tight microsyntenic clusters were applied to establish a final version of consensus map.
Quantitative trait loci mapping. Interval mapping 56 was performed using MapQTL 6 (Kyazma, Wageningen, Nederlands) (Supplementary Table S12). In order to determine the significance threshold of the LOD score under the null-hypothesis (no QTL present), a permutation test (1000×) was performed. Based on the results of this test (Supplementary Table S13), the LOD threshold of 3.5 was used for QTL determination. QTLs were localized at positions with the highest LOD values. QTL boundaries were assigned according to LOD values: outer, 2.0 below maximum; inner, 1.0 below maximum. Composite interval mapping was performed in Windows QTL Cartographer V2.5 (North Carolina State University, Raleigh, USA) using 20 background control markers window size 10 cM and walk speed 0.5 cM. Genome-wide composite interval mapping was executed in the mrMLM package 57 . SNP genotypes of candidate markers linked to the QTLs were validated by using Sanger sequencing of amplicons for further PCR-based genotyping. Primer sequences and annealing temperatures used for PCR during SNP validation were provided in the Supplementary Table S14.
Visualization. Linkage groups were drawn in MapChart 93 . Sequence collinearity blocks were visualized using Genome Synteny Viewer 94 , Circos 95 and Legume Information System 96 .
Data availability. All data generated during this study are included in this published article, its Supplementary Information files and in public repository (NCBI short read archive PRJNA380248 and accessions MF630883-MF630916).