Expansion of the phosphatidylethanolamine binding protein family in legumes: a case study of Lupinus angustifolius L. FLOWERING LOCUS T homologs, LanFTc1 and LanFTc2

Chromosomal localization of LanFTc genes

In the narrow-leafed lupin genome two copies of FTc homolog exist, LanFTc1 assigned to linkage group NLL-10 and LanFTc2 localized in NLL-17 [13]. Based on the sequence annotation, BAC clones carrying these genes were grouped into two contigs, 12 in contig 1 carrying LanFTc1 and 5 in contig 2 carrying LanFTc2. To visualize their chromosomal localization, all clones derived from these contigs were used as molecular probes and subjected to fluorescent in situ hybridization (BAC-FISH) to mitotic chromosomes of L. angustifolius. Five BACs from contig 1 (006C24, 015A19, 075P11, 082M07, 133N08) and one from contig 2 (042F24) gave single locus signals, whereas the remaining ones produced repetitive signals dispersed over numerous chromosomes. BAC clones hybridizing to single loci were used in various combinations in two-colored BAC-FISH experiments. To verify BAC localization, available chromosome-specific landmarks for linkage groups NLL-10 (clones 057K22 and 077C13) and NLL-17 (003B18 and 136C16) [23–25] were included in the survey. In total, 15 BAC pairs were tested: 5 for contig 1 itself, 4 for contig 1 vs 2, and 6 for contig vs linkage group (Table 1). BAC clones originating from the same contig or linkage group produced overlapping signals on one pair of homologous chromosomes whereas those from different contigs or linkage groups yielded single locus signals on two different chromosome pairs (Fig. 1). BAC clones which hybridized to single loci on mitotic chromosomes were considered as new cytogenetic markers. Together with those already published [23–27], they constituted physical anchors for the integration of genetic and cytogenetic maps of the L. angustifolius genome.

	006C24, contig 1	015A19, contig 1	075P11, contig 1	082M07, contig 1	057K22, NLL-10	042F24, contig 2
015A19, contig 1	Y	-	-	N	Y	N
082M07, contig 1	Y	N	Y	-	-	N
133N08, contig 1	-	Y	Y	-	Y	N
042F24, contig 2	N	N	-	N	-	-
003B18, NLL-17	-	-	-	-	-	Y
136C16, NLL-17	-	-	-	-	-	Y
077C13, NLL-10	-	-	Y	Y	-	-

Y: BAC-FISH signals for clones from this pair were observed in the same chromosome pair

N: BAC-FISH signals for clones from this pair were observed in different chromosome pairs

Anchoring L. angustifolius genome scaffolds to the updated reference linkage map

The reference map of the narrow-leafed lupin genome [17] was supplemented with recently published BAC-end sequence (BES)-derived markers from contig 1 (133N08_5, 080K01_3D, 008L15_5, 082M07_3, 130O02_5) and 2 (137O24_5, 071N21_3, 042F24_5, 092M12_5, 042F24_3) [13]. Moreover, the existing BES markers from linkage groups NLL-10 (077C13_3, 057K22_3F2) and NLL-17 (107M16_3, 111G03_5, 024F12_5D, 111L22_5, 141C03_5D, 003B18_3) [23–25] not previously incorporated to the reference map [17] were mapped here. LanFTc1 and LanFTc2 regions differed in their estimated recombination frequencies. Between all BAC and gene-based markers from contig 1 no recombination in RILs was observed, whereas for such markers from contig 2 two recombination events were observed (see Additional file 1). 44 sequence-defined markers from linkage groups NLL-10 and 59 from NLL-17 were subjected to repeat masking and L. angustifolius genome mapping. Markers from LanFTc1 and LanFTc2 loci were not included in this alignment as there were BAC clone sequences available. High confidence alignments to genome sequences were constructed for 29 markers from the NLL-10 and 37 from NLL-17. As some redundant or closely linked molecular markers tagged the same scaffolds, the numbers of unique scaffolds anchored to particular linkage groups were lower: 28 for NLL-10 and 31 for NLL-17 (see Additional file 2). The total length of sequence assigned to the linkage map was calculated as 757 164 nucleotides for NLL-10 and 1 226 451 for NLL-17.

Assembly and functional annotation of superscaffolds carrying LanFTc genes

BES markers from both BAC contigs were aligned to their respective BAC sequences as well as to the narrow-leafed lupin genome assembly. 23 BES markers from contig 1 were anchored in clone 133N08 whereas one BES marker mapped to scaffold28512 (accession KB416309.1). The 5′ part of 133N08 was found to overlap with the scaffold14655 (KB408211.1). 8 BESs from contig 2 were mapped to clone 137O24 but two remained unlinked. 5′ part of 137O24 overlapped with the scaffold72960 (AOCW01133146.1). Based on these results, superscaffolds were constructed, one for LanFTc1 encompassing 208 734 nt, and the other for LanFTc2 carrying 189 931 nt.

Hidden Markov Model gene prediction based on reference protein sequences identified four and six genes in the superscaffolds, respectively (see Additional file 3). These genes included, besides previously annotated LanFTc1 and LanFTc2, apyrase, acetyl-CoA carboxylase biotin carboxylase subunit, ethanolamine-phosphate cytidylyltransferase, galacturonokinase, and four not functionally characterized sequences (Table 3). Alignment of predicted coding sequences (CDSs) to the transcriptome assemblies of L. angustifolius [17], L. albus [41] and L. luteus [42] provided evidence of expression for 9, 9 and 6 genes, respectively (see Additional file 4). Based on this annotation data, gene density for LanFTc1 region was calculated as 1.9 genes/100 kbp, whereas for LanFTc2 as 3.2 genes/100 kbp. Graphical visualization of BAC, BES, linkage marker, repeat and predicted gene model localization in the superscaffolds sequences is presented at Fig. 2.

Superscaffold	Gene No.	Gene (protein) name	CDS length	Reference accession	Fgenesh + score	Positives (%)	Transcriptomic evidence
1	1	DNA-binding domain protein	684	XP_003624562.1	975.49	77.98	ang, alb
1	2	LanFTc1	525	LanFTc1	966.60	100.00	ang
1	3	uncharacterized protein	492	XP_006603704.1	498.86	83.07	ang, alb, lut
1	4	apyrase 2-like	1407	XP_006598985.1	1916.76	95.51	ang, alb, lut
2	1	U-box domain-containing protein	2208	XP_006580498.1	2276.42	87.66	alb
2	2	acetyl-CoA carboxylase biotin carboxylase subunit	1614	XP_003630608.2	2159.24	91.81	ang, alb, lut
2	3	uncharacterized protein	2934	XP_006580493.1	2851.19	79.84	ang, alb
2	4	ethanolamine-phosphate cytidylyltransferase	1191	XP_006580492.1	1552.44	88.81	ang, alb, lut
2	5	galacturonokinase-like	1317	XP_006580491.1	1561.85	88.13	ang, alb, lut
2	6	LanFTc2	528	LanFTc2	978.98	100.00	ang, alb, lut

Abbreviations: ang, L. angustifolius; alb, L. albus; lut, L. luteus

Insight into micro- and macrosynteny of LanFTc regions

Both superscaffolds with low complexity regions and repeats masked by RepeatMasker and Censor were aligned to genome sequences of nine legume species: A. duranensis, A. ipaensis [30], C. cajan [31], C. arietinum [32], G. max [33], L. japonicus [34], M. truncatula [35], P. vulgaris [36] and V. radiata [37]. Both regions revealed distinct microsyntenic patterns (see Additional files 5 and 6 where coordinates and statistics of alignments were provided). Superscaffold-1 showed links of sequence collinearity to A. duranensis chromosome 10 (region between 72.7 and 73.5 Mbp), A. ipaensis chr. 6 (113.2–114.2 Mbp), C. arietinum chr. 3 (26.3–26.5 Mbp), C. cajan scaffold 54 (0.7–1.0 Mbp), G. max chr. 16 (3.6–4.1 Mbp) and 19 (35.6–36.6 Mbp), L. japonicus chr. 1 (48.8–49.9 Mbp), M. truncatula chr. 7 (32.8–33.0 Mbp), P. vulgaris chr. 1 (20.5–22.8 Mbp), and V. radiata scaffold 7 (2.7–3.3 Mbp) (Fig. 3). In all these regions but C. cajan and L. japonicus one copy of FTc gene was identified. In the C. cajan region showing synteny to LanFTc1 no FT gene was recognized whereas in the L. japonicus syntenic region one copy of FTa gene was annotated. Superscaffold-2 was found to have higher level of synteny to genome regions of all studied species than superscaffold-1, however, the patterns of preserved sequence collinearity did not include the 3′ end of the superscaffold-2, containing the LanFTc2 gene (Fig. 4). Syntenic links observed for this scaffold matched different chromosomes than those found for superscaffold-1. Namely, they were as follows: A. duranensis chr. 3 (60.3–99.6 Mbp), A. ipaensis chr. 3 (70.9–85.5 Mbp), C. arietinum chr. 6 (4.1–4.4 Mbp), C. cajan chr. 5 (1.3–1.5 Mbp), G. max chr. 5 (40.2–40.3 Mbp) and chr. 8 (2.2–2.3 Mbp), L. japonicus chr. 4 (40.7–40.8 Mbp), P. vulgaris chr. 2 (46.5–46.7 Mbp), V. radiata chr. 7 (52.8–53.0 Mbp). No copy of FT gene was annotated in any of these syntenic regions.

Comparative mapping of sequence defined markers from linkage groups NLL-10 and NLL-17 resulted in the identification of complex networks of syntenic links to several chromosomes in every species studied (see Additional file 7). Regions syntenic to NLL-10 included those showing microsynteny to superscaffold-1. LanFTc1 gene appeared to be located in the middle of a large conserved region, marked by 9 linkage map loci spanning a distance of 15.3 cM. Despite some inversions and/or insertions/deletions, the arrays of macrosyntenic links observed in this section between sequences from L. angustifolius and studied species were visually similar to each other.

The most conserved arrays of collinearity between NLL-17 and sequenced genome sequences were observed for A. duranensis chr. 1, 3 and 5, A. ipaensis chr. 1, 3 and 5, C. arietinum chr. 6, 7 and 8, C. cajan chr. 6, G. max chr. 1, 5, 8 and 11, L. japonicus chr. 2 and 4, M. truncatula chr. 5 and 8, P. vulgaris chr. 2 and 3, and V. radiata chr. 7 and 11. Those arrays included also all regions (but C. cajan chr. 5) matching superscaffold-2 LanFTc2 (Fig. 5). The syntenic block adjacent to LanFTc2 gene was flanked by markers 111G03_5 (53.1 cM) and 042F24_5 (73.5 cM). It should be noted that the superscaffold-2 was located at the edge of the large syntenic region (73.0–74.1 cM), with a LanFTc2 gene itself excluded from this pattern of collinearity.

Phylogenetic inference of PEBP-family sequences in sequenced legume species

The majority rule consensus tree based on the alignment of 122 PEBP-family sequences (Figs. 6, 7 and 8) shed light on evidences of numerous duplications, visualized by sister branches of shared nodes. In lineages leading to the following species or clades such duplications were observed: C. arietinum: two duplications of TFL1 genes, G. max: two TFL1, one BFT, FTa and FTc; V. radiata: one FTb; M. truncatula: one FTb; L. angustifolius: one FTa, FTc and MFT, and 2 TFL1, giving three closely related final copies. Duplications occurred also in genomes located at ancestral nodes linking closely related species, like the ancestor of A. duranensis and A. ipaensis: FTa, FTb, and MFT, and ancestor of C. arietinum and M. truncatula: FTa. The topology of the tree confirmed close relationships between sister Arachis species, between Lotus, Medicago and Cicer, as well as between Glycine, Phaseolus, Cajanus and Vigna. It also reflected directly more distant relations of early diverged Arachis and Lupinus to other succeeding clades. To facilitate further sequence analysis, gene names were assigned following the topology of the tree. As M. truncatula and G. max PEBP homologs had already common names in use [8] reference to those was provided (see Additional file 8).

Promoter analysis of FT homologs

The presence of multiple copies of FT in regions with very well preserved sequence collinearity raised a question about the sequence conservation in regulatory regions of these homologs. FT upstream sequences of 8000 nt were extracted from all homologs except for L. angustifolius FTa1a and FTa1b (which are localized in short scaffolds). Annotation revealed no significant similarity to CORE1, CORE2 and S4 elements. Despite numerous alignments, no evidence for conserved collinearity of P2-S3-P1 blocks was identified. Conserved positions and orientations were identified for RE-alpha and CCAAT-box sequences.

Assay of substitution rates and selective pressure of PEBP-family genes in legumes

According to the results of phylogenetic inference, PEBP sequences were grouped to six subclades, namely BFT, FTa, FTb, FTc, MFT, and TFL1. Based on the tree topology, 61 species-specific pairs of duplicated sequences were selected, including those located at sister and quasi-sister branches (i.e. paralogs) and those of the most recent common ancestor origin (i.e. orthologs). The nonsynonymous to synonymous substitution rate (Ka/Ks) ratio analysis revealed that all pairs but C. arietinum TFL1d/TFL1e were under strong purifying selection, with Ka/Ks values from 0 to 0.53 (see Additional file 9). The outlier obtained for C. arietinum TFL1d/TFL1e pair may result from the incorrect assembly of C. arietinum TFL1e gene as this sequence was derived from three short contigs (Cascaffold1275, Cascaffold127513210 and Cascaffold127511243) and such a construct may be artificial.

FT-like genes are located in regions of shared cross-genera legume synteny

The analysis based on linkage mapping, micro- and macrosynteny survey, and phylogenetic inference provided clear lines of evidence to support the hypothesis that LanFTc1 gene, involved in vernalization pathway in Lupinus angustifolius L. [13, 18], is a direct descendant of an ancestral FTc, whereas LanFTc2, not responsive to vernalization, is a more recent duplicate.

The availability of L. angustifolius draft genome sequence [29] and reference sequences of nine legume species A. ipaensis, A. duranensis [30], C. cajan [31], C. arietinum [32], G. max [33], L. japonicus [34], M. truncatula [35], P. vulgaris [36] and V. radiata [37] greatly contributed to lupin molecular studies [23, 24, 26, 28]. The cross-clade comparative mapping of genome sequences provided significant improvements to the phylogeny reconstruction of FT-like and other genes from the PEBP family in legumes. Based on the model A. thaliana and the related crop species Brassica rapa it was demonstrated that synteny can be exploited to order gene models by homology providing alignment targets for mapping-by-sequencing [43]. It is a case of many species with missing full reference sequence but partial genome information available. In narrow-leafed lupin, synteny was harnessed to track the evolutionary history of duplicated chalcone isomerase like genes [26] as well as to select candidate gene family (FT-like) underlying vernalization insensitivity locus Ku, conferring early flowering phenotype [16]. The alignment of the first gene-based linkage map of L. angustifolius to the genome sequence of model legume M. truncatula revealed the existence of regions of conserved synteny in 14 lupin chromosomes [16]. The section of L. angustifolius linkage group NLL-10 containing the Ku locus was syntenic with a region of M. truncatula chromosome 7 containing three homologs of the FT gene. Marker development, followed by linkage mapping of FT-like genes identified by Lupinus vs Medicago synteny, turned out to be a successful approach, and revealed the tight linkage of one of FTc-derived markers with locus Ku [18]. Moreover, a synteny survey supported the identification of affinities within FTa, FTb and FTc genes in M. truncatula and G. max genomes [8]. A synteny based approach was also used to track remnants of ancient duplications in G. max and resulted in identification which collinear relationships among blocks containing CO-like and FT-like genes arose after the whole-genome duplication (WGD) events [44]. In the present study, comparative analysis revealed that despite a complex evolution pattern of PEBP family genes, arrays of collinearity in genome regions carrying FT-like genes were not substantially disrupted. The survey of sequence collinearity links surrounding PEBP genes considerably supported phylogenetic inference of these genes.

Copy number of PEBP-family genes in legumes

The research revealed that a majority of PEBP-family genes are present in multiple copies in analyzed legume species. The PEBP family consists of six subclades, namely FTa, FTb, FTc, MFT, TFL1, and BFT. According to our knowledge, presented research is the first whole-genome survey of PEBP gene family across legume species representing several Papilionoideae lineages. It should be emphasized that the estimation of exact number of copies in particular species depends on two main constraints: the quality of the genome assembly and the availability of comprehensive reference transcriptome datasets. In the legume clade, the number of PEBP family genes has already been evaluated for G. max, M. truncatula, Pisum sativum and P. vulgaris genomes [8, 45, 46]. PEBP gene copy number was extensively studied in M. truncatula, and resulted in identification two FTa and FTb genes and one FTc, TFL1, BFT and MFT gene [11, 47–49]. With the use of the new assembly of the M. truncatula genome we were able to identify one more FTa sequence and two other copies of TFL1. Four FTa, four FTb, two FTc, four TFL1, two BFT and two MFT sequences were identified by BLAST searches of the G. max genome [8]. Those results mostly converge with ours but we found one more TFL1 and one less MFT copies. The lacking MFT sequence, Gm08g05650, appeared only in the first genome assembly which contained also questionable FTa Gm02g07650), two FTb (Gm08g28470, Gm18g53670) [46]. Gm02g07650, Gm08g05650, Gm18g53670 and Gm08g28470 are putatively contig-assembly artefacts of the first draft G. max genome release because all soybean PEBP gene models but those Gm02g07650, Gm08g05650, Gm18g53670 have conserved structure, and all but Gm08g28470 have duplicated counterparts generated from the last WGD [46]. Moreover, other genome-wide screening resulted in detection of all G. max copies reported here except for TFL1 (Gm13g39360) and FTb (Gm08g47820), the loci removed from the new genome assembly [45]. In the present study, novel homologs were extracted from P. vulgaris and M. truncatula genome sequences, which were not identified in previous research [46], namely P. vulgaris FTb Phvul.008G003700.1 and M. truncatula FTa (Medtr6g033040.1 and Medtr7g085020.1), BFT (Medtr0020s0120.1), and TFL1 (Medtr1g060190.1). The molecular and in silico survey of P. sativum PEBP genes resulted in the deciphering of coding sequences of two FTa, two FTb, one FTc, and three TFL1 homologs [8]. The revealed inconsistency in determination of precise number of PEBP homologs in legumes among recent reports is a minor issue which does not undermine the general conclusion that FTa, FTb and TFL1 are present in higher copy number than FTc, MFT and BFT.

Conserved motifs of FT promoter sequences

Typically, promoters in A. thaliana are short with an average length of ~500 bp [50]. However, the FT promoter is much longer and contains four major blocks: A, located from ~0.4 kb to the first codon; B, at ~ -2 kb; ID, at ~ -3.7 kb; and C, at ~ -5.5 kb. Block A contains conserved sequences: CO-responsive elements (CORE1 and CORE2) directly bound by CO, and S1, P1, S3, P2, S4 with unknown functions [51, 52]. Block B carries E-box binding site for basic-helix-loop-helix proteins promoting FT expression under blue light [1, 51]. Block ID has no annotated element, however, it shows length variation in Arabidopsis natural populations caused by at least two independent indel events. There were some differences between these types in fitness of plants grown over winter but not in flowering time [53]. Taking into consideration the position of a large indel in the LanFTc1 promoter (~-2.8 to ~ -4.2 kb) [13], it may correspond to Arabidopsis block ID, however, it shows lack of sequence conservation. Block C contains RE-alpha and CCAAT boxes. CCAAT box is a binding site for the NUCLEAR FACTOR Y (NFY) complex participating in making long distance DNA loop and activation of FT by bringing CO to CORE1 and CORE2 [52]. In general, RE-alpha and CCAAT boxes were found to be present in legume FT promoters at expected positions (at ~ -5 kb to ~ -7 kb in legumes, ~-5.3 kb in A. thaliana) indicating that FT promoter lengths in legumes are at least as big as in A. thaliana. It may also indicate that the NFY-mediated mechanism is preserved in legumes.

The lack of annotation of several important blocks, like CORE1 and CORE2, may relate to evolutionary sequence divergence rather than a true absence. However, it is well known that vernalization pathways evolved in temperate Cenozoic when global cooling occurred and differ between plant families which were already separated at that time [54–56]. As FT is a vernalization-responsive gene, different regulatory elements in FT promoter sequence might be present in legume clade. It was evidenced for A. thaliana that vernalization causes the accumulation of repressive histone modification marks at FLC locus, like H3K27me3, silences this gene and, therefore, unlocks FT expression [57]. Differences in vernalization requirement between winter and summer A. thaliana accessions are conferred by allelic variation at FRIGIDA and FLC loci, causing loss or considerable reduction of function [58–60]. In L. angustifolius different mechanism of vernalization independence exists, based on large deletion in the promoter of FT, a target gene for FLC [13].

Duplicated genes as traces of ancient WGD events

In the present research, novel patterns of species-specific PEBP family gene duplications in all but L. japonicus genomes were identified, as well as remnants of preceding ancient duplications, which are thought to have occurred at ancestral nodes shared by several genera. Are these multiple copies traces of ancient WGDs or appeared independently to these events?

Episodes of WGDs have been frequent through angiosperm history [61]. It was proposed that a WGD event occurred in the common ancestor of all extant seed plants (dated ~319 million years ago, mya), which was later followed by the another WGD in the common ancestor of all extant angiosperms (~192 mya) [62]. Moreover, a phylogenomic approach harnessed to investigate the timing of gene set duplications located on syntenic gamma blocks resulted in identification of traces of old genome triplication (~117 mya), associated with early diversification of the core eudicots [63]. It is anticipated that ancestral WGDs facilitated the diversification of regulatory genes important to seed and flower development and therefore ultimately contributed to major innovations which promoted the rise and eventual dominance of seed plants [64]. Taking into consideration all newly discovered and those well-documented and widely accepted genome duplications, even the relatively small genome of Arabidopsis should carry the traces of at least five WGDs [64]. However, gene pairs formed by duplication usually have a relatively short life span, as some copies may be lost, others pseudogenized, and only a limited number survive after duplication [46]. In the legume clade, an ancient WGD occurred in the progenitor line of papilionoids, the remnants of which can be seen in extant clades including the genistoids (e.g. L. angustifolius), dalbergioids (Arachis spp.), millettioids (P. vulgaris, G. max, C. cajan, V. radiata), galegoids (M. truncatula, L. japonicus, C. arietinum), Xanthocercis and Cladrastis [40, 65–68]. This event was dated to about 44–65 mya and directly predated the divergence of ancient lineages of Papilionoideae [32, 40, 65, 69, 70]. It should be noted that evidence was found for several independent WGDs which occurred near the base of other important legume lineages, i.e. Mimosoideae-Cassiinae-Caesalpinieae, Detarieae, Cercideae and Lupinus clades, dated roughly ~30–55 mya [40]. Lupin WGD is believed to have occurred before the divergence of New World and Old World clades [26, 40]. Some legume WGD events occurred in evolutionary recent times, like 13 mya in G. max and several mya in Arachis [33, 40]. The Ks analysis of 110 duplicated C. arietinum genome blocks showing synteny to four legume (M. truncatula, L. japonicus, G. max, and C. cajan) and two non-legume dicot (A. thaliana and Vitis vinifera) genomes indicated a papilionoid divergence time of 58 mya ago [32]. Galegoid and millettioid clades separated ~54 mya [70]. The survey of galegoid genome sequences using genetic distance–transversion rates at fourfold degenerate sites dated the divergence of C. arietinum from L. japonicus ~20–30 mya and from M. truncatula ~10–20 mya [32]. Estimated dates of major WGD events converges not only with the divergence of the particular lineages but also with the appearance of new PEBP gene copies. Analysis of synonymous substitution rates in soybean genome indicated that the most of PEBP duplicated gene pairs originated from two WGD episodes, dated as 59 and 13 mya, and only few sequences appeared relatively recently, about 3–7 mya [45]. The survey of collinear relationships of homologous G. max blocks containing FT-like genes revealed that FTa and FTc all experienced WGDs as well as tandem duplications [44]. In general, the G. max genome has a very high number of retained duplicate genes, including terminal duplicates (visualized as sister genes in the gene trees), derived from the WGD event. Similar observations were made for Arachis and Lupinus genomes [40].

Contribution of duplications to the evolution of PEBP genes

It is anticipated that the basal clade among PEBP genes are MFT-like genes, which are present in angiosperms, gymnosperms, lycophytes and bryophytes [71]. The most likely scenario is that the ancestral FT/TFL1 genes originated from a duplication of an MFT-like gene after divergence of the basal plant lineage, lycophytes [46] and could have contributed to the radiation of seed plants [72]. The second duplication resulting in the production of FT and TFL1 clades, encoding the proteins conferring antagonistic functions, likely occurred with the evolution of flowering plants and contributed to the development of angiosperm lineage [44, 71]. The TFL1 ancestor underwent two separate duplication episodes in the common ancestor of angiosperms, which created BFT and TFL1lineages, followed by another TFL1 duplication in the lineage leading to the Papilionoideae [46]. The pattern of FT-like duplications is more complex: they were at least three subsequent and two parallel duplications of FT genes during the evolution of legumes (including ancient Eurosids I specific tandem duplication) [46]. It may be summarized that all these duplications resulted in the origin of subclades FTa, FTb and FTc as well as in the multiplication of gene copy number in FTa and FTb clades.

In the present study numerous lineage-specific duplications of legume PEBP-family genes were revealed (Figs. 6, 7 and 8). Gene pair divergence timing based on Ka/Ks analysis presented here (Table 6) strongly evidenced that major MFT and TFL1 lineages originated before the well-documented WGD event around the time of the papilionoid diversification (~157 and ~117 mya vs ~48–65 mya). The estimated date of expansion of legume-specific orthologous copies in FTa and FTb converges with the papilionoid ancestor WGD date (~76 and ~65 mya vs ~48–65 mya). The value of 76 mya obtained for FTa resulted putatively from the existence of remnants of ancient Eurosids I specific tandem duplication in this subclade. The present study provided also novel evidence of lineage-specific WGDs, demonstrated by orchestrated convergence of averaged sister-branch gene pair duplication dates for FTa, FTb and FTc, evaluated as 24, 29 and 24 mya, respectively. Taking into consideration the total number of WGD events in the history of legumes and their ancestral lines, evolutionarily ancient MFT and BFT subclades were subjected to at least three WGDs. Such observations considerably complement the knowledge on the complex evolutionary scheme of PEBP family in legumes and validates the conclusion that the WGD events were major mechanisms launching the divergence of two large PEBP subfamilies, FT and TFL1, and to a lesser extent, BFT and MFT. Thus, despite so many WGD episodes, evolutionary old MFT and BFT subclades did not accumulate as many legume paralogs as FT and TFL1. Indeed, the vast majority of hypothetical WGD-derived gene copies must have been completely annihilated in these subfamilies.

Functional divergence of duplicated PEBP genes

Different evolutionary patterns are possible when duplicated gene copies appear in the genome. Homologous copies of a particular gene may acquire different functions, by several possible mechanisms: pseudogenization, subfunctionalization and neo-functionalization, which can lead to the evolutionary divergence [9, 10, 73]. Indeed, functional diversification of the surviving duplicated genes is a major feature of the long-term evolution of polyploids [9]. Genes from the PEBP family, as accumulated in numerous copies in legume genomes, may be preferential targets of natural selection. Insights from patterns of molecular adaptation based on nucleotide variation among large subset of M. truncatula accessions from natural populations justified the conclusion that existing polymorphism of FT genes in this species has probably been shaped by recent or ongoing positive selection [74]. However, our study provided well-supported evidence of strong purifying selection of PEBP genes in legumes.

Mutations are considered to play key roles in coding sequence evolution. Mutations are powerful mechanisms to facilitate the evolution of new functions as by providing the increased genetic diversity to drive the divergence of the duplicated proteins and creating discrete folding pathways contributing to the conformational diversity of newly-emerged proteins [75]. To exemplify, A. thaliana TFL1 and FT are key controllers of flowering but have opposite functions: TFL1 is a repressor, FT is an activator [76]. Although these sequences are less than 60 % identical, just the substitution of single amino acid (Tyr85/Gln140 in FT and His88/Asp144 in TLF1) is sufficient to functionally convert TFL1 to FT and vice versa [77, 78]. However, analysis of legume and other FT genes revealed that although these sequences are highly conserved in rosid FT genes, the legume FT genes differ from the consensus at 3 to 6 positions, with FTc proteins having a His rather than Gln at 140 [11]. The complete set of amino acids that drove functional divergence of all PEBP subfamilies has not yet been identified. Up to now, 46 amino acid residues that exhibited high conservation with 95 % identity in the alignment encompassing FTa, FTb, FTc, BFT, TFL1 and MFT subclades were identified [46]. About 20 % of these residues have been shown to cause loss of function when mutated [76, 79–81]. Such observations may suggest that these conserved amino acids are the prerequisites of maintaining basic functions of PEBP proteins.

Chromosome localization of LanFTc genes (BAC-FISH)

DNA from clones selected from the cv. Sonet BAC library by LanFTc probe hybridization was isolated using the QIAprep Spin Miniprep Kit (Qiagen, Velno, Nederlands) and verified by PCR using insert DNA template and LanFTc1 and LanFTc2 gene-based primers [13]. Agarose gel (1 %) electrophoresis was performed for 3 h at 9 V/cm to determine the quality of DNA inserts. Size marker GeneRuler™ 1 kb Plus (Fermentas Waltham, MA, USA) was used.

Cytological preparations of mitotic metaphase chromosomes were made from dissected root tip meristematic tissue. Chromosome squashes and the BAC-FISH procedure were performed according to the protocol [25], with few minor modifications. These changes included: probe labelling with digoxygenin-11-dUTP and tetramethyl-rhodamine-5-dUTP by incubation at 15 °C for 110 min followed by inactivation at 65 °C for 15 min using Sensoquest Labcycler (Göttingen, Germany) and hybridization at 37 °C for 22 h.

Linkage mapping of sequence-defined markers

The mapping population comprised 89 F₈ recombinant inbred lines (RILs) developed from the cross combination 83A:476 (domestic) × P27255 (wild-type) L. angustifolius [15]. The reference genetic map of L. angustifolius [17] containing skeleton and redundant markers was imported to MapManager v. QTXb20 [82]. Based on the segregation pattern, new markers were distributed in the existing linkage groups (map function Kosambi, linkage criterion 1e-4). MapChart software [83] was used to draw the Lupinus linkage groups. The approximate positions of attached markers were calculated by linear interpolation of adjacent markers.

Anchoring genome sequence to linkage map

Marker sequences from linkage map groups NLL-10 and NLL-17 were used to screen the collection of L. angustifolius whole-genome shotgun contigs and scaffolds [29] obtained from the NCBI sequence database (Project No. PRJNA179231, assembly version GCA_000338175.1, accessions AOCW01000001–AOCW01191454). The BLAST algorithm was optimized for highly similar sequences (word size, 28; match/mismatch scores, 1/-2; and gap costs, linear). Scaffolds producing alignments with sequence identity value above 99 % were integrated to the map at corresponding marker positions. Alignments with sequence identity value of 95–99 % were manually checked for the possibility of sequencing or assembly errors and appropriate scaffolds were placed on map in case of positive verification.

Assembly of superscaffolds carrying LanFTc genes

BAC-end sequences (BESs) from LanFTc contigs, 24 for contig 1 and 10 for contig 2 [13], and BES-derived marker sequences were aligned to sequenced LanFTc BAC clones, 133N08 and 137O24 (accessions LN851864 and LN851865, respectively) as well as to L. angustifolius whole-genome shotgun assembly [29]. The BLAST (https://blast.ncbi.nlm.nih.gov/Blast.cgi) algorithm [84] was optimized for highly similar sequences (word size, 28; match/mismatch scores, 1/-2; and gap costs, linear). Scaffolds producing alignments with sequence identity value above 99 % to particular BACs or BESs were used to assembly longer, consensus sequences of LanFTc1 and LanFTc2 regions. Based on the alignments, BESs were positioned in the constructed superscaffolds.

Functional annotation of superscaffolds

Repetitive elements were annotated and masked using RepeatMasker Web Server version 4.0.3 with implemented repeat libraries RMLib 20140131 and Dfam 1.4 (A.F.A. Smit, R. Hubley & P. Green, unpublished data, http://www.repeatmasker.orgo). The following options were selected: search engine, cross_match; speed/sensitivity, slow; DNA source, Arabidopsis thaliana. Preliminary masked DNA sequences were then compared to a database of transposable element encoded proteins. A third round of masking was performed using Censor maintained at Genetic Information Research Institute [85] with the following settings selected: sequence source, Viridiplantae; force translated search; mask pseudogenes. The annotation of the genetic elements was based on comparative analyses with known sequences. Integrated, non-redundant, sequences of genomic DNA, transcripts, and proteins in the RefSeq database (http://www.ncbi.nlm.nih.gov/refseq) were examined by BLAST. Moreover, scaffolds were subjected to sequence homology searches against the transcriptome sequences of L. luteus [42], L. albus [41], and L. angustifolius [17]. The following sequence repositories were used: L. albus, http://comparative-legumes.org (gene index LAGI 1.0), L. luteus, http://sra.dnanexus.com/studies/SRP014198/runs (project PRJNA170318, archive SRX159101), L. angustifolius (NCBI accession: GBRP00000000.1). CoGe BLAST algorithm [86] was used with the following parameters: e-value cut-off, 1e-10; word size, 8; gap existence cost, 5; gap elongation cost, 2; nucleotide match/mismatch scores, 1/-2. Sequences producing alignments with the lowest e-values were selected as references for protein-based Hidden Markov Model gene prediction in Fgenesh + [87]. Annotation data were exported to European Molecular Biology Laboratory (EMBL) feature table format and visualised in Geneious v8.1.5 [88].

Micro- and macrosynteny survey

Both superscaffold sequences, masked for repetitive contents and low-complexity regions, were aligned to the following genome sequences: A. duranensis (accession V14167) and A. ipaensis (accession K30076) (Peanut Genome Project), C. cajan [31] (project PRJNA72815, v1.0), C. arietinum [32] (v1.0 unmasked, http://comparative-legumes.org), G. max [33] (Phytozome v9.0 unmasked, http://www.phytozome.net), L. japonicus [34] (pseudomolecules v2.5 unmasked, http://www.kazusa.or.jp), M. truncatula [35] (strain A17, JCVI v4.0 unmasked, http://www.jcvi.org/medicago/), P. vulgaris [36] (v1.0, https://phytozome.jgi.doe.gov/pz/portal.html), and V. radiata [37] (project PRJNA243847). The CoGe BLAST algorithm [86] was used to make sequence similarity analyses with the following parameters: e-value cut-off, 1e-10; word size, 8; gap existence cost, 5; gap elongation cost, 2; nucleotide match/mismatch scores, 1/-2. Sequences producing alignments to numerous loci dispersed over many chromosomes were marked as repetitive and excluded from further analysis. Sequence collinearity blocks were visualized using the Web-based Genome Synteny Viewer [89] and Circos [90].

All marker sequences from linkage groups NLL-10 and NLL-17 [17, 23, 24] were repeat-masked by RepeatMasker and Censor [85]. Masked sequences were then aligned to the L. angustifolius whole-genome shotgun assembly [29] using BLAST algorithm optimized for highly similar sequences. Alignments with sequence identity value below 95 % were discarded, whereas alignments with sequence identity value equal or above 95 % were manually checked for the distribution of mismatches, and if accepted, appropriate markers were replaced by scaffolds. Obtained set of markers and scaffolds with assigned genetic distance positions was subjected to RepeatMasker and Censor [85] masking, followed by CoGe BLAST [86] mapping to sequences of nine sequenced legumes (as described).

Identification of PEBP-family homologs in sequenced legume genomes

Selected A. thaliana reference sequences (FT (AT1G65480); TWIN SISTER OF FT, TSF (NM_118156.1), BROTHER OF FT, BFT (AT5G62040), MOTHER OF FT, MFT (AT1G18100), and TERMINAL FLOWER 1, TFL1 (AT5G03840) were aligned to the G. max genome sequence by custom BLAST implemented in Geneious [88] under the following parameters: e-value cut-off, 1e-6; word size, 7; gap existence cost, 5; gap elongation cost, 2; nucleotide match/mismatch scores, 2/-3. Genome sequence regions, extended by 10000 nt in both directions from alignment locus, were extracted and submitted to Fgenesh + [87] gene prediction using original A. thaliana protein sequences as references. Coding (CDS) and protein sequences obtained from Fgenesh + were compared with those deposited in Phytozome G. max annotated genome assembly at appropriate nucleotide positions (https://phytozome.jgi.doe.gov/pz/portal.html). To facilitate further tracking of these sequences, Phytozome accession numbers were added to sequence names. The set of G. max PEBP-family protein sequences was subsequently used to screen the genome assemblies of L. angustifolius [29] and eight other legume species using the same settings as those applied for G. max. The transcriptome unigene assembly of L. angustifolius was also searched [17]. Fgenesh + gene predictions were made using G. max protein sequences as references. Predicted PEBP-family protein sequences were compared with annotation data deposited at PeanutBase http://www.peanutbase.org (A. duranensis and A. ipaensis) and LegumeIP database [91] (C. cajan, C. arietinum, L. japonicus, M. truncatula, P. vulgaris). Accession numbers from these databases were assigned to predicted sequences if applicable.

Phylogenetic inference of PEBP-family sequences

The set of 122 PEBP-family protein coding sequences, including 4 reference A. thaliana sequences, was selected (see Additional file 10). Multiple translation sequence alignment was achieved in MAFFT v7.017 [92], using the following parameters: standard genetic code, BLOSUM 62 substitution matrix, gap open penalty 1.25. Two algorithms were tested, G-INS optimized for global homology, and E-INS for multiple conserved domains and long gaps. Bayesian inference of phylogeny was performed in MrBayes 3.2.2 [93]. A phylogenetic tree was drawn in Geneious [88] using specified parameters and settings (see Additional file 11). Based on the topology of the tree, paralogous and orthologous pairs of sequences were selected. Pairwise translation sequence alignments were done in MAFFT v7.017, using E-INS, BLOSUM 80 and gap open penalty 1.25. Selection pressure parameters, Ka (the number of nonsynonymous substitutions per nonsynonymous site), Ks (the number of synonymous substitutions per synonymous site), and Ka/Ks ratios were calculated in DnaSP 5 [94]. The estimation of dates of duplication events was based on the Ks parameter [45].

Identification of conserved motifs in FT promoter sequences