Phylogenetic Placement and Phylogeography of Large-Flowered Lotus Species (Leguminosae) Formerly Classified in Dorycnium: Evidence of Pre-Pleistocene Differentiation of Western and Eastern Intraspecific Groups

The Mediterranean region is a center of species and genetic diversity of many plant groups, which served as a source of recolonization of temperate regions of Eurasia in Holocene. We investigate the evolutionary history of species currently classified in Lotus sect. Bonjeanea in the context of the evolution of the genus Lotus as a whole, using phylogenetic, phylogeographic and dating analyses. Of three species of the section, L. rectus and L. hirsutus have wide Mediterranean distribution while L. strictus has a disjunctive range in Bulgaria, Turkey, Armenia, Eastern Kazakhstan, and adjacent parts of Russia and China. We used entire nuclear ribosomal ITS1-5.8S-ITS2 region (nrITS) and a plastid dataset (rps16 and trnL-F) to reconstruct phylogenetic relationships within Lotus with an extended representation of Bonjeanea group. We analyzed the phylogeographic patterns within each species based on the plastid dataset. For divergence time estimation, the nrITS dataset was analyzed. Our results confirmed the non-monophyletic nature of the section Bonjeanea. They indicate that Lotus is likely to have diverged about 15.87 (9.99–19.81) million years ago (Ma), which is much older than an earlier estimate of ca. 5.54 Ma. Estimated divergence ages within L. strictus, L. rectus, and L. hisrutus (6.1, 4.94, and 4.16 Ma, respectively) well predate the onset of the current type of Mediterranean climate. Our data suggest that relatively ancient geological events and/or climatic changes apparently played roles in early diversification of Lotus and its major clades, as well as in formation of phylogeographic patterns, in at least some species.


Introduction
is the largest genus of the tribe Loteae (Leguminosae-Papilionoideae) containing ca. 130 species of annual and perennial herbs, semishrubs, and shrubs or dwarf shrubs widely distributed in Eurasia, Africa, Australia and several islands of the Atlantic, Indian, and Pacific oceans [1]. Kramina and Sokoloff [2] and Sokoloff [3] divided Lotus into 14 sections, and this classification was used as a base for molecular phylogenetic studies [4][5][6][7]. In the phylogeny of the genus Lotus, whose major center of species diversity is located in the Mediterranean region, an early split into "southern" and Position of sterile bract (see Sokoloff et al., 2007) Often separated from partial inflorescence by a stalk Often separated from partial inflorescence by a stalk Often separated from partial inflorescence by a stalk At the base of partial inflorescence At the base of partial inflorescence At the base of partial inflorescence (Hammatolobium, Tripodion) or absent (Cytisopsis)

Flowers and partial inflorescences
Partial inflorescences more commonly many-flowered (6-30 flowers often arranged on inflorescence axis in more than one whorl) Partial inflorescences commonly few-flowered (a whorl of 1-8 flowers). An important progress has been achieved in dating main evolutionary events in the phylogeny of legumes. The family Leguminosae diversified during the Early Tertiary. First definitely determined legumes appeared in the Late Paleocene, about 56 Ma [23]. The legume stem clade age was estimated at 59.9 Ma [20]. The Papilionoideae stem clade age was set to a minimum of 55 Ma (Late Paleocene) on the basis of fossil records. The robinioid crown clade age, i.е., the MRCA (the most recent common ancestor) of Hebestigma Urb. and Robinia L., was evaluated as 45.4 Ma or 48.3 Ma, according to rbcL and matK based phylogenetic reconstructions, respectively [20]. The Robinia stem clade was directly dated at a minimum of 33.7 Ma based on the fossil wood of Robinia zirkelii (Platen) Matten, Gastaldo & Lee [24][25][26]. The age of the clades Loteae-Sesbania and Loteae was estimated on the basis of previous dating analyses as 36.7 ± 1.8 and 21.5-24.6 Ma, respectively [20,22]. Based on the dating of Lavin et al. [20], Jaén-Molina et al. [7] estimated the time when Lotus separated from its sister group (i.e., Hammatolobium Fenzl and Cytisopsis Jaub. & Spach) at 7.86 Ma, and the time of diversification of extant clades of the genus at 5.54 Ma.
The main center of species diversity of the genus Lotus is located in the Mediterranean region [6]. There is also another important center in Macaronesia, but it is presumably much younger than the Mediterranean one [7]. In recent decades, phylogeography of Mediterranean plants has been extensively studied, most of the studies being focused on sclerophytes and woody plants typical for vegetation of this region (reviewed by References [27][28][29]). The vast species richness of the Mediterranean basin, one of the world biodiversity hotspots [30][31][32], is explained by biogeographic patterns associated with highly heterogeneous landscape, complex geological and climatic history, and long-term human activity [29,33]. The three large Mediterranean peninsulas, i.e., Iberian, Apennine, and Balkan, as well as Anatolia, accumulate much genetic and species diversity, which generally decreases towards higher latitudes. The peninsulas acted as refugia for many species through Pleistocene climatic oscillations [28]. The Mediterranean area was less affected by the latest glaciations [27], and the genetic structure of some Mediterranean species may be the result of older processes (e.g., Reference [33]). Presumably existing from the Tertiary period, Mediterranean refugia are climatically stable areas for long-term conservation of species and genetic diversity [27]. They served as sources for recolonization of central and northern Europe during interglacials [28].
In this connection, the biogeography of the northern evolutionary lineage of Lotus is of particular interest, since each of the sections of the northern lineage includes species common in the Mediterranean, as well as species that have advanced farthest into temperate latitudes. In this paper, we investigate the evolutionary history of three species of Lotus sect. Bonjeanea in the context of the evolution of the genus Lotus as a whole, using phylogenetic, phylogeographic and dating analyses. Two currently recognized species of Lotus sect. Bonjeanea (L. hirsutus L. and L. rectus L.) have wide Mediterranean ranges, and the third species (L. strictus) is spread from east Mediterranean to south-western Siberia and Northern Dzungaria. We also study phylogeography of L. graecus L., formerly included in the section Bonjeanea, as well as related Turkish endemic species L. axilliflorus (Hub.-Mor.) D.D. Sokoloff and L. sanguineus (Vural) D.D. Sokoloff. Sampling of species for the present study is largely dictated by their conflicting positions in earlier analyses of morphological, nuclear and plastid data [4,6,17]. Basically, we deal with all species placed in Dorycnium by Lassen [15] and Greuter et al. [34] except members of the highly polymorphic complex of Lotus dorycnium sensu lato. The L. dorycnium complex includes a few small-flowered species whose limits, diagnostic characters, and phylogeographic patterns will be discussed separately.
Our objectives are: (1) to detail genetic diversity in the species of the Bonjeanea group (i.e., Lotus rectus, L. strictus, L. hirsutus, and L. graecus) across the distribution ranges using molecular data and phylogenetic analyses of the genus Lotus with extended representation of the Bonjeanea group; (2) to assess the phylogeographic patterns and infer the main drivers of differentiation among studied species; and (3) to estimate the ages of major subclades of the genus Lotus using molecular phylogenetic analyses.

Analysis of nrDNA ITS
The ITS1-5.8S-ITS2 region (nrITS) dataset included 107 accessions (Supplementary Materials, Dataset S1), 104 in the ingroup (i.e., Lotus) and three in the outgroups (i.e., Hammatolobium kremerianum, Cytisopsis lotoides and Tripodion tetraphyllum). One hundred percent identical sequences of the same species were combined into one accession. The ingroup covers 31 species of Lotus and 13 of 14 sections of the genus with enlarged representation of Lotus rectus, L. hirsutus, L. strictus and L. graecus. The total alignment length was 667 bp (617 bp after the exclusion of gap-rich and ambiguous positions).
The genus Lotus clade is well supported by Bayesian and maximum likelihood (ML) analyses (posterior probabiliy PP 1.00, bootstrap support BS 99%) and further splits into several branches (Figure 1) Both Bayesian and ML analyses strongly support the separation of L. rectus into Western (PP 1.00, BS 99%) and Eastern (PP 1.00, BS 100%) subclades (Figures 1 and 2A). A similar pattern of clusterization is observed in L. hirsutus. It also splits into Western (PP 1.00, BS 100%) and Eastern (PP 1.00, BS 95%) branches, which are however shorter than the corresponding branches of L. rectus (Figures 1 and 2C). Turkish specimens of L. strictus form a clade (PP 1.00, BS 97%) which is consistently combined first with the Bulgarian, and then with the Altai-Kazakhstan samples (Figures 1 and 2B). Lotus graecus forms a well-supported clade with L. axilliflorus, but L. sanguineus is more distantly related. Bayesian and ML analyses undoubtedly classified three species of L. sect. Canaria (i.e., L. broussonetii Choisy ex Ser., L. spectabilis Choisy ex Ser., and L. eriophthalmus Webb et Berth.) within the Lotus Southern clade and did not support their close relationships with any member of the former genus Dorycnium.

Analysis of Plastid DNA Dataset
The plastid DNA dataset included 107 sequences representing the same outgroup and ingroup as in nrITS analysis plus the sequence of L. hirsutus 03052324 (Supplementary Materials, Dataset S2). The total alignment length of the combined dataset was 1838 bp (incl. rps16 intron of 891 bp and trnL-F of 947 bp), but, after the exclusion of gap-rich and ambiguous positions, the alignment length was reduced to 1703 bp.
The genus Lotus is highly supported by Bayesian phylogenetic analysis (PP 1.00) and slightly weaker by ML analysis (BS 90%) ( Figure 3). Both analyses demonstrate the split of Lotus into the Northern (PP 1.00, BS 97%) and Southern (PP 1.00, BS 94%) clades. The Southern clade combines Lotus species, which belong to nine sections, including L. sect. Canaria and Chamaelotus. The Lotus Northern clade includes several rather well supported branches, i.e., L. dorycnium complex plus L. hirsutus clade, L. graecus plus related species clade, two separate clades and a branch of L. parviflorus of L. sect. Lotus, L. rectus clade and L. strictus clade. Monophyly of none of the sections of the Northern clade (i.e., L. sect. Lotus, Dorycnium, and Bonjeanea) was confirmed by both methods of analysis. The analysis of the plastid dataset revealed a well-supported clade comprising all accessions of L. graecus, L. axilliflorus, and L. sanguineus, but failed to support monophyly of L. graecus. Both L. axilliflorus and L. sanguineus are revealed as well-supported clades on relatively long branches.
The clade of L. rectus includes a basal grade of samples from the western part of its range, an intermediate grade of specimens from Morocco and Mediterranean islands (Crete and Corse) and an Eastern clade. The clade of L. strictus includes a basal grade formed by Turkish and Bulgarian samples and a clade of Altai-Kazakhstan specimens.
Within the clade (L. dorycnium complex plus L. hirsutus), the sequences of L. hirsutus are included in three subclades: some of them form a fairly well supported Western subclade (PP 1.00, BS 92%), the other part is combined in a less supported Eastern 1 subclade (PP 0.76), and the third part (Eastern 2 group) is intermingled with members of L. dorycnium complex in a common subclade.
Similar to nrITS analysis, the shorter branches were observed: in Turkish, rather than in Altai-Kazakhstan, samples of L. strictus; in western, rather than in eastern, samples of L. rectus; in eastern, rather than in some western, samples of L. hirsutus.

Lotus rectus
We analyzed 17 sequences, twelve in the ingroup (i.e., L. rectus) and five in the outgroup, represented by L. strictus, L. hirsutus, L. dorycnium, L. graecus, and L. corniculatus ( Figure 4). The program calculated the parsimony limit of 30 steps and collapsed sequences into 15 haplotypes, five in the outgroup and ten in L. rectus (H1-H10). Haplotype diversity (Hd) and nucleotide diversity (pi) in L. rectus sequences are rather high (0.909 and 0.00436, respectively). On a parsimony network constructed using TCS software ( Figure 4B), the internal haplotype H1 from Spain is the closest to the hypothetical haplotype X, which represents the connection point with sequences of the outgroups (L. strictus, L. graecus, and other species of Lotus). The haplotype H1 differs from X by six mutations. The haplotype H1 is further connected to a branch of haplotypes H2-H4, which differ from it by two, three and four mutations, respectively. The haplotypes H1-H4 are distributed in Spain, France, and north-western Italy and represent a western group of L. rectus ( Figure 4A). A long series of mutations (12 or 13) connects the haplotype H1 with haplotypes H8-H10 from the L. rectus eastern group, which includes samples from Lebanon, Israel, and Turkey. The haplotypes H5, H6, and H7 occupy intermediate positions between western and eastern groups of haplotypes, being closer to the first (H5 from Crete) or second (H6 from Morocco and H7 from Corse). The distribution of pairwise differences between studied sequences is multimodal (Supplementary Figure S1A), which may indicate that L. rectus has been widespread in this area for a long time and has repeatedly expanded and reduced in number. The preliminary morphological analysis of herbarium specimens did not reveal clear morphological differences between western and eastern populations of L. rectus.

Lotus strictus
Of the 22 plastid sequences studied, 17 belonged to L. strictus, and five belonged to the outgroup, represented by L. rectus, L. hirsutus, L. dorycnium, L. graecus, and L. corniculatus ( Figure 5). The program calculated the parsimony limit of 30 steps and collapsed sequences into 10 haplotypes, five in the outgroup and five in L. strictus (H11-H15). L. strictus is characterized by lower gene (Hd = 0.596) and nucleotide diversity (pi = 0.00141) compared to L. rectus. On a parsimony TCS network of L. strictus haplotypes ( Figure 5), the haplotype H11 is the closest to the hypothetical haplotype X, differing from the latter by the only mutation. The haplotype H11 is the most frequent haplotype present in 10 samples from Turkey and Bulgaria. H11 is connected by a short branch of two mutations with the singleton H12 from Bulgaria. The other longer branch connects the haplotype H11 with a group of Altai-Kazakhstan haplotypes H13, H14, and H15, which are consecutively connected to each other and differ from H11 by four, five, and six mutations, respectively. Among this group of haplotypes, the haplotype H14 is more frequent (was found in four specimens). According to the distribution of pairwise differences (Supplementary Figure S1B), the species L. strictus also experienced fluctuations in population size and went through the bottlenecks. Now, it is not currently expanding. Studied herbarium specimens of L. strictus are variable by morphology, but our preliminary morphological analysis did not reveal clear morphological differences between Anatolian-Bulgarian and Altai-Kazakhstan populations.

Lotus hirsutus
We analyzed 37 sequences, including 24 in L. hirsutus, nine in L. dorycnium complex, and four in the outgroup, represented by L. strictus, L. rectus, L. graecus, and L. corniculatus ( Figure 6). The program calculated the parsimony limit of 30 steps and collapsed sequences into 30 haplotypes, four in the outgroup, five in L. dorycnium complex, 20 in L. hirsutus (H16-H27, H29-H36), and one haplotype (H28) shared by L. hirsutus and L. dorycnium s.l. (Figure 6). The haplotype diversity within L. hirsutus (0.960) was the highest among all studied species, while nucleotide diversity (0.00392) was lower than that of L. rectus, but higher than in L. strictus. On a TCS network of haplotypes ( Figure 6B), the hypothetical haplotype Y is a first putative center of divergence of L. hirsutus, which differs from the hypothetical haplotype X by eight mutations. Two main subclades of L. hirsutus (i.e., western subclade and eastern 1 subclade) may originate from the haplotype Y. The eastern 1 subclade includes haplotypes H16-H25 which are distributed from Turkey to Italy. This group of haplotypes starts to diverge from the haplotype H16 from Turkey which differs by the only mutation from the hypothetical haplotype Y. The eastern 1 sublclade of L. hirsutus is the most variable group by Hd (0.927) within the species. The western subclade includes haplotypes H30-H36 from the Iberian peninsula. It presumably starts to diverge from the hypothetical haplotype N, differing from the haplotype Y by four mutations. This group includes mainly singletons (except for H35), which are either internal, or tip by position. Two other branches are connected to the point Y, these are the branch towards the third group of haplotypes, L. hirsutus eastern 2 group, and the branch of haplotypes of L. dorycnium complex. The eastern 2 group of L. hirsutus haplotypes, like several haplotypes of the L. dorycnium complex, possibly originates from the point Z, a hypothetical haplotype which differs from the haplotype Y by two mutations. The closest haplotype to Z is haplotype H28 shared by two Turkish specimens of L. hirsutus (samples 7 and D11), as well as three specimens of L. germanicus (samples D1, D4, and D5), a member of L. dorycnium complex. Samples D1, D4, and D5 are geographically located in S. Germany and in the Balkans. Three other haplotypes (H26, H27, and H29) from the Balkans and France also belong to the eastern 2 groups of L. hirsutus. They are connected to either Z point or the haplotype H28, differing by several mutations.
To summarize, the haplotypes of L. hirsutus can be divided into three groups, i.e., western, eastern 1 and eastern 2, and the eastern 2 group of haplotypes is combined with haplotypes of L. dorycnium complex, moreover, there is one haplotype (H28) shared by two taxa. The distribution of pairwise differences in a plastid set of L. hirsutus is multimodal, which indicates a long-term existence of the species in its distribution area and occurrence of demographic fluctuations (Supplementary Figure S1C).

Lotus graecus and Related Taxa
We studied 27 plastid sequences, including 17 sequences of L. graecus, three sequences of L. axilliflorus, two sequences of L. sanguineus and five sequences of the outgroup, represented by L. rectus, L. hirsutus, L. strictus, L. dorycnium, and L. corniculatus. The program calculated the parsimony limit of 30 steps and collapsed sequences into 16 haplotypes, eight in L. graecus, two in L. axilliflorus, one in L. sanguineus, and five in the outgroup (Figure 7). The haplotype diversity of L. graecus sequences is comparatively low and equal to that of L. strictus (Hd = 0.596), and its nucleotide diversity (0.00055) is the lowest among studied species. The central haplotype H37, widely distributed in Turkey, Greece and the Caucasus, is the most frequent (observed in 9 samples) (Figure 7). Other haplotypes (H38-H44) differ from it by one or two mutations and have narrow geographic distribution. They are mainly singletons except for H40, which was recorded in two samples from Turkey and the Caucasus. Two haplotypes derived from H37 were identified in two Crimean samples. The distribution of pairwise differences is unimodal with maximum at 0, which may indicate that this species has recently undergone demographic expansion, while low values of diversity coefficients may be evidence of its relatively small age (Supplementary Figure S1D). The haplotypes of both Turkish endemic species (L. sanguineus and L. axilliflorus) seem to originate from the haplotype of L. graecus. Both samples of L. sanguineus belong to the same haplotype H45, which is related to the central haplotype of L. graecus, differing from the latter by four mutations. L. axilliflorus is characterized by two haplotypes, the internal haplotype H46 and the tip haplotype H47, which differ from the central haplotype of L. graecus by four and six mutations, respectively.

Dating Phylogeny of Lotus
Results of Bayesian dating analyses with different speciation priors were very similar, so we present and discuss only results obtained using the Birth-Death branching pattern ( Figure 8). The analyses revealed a medium level of rate heterogeneity with the mean of the coefficient of rate variation 0.53 (95% highest posterior density HPD interval 0.3744-0.6978) and mean clock rate 4.22 × 10 −9 substitutions/site/year. The mean age of the Loteae-Sesbania clade is estimated from nrITS sequences at  Ma), whereas the Loteae crown age is 29.11 (18.28-34.82) Ma ( Figure 8, Table 2). Our results indicate that Lotus likely diverged about 15.87 (9.99-19.81) Ma. Within Lotus, the mean ages of large clades (i.e., Lotus Northern clades 1 and 2 and Lotus Southern clade) varied within 10.21-12.47 Ma, while species of Lotus section Bonjeanea, L. strictus, L. rectus, and L. hisrutus, likely start to diverge at mean ages estimated as 6.1, 4.94, and 4.16 Ma, respectively.  Table 2.

Diversification of the Genus Lotus Much Pre-Dates Formation of the Extant Mediterranean Biome
Our data provide an updated framework for understanding tempo of evolution of the genus Lotus. Recently, Jaén-Molina et al. [7] published a dated phylogeny of Lotus based on nrITS sequences of 116 species, which is about 94% of the total species diversity in the genus [4], plus nine species from five closely related genera within Loteae and two species of Sesbania as outgroups. The analysis of Jaén-Molina et al. [7] was focused on colonization of Macaronesian islands, representing a collection of relatively recent diversification and dispersal events mostly confined to a particular clade (the Pedrosia clade). Our dating analysis is also based on nrITS sequences, but it is focused on several lineages that are closer to the root of Lotus than the Pedrosia clade. Therefore, our taxon sampling scheme involved a much broader set of outgroups that included most currently recognized genera of the tribe and major clades of Sesbania and Robinieae. While Jaén-Molina et al. [7] used an estimate age of Loteae inferred from an analysis of a matK dataset [20], we found it reasonable to use the robinioid legumes age as a secondary calibration point, which is less-dependent on the analyzed dataset [20]. Jaén-Molina et al. [7] estimated an age of the crown group of Lotus at an average of 5.54 (3.45-7.90) Ma (i.e., close to the Pliocene to Miocene boundary), but our analysis revealed an about three times older the Miocene estimate. Consequently, ages of other clades within the genus are also older than those in the study of Jaén-Molina et al. [7].
As the genus Lotus and the tribe Loteae have major diversity centers in the Mediterranean, dated phylogenies have to be discussed in the context of important events in climate and geology of the Mediterranean region. The most important series of events of the Late Miocene is the Messinian salinity crisis (MSC) (5.96 to 5.33 Ma), when the Mediterranean Sea had no stable connection to the Atlantic Ocean and exhibited a great desiccation. This period ended with the appearance of the Gibraltar strait (the Zanclean flood) [35]. The Messinian salinity crisis played important roles in plant evolution and distribution patterns [28]. The influence of MSC on the formation of the diversity of Mediterranean plant lineages has been extensively discussed. In the study of Stauracanthus Link (Leguminosae), Pardo et al. [36] documented the negative effects of MSC causing plant retractions and local extinctions. The studies of Jabbour and Renner [37] on Consolida Gray (Ranunculaceae) and Lledó et al. [38] on Limonium Mill. (Plumbaginaceae), on the opposite, demonstrated the role of MSC as a factor facilitating geographic expansion and promoting radiations. The data of Jaén-Molina et al. [7] suggest that major contemporary lineages of the genus (including the Northern and Southern groups [6]) appeared close to the Messinian/Zanclean border. In contrast, our analysis suggests that both Northern and Southern groups were well diversified by the beginning of the Messinian crisis. Moreover, our inferred mean crown group age of apparently the most enigmatic species of the genus, L. strictus, pre-dates the Messinian crisis and our estimates for two other members of the section Bonjeanea, L. hirsutus and L. rectus, are close to the time of the crisis, at least with respect to their stem groups. Importantly, our inferred crown clade ages of all three members of the section Bonjeanea exceed existing estimates of the age of the contemporary Mediterranean summer-drought regime (c. 2.8 Ma [39]).
Phylogeographic differentiation between western and eastern parts of distribution ranges is documented in several Mediterranean [29,[40][41][42][43] and submiditerranean [44] angiosperm species. Apparently, such differentiation did not appear synchronously in various taxa. Our data suggest that differentiation between western and eastern lineages of L. hirsutus and L. rectus most likely took place as early as in Pliocene. In contrast, a similar differentiation in another member of the tribe Loteae, the submediterranean species Anthyllis montana L. has been dated as Late Quaternary [44]. Despite the apparently much younger age of geographical separation, morphological differentiation between the western and eastern groups of A. montana [45] is much more pronounced than between the groups within L. rectus and L. hirsutus. The very limited incongruence between molecular and morphological patterns of differentiation in A. montana somewhat resembles the incongruence between the plastid and molecular marker we observed in L. hirsutus and suggest complex evolutionary histories of these lineages. Kropf et al. [44] revealed a western Mediterranean origin followed by an eastward migration in A. montana, which is similar to the direction we hypothesize for L. rectus, though with much younger dating.
An example of molecular clock estimate suggesting an ancient Western/Eastern Mediterranean differentiation is provided by Casimiro-Soriguer [43] for the monospecific herbaceous legume Erophaca Boiss. (Leguminosae). Casimiro-Soriguer [43] concluded that Erophaca is one of the many Tertiary relicts that form part of the present Mediterranean flora. The most recent common ancestor of the eastern and western Mediterranean subspecies of Erophaca baetica Boiss. was dated using relaxed molecular clock as 11.9 Ma, although with a large 95% confidence interval, 3.3-19.1 Ma, from the Late Miocene to the Early Pliocene. This estimate approaches our estimates of differentiation in L. rectus and L. hirsutus. Similar temporal patterns of infraspecific differentiation was revealed in the Mediterranean palm Chamaerops humilis L. [46]. As revealed by Hardion et al. [47], the origin of the Mediterranean thorny cushion-like xerophytes from Astragalus L. sect. Tragacantha DC. (Leguminosae) takes root in the middle of the Pliocene (4.4 Ma), between the Messinian salinity crisis and the onset of the Mediterranean climate.
Ackerly [48] suggested that the distinction between the age of traits and the age of taxa should be made while analyzing historical biogeorgaphy of Mediterranean floras. In this respect, it is puzzling that it is not easy to formulate which morphological traits are associated with diversification of Lotus in the Mediterranean region. For example, plants of L. strictus s.l. from non-Mediterranean part of the range (described as L. strictus Fisch. et C.A.Mey. s. str. from Armenia and L. albus Janka from Bulgaria), are morphologically similar to those from the Mediterranean part (described as L. thermalis Boiss.). At least, it is unlikely that characters once used to segregate L. albus (e.g., white corolla color) can be used as markers of any climatic adaptation. Apparently, there is only one species of Lotus that clearly fits the criteria of a Mediterranean habit, namely L. fulgurans. This is a dwarf thorny shrub from the Balearic islands. The climatic-based nature of its habit is further supported by the occurrence of an externally similar species of the same tribe, Anthyllis hystrix (Willk. ex Barceló) Cardona, Contandr. & Sierra that occurs on the same archipelago [49,50]. Our inferred divergence time estimate of L. fuglugans and L. germanicus is well below 2.8 Ma. Thus the estimated age of L. fulgurans agrees with the idea that its characteristic habit evolved as an adaptation to the current Mediterranean type of climate.
Our sampling of L. strictus covers Anatolian, Bulgarian and Altai-Kazakhstan parts of its distribution range. Transcaucasian and Dzungarian parts of the range were not sampled. The Transcaucasian (Armenian) localities lie between the Anatolian and Altai-Kazakhstan parts, so we may suppose that they presumably have plastid haplotypes intermediate between the two groups of haplotypes revealed in Anatolian-Bulgarian and Altai-Kazakhstan groups. This hypothesis needs testing using plant material form the Transcaucasian part of L. strictus range. The Dzungarian part of the distribution range is located very close to the Altai-Kazakhstan part, so we may suppose that Dzungarian plants of L. strictus may possess plastid haplotypes of the same group. This assumption also needs verification using additional sampling from Northern Xinjiang.
The results obtained in the present study suggest that the distribution range of L. strictus was large and then experienced reduction and fragmentation. These processes led to population size decline and extinction of the haplotypes intermediate between two haplogroups. Our data allowed to suppose that the ancestral haplotypes of L. strictus putatively originated from the Anatolian part of the range and then the species spread to Transcaucasian and Altai-Kazakhstan-Dzungarian parts.
The extant range of L. strictus is intriguing. There is an impressive gap of almost 3000 km between the localities in Armenia (that closely approach those in eastern Turkey [58]) and closest localities of the species in the Eastern part of Kazakhstan (to the East of Irtysh valley in the North and Lake Alakol in the south [59]). It is unlikely that appropriate habitats are lacking within this gap, as well as in NW Kazahkstan and SW European Russia, where saline lakes are rather common. Based on our dated phylogeny, one may speculate that dispersal of L. strictus to the East took place in Late Miocene along the southern shore of Eastern Paratethys, a marine basin that still existed by that time (e.g., Reference [60]) and apparently prevented a direct migration northwards. The extant range of L. strictus may be related to margins of the Eastern Paratethys.
Lotus rectus is a perennial that grows along the margins of water courses, damp and bushy places, in preferably basic substrates at altitudes from 0 to 1300 m [61,62]. The species has a wide Mediterranean distribution. It occurs in Portugal and in all countries around the Mediterranean Sea, except for Egypt and the countries of the former Yugoslavia (Slovenia, Croatia, Montenegro) [34]. Our sampling rather well covers the European part of L. rectus distribution range (except for Portugal, Sardinia, Sicily, and Greece) and its Asian part. The African part of the range is insufficiently represented in our study and includes a sample from Morocco. Our results imply the east-west phylogeographic differentiation of L. rectus, which presumably started at 4.94 Ma, soon after the Messinian crisis. The nature of genetic variability suggests that L. rectus has repeatedly experienced fluctuations in abundance. The present phylogeographic pattern may be a result of the range reduction, which led to decrease of gene flow between its western and eastern parts. The western and eastern populations are apparently not completely genetically isolated, as evidenced by the Moroccan and two island samples, which demonstrate contrasting genetic relatedness by nrDNA ITS and plastid data. Such a longitudinal pattern corresponds to phylogeographical or even geographical break described for many Mediterranean plant species (e.g., see review in Reference [29]). The present analysis of plastid dataset of L. rectus allows to suppose the location of the ancestral area of this species in the Western part of its range (i.e., in the Iberian peninsula) and its further distribution from the Western to Eastern Mediterranean.
Lotus hirsutus is a sub-shrub or small shrub growing on calcareous slopes and cliffs, dry hills, macchie, pine forests at altitudes from 0 to 900 m. This species, like L. rectus, has an extensive Mediterranean range. It occurs in Portugal and in all countries around the Mediterranean Sea, except for Egypt, Tunisia and Morocco [34]. Our sampling covers well the European part of its distribution range, but only samples from Anatolia were available for the Asian part. Samples from Africa (Libya and Algeria) and large islands (Corse and Sardinia) were not available. According to our dated phylogeny, based on ITS, separation between L. dorycnium complex and L. hirsutus may happen at about 5.55 Ma and the divergence of L. hirsutus into western and eastern geographical groups could occur about 4.16 Ma. Genetic variation pattern of L. hirsutus plastid dataset implies that diversification of the species likely started in the eastern part of the range, in Anatolia, from where the species spread to the western part (the Iberian peninsula), where it also experienced a later divergence. The third plastid evolutionary branch of L. hirsutus is scattered within the eastern part of the range (from Turkey to France) and demonstrates incomplete lineage sorting with L. dorycnium complex. It can be assumed that L. hirsutus, like L. rectus, exhibits longitudinal phylogeographic differentiation, but with the opposite direction of distribution from east to west.
Lotus graecus is a perennial herb with root suckers [17], growing on roadsides, open slopes, macchie, coniferous and deciduous forests at altitudes from 0 to 2000 m [52,61,63]. It is distributed in Greece, Bulgaria, Turkey, the Caucasus, and the Crimea. The main part of the range is located around the Black Sea reaching the East Mediterranean in the South. Our samples of L. graecus sufficiently well cover its entire range. The phylogeographic analysis of the L. graecus plastid dataset suggested a relatively young age of the species, which has recently undergone demographic expansion, as evidenced by the presence of one widely distributed central haplotype and several derived haplotypes. Both ITS and plastid data clearly indicate that two Turkish endemics, L. axilliflorus and L. sanguineus, are related to L. graecus. Lotus axilliflorus is a perennial occuring in oak scrub on marly soil in SW Anatolia (Burdur prov.) [61]. Lotus sanguineus is a perennial herb with a very restricted distribution in Southern part of Central Anatolia, near Bucakkışla in Karaman. It grows on south-facing slopes on calcareous substrate, in open Pinus brutia forests and maquis, at altitudes from 350 to 900 m [64]. According to the results on ITS dating phylogeny obtained in the present study, the separation of L. sanguineus from L. graecus and L. axilliflorus can be estimated at about 6.38 (2.715, 9.641) Ma, before the Messinian crisis, and two latter species diverged from each other much later, at about 2.81 (0.776, 5.042) Ma.

Plant Material
The plant material used for the present study was sampled from herbarium specimens stored in several large Herbaria (ALTB, ANK, GAZI, LE, MA, MHA, MW, NS, P, SO, SOM, WAG, ZA). In total, we analyzed 104 specimens, which belong to 32 species and thirteen sections of the genus Lotus with expanded representation of species of Bonjeanea group (i.e., L. hirsutus, L. rectus, L. strictus, and L. graecus). Cytisopsis pseudocytisus (Boiss.) Fertig, Hammatolobium kremerianum (Coss.) Müll. Berol., and Tripodion tetraphyllum (L.) Fourr. were used as an outgroup. For 70 specimens, the sequences of all the studied DNA regions or some regions were obtained in this study, the rest of the sequences were taken from GenBank. Voucher information and GenBank accession numbers are presented in Appendix A. Distribution maps of the studied specimens were constructed using SimpleMappr [65].

DNA Extraction, Amplification and Sequencing
DNA was extracted from dry leaves taken from herbarium (ca. 20 mg of leaf tissue) with NucleoSpin Plant II kit (Macherey-Nagel, Germany) according to the manufacturer's instructions or using the CTAB (cetyl trimethylammonium bromide) method [66].
The sequences of the nrITS were amplified with primers NNC-18S10 иC26A [67] and universal primers ITS2 and ITS3 [68]. The sequences of trnL-trnF intergenic spacer (IGS) and trnL intron of plastid DNA were amplified using standard primers 'c', 'd', 'e' and 'f' [69]. The sequences of rps16 intron of plastid DNA were amplified using primers rpsF and rpsR2 [70]. We also used forward internal primer rps16internF [71] in a combination with rpsR2 for amplification of the second (3 ) part of rps16 intron. However, the internal reverse primer rps16internR [71] in a combination with rpsF gave no amplification in Lotus samples. To amplify the first (5 ) part of rps16 intron, we designed a pair of primers specific to Lotus, and Lot-rps16-intR . The primers were developed using the Primer-BLAST [72] software based on the complete sequence of Lotus japonicus (Regel) K.Larsen plastome (GenBank accession number NC_002694.1). Newly developed primer Lot-rps16-intR in a combination with either rpsF or Lot-rps16-F produced good results in many Lotus species and allowed to amplify sequences about 500 bp long.
PCRs PCR products were purified using Cleanup Mini kit (Evrogen, Moscow, Russia) following the manufacturer's instructions. Direct sequencing was performed on the ABI PRISM 3100 genetic analyzer (Applied Biosystems, Foster City, CA, USA), using ABI Prism BigDye Terminator Cycle Sequencing Ready Reaction Kit v. 3.1 for cycle sequencing reactions following the manufacturer's instructions. Forward and reverse strands of all samples were sequenced. The polymorphism of ITS within one specimen was detected by direct sequencing (without cloning), by the presence of double peaks on electropherogram.
The sequences were aligned using MAFFT version 7.215 [73,74] and then adjusted manually in BioEdit version. 7.2.5 [75]. The matrices of rps16 intron and trnL-F cpDNA regions were combined into a single matrix. Gap-rich and ambiguous positions were excluded from the analyses. The aligned data matrices are presented in on-line Supplement (Datasets S1-S2).

Phylogenetic Analyses
The Maximum Likelihood analyses were performed with MEGA X [76] with GTR+Г+I model of nucleotide substitutions (the general time-reversible model with the presence of invariable sites and substitution rate heterogeneity) for the plastid dataset sequences and GTR+Г model for nrITS sequences. The models were determined as the best choice for corresponding datasets following the Model Selection option implemented in MEGA X based on the corrected Akaike information criterion (AICc). Nonparametric bootstrap method with 500 replications was used for branch support assessment.
The Bayesian inference was performed using MrBayes v. 3.2.6 [77] considering the optimal model of nucleotide substitutions selected by AICc in PAUP version 4.0a [78] for each marker: SYM+Г (symmetrical model with substitution rate heterogeneity) for nrITS, and GTR+Г for plastid data. The Bayesian analysis used two independent runs of 20 million generations and four chains sampling every 1000th generation. Non-convergence assessment and burn-in estimation was carried out in VMCMC ver. 1.0.1 [79]. The first two million generations were discarded as burn-in and the remaining trees from both runs were combined in a 50% majority-rule consensus tree.
Phylogenetic relationships among the cpDNA haplotypes were reconstructed using statistical parsimony analysis as implemented in TCS v1.2 [80]. Long indels were reduced to one character, then gaps were treated as fifth state. TCS networks of haplotypes were constructed separately for L. rectus, L. strictus, L. hirsutus, and L. graecus. L. corniculatus and L. dorycnium were used as outgroups. Haplotype networks were then visualized using the online program tcsBU [81]. Parameters of genetic variability were calculated using DnaSP 6 software [82].

Dating Analyses
For dating the Lotus phylogeny, we used partially reduced dataset of Lotus nrITS sequences (45 accessions representing 30 species) but with an enlarged outgroup, including 13 genera of the tribe Loteae, two species of Sesbania, three species of Robinia and Coursetia glandulosa A.Gray from Robineae (nrITS sequences of the outgroup were retrieved from GenBank; accession numbers are presented in Appendix B). The aligned data matrix is presented in on-line Supplement (Supplementary Materials, Dataset S3). For divergence time estimation, the nrITS dataset was analyzed using BEAST version v.2.6.3 [83]. The best performing substitution model (SYM + Г) was implemented along with the Yule [84] and Birth-Death [85] speciation priors. Two clock models (strict clock and uncorrelated relaxed clock with a log-normal distribution) were tested using the nested sampling approach [86] implemented in the NS package of BEAST 2, and the strict clock model was rejected basing on marginal likelihood estimations (Table 3). We used two calibration points in phylogenetic dating: a normal distribution (mean = 48.3, SD = 1) was applied to the root age according to the Robinioid crown age estimated by Lavin et al. [20], and a log-normal distribution (mean = 5, S = 0.7) with a minimum age constraint (offset) at 34 Ma to the Robinia stem age according to Robinia L. wood fossil (discussed in [20,26]). Two Markov chains Monte Carlo were run for 40 million generations, trees and parameters were sampled every 4000 generations. Tracer ver. 1.7.1 [87] was used to assess the chain convergence, burn-in and the effective sample sizes. The effective sample size exceeded 200 in all analyses, and burn-in was set to 10%. Before conducting the dating analyses, sampling from prior distributions only was performed to be sure that the marginal distribution of the priors reflected our intended settings.

Conclusions
The results of the present phylogenetic study of Lotus section Bonjeanea, like previously obtained data (Degtjareva et al., 2006(Degtjareva et al., , 2008Kramina et al., 2016), clearly demonstrated the non-monophyletic nature of the section. Our data suggest that L. rectus and L. strictus differentiated at the early stages of evolution of the genus Lotus and represent two separate evolutionary lineages, which can probably be considered as two distinct sections. Our preliminary morphological analysis of herbarium specimens did not reveal clear morphological differences between western and eastern populations of L. rectus, as well as between Anatolian-Bulgarian and Altai-Kazakhstan populations of L. strictus, despite significant genetic differences. Currently, the differentiation observed within each of these two species is probably consistent with the concept of cryptic species; however, both species require more careful morphological analysis. The present study further supports our earlier-quite unexpected-findings (Kramina et al., 2016) on close relationships between L. hirsutus and members of the L. dorycnium complex with strong incongruence between plastid and nuclear data.
Author Contributions: T.E.K. contributed to the field collection of samples, all the laboratory work including DNA extraction, PCR, sequencing, analyses of the data, figure preparation, manuscript writing and editing. M.V.L. contributed to the field collection of samples, all the laboratory work including DNA extraction, PCR, and sequencing. T.H.S. contributed to phylogenetic analyses of the data and preparation of the manuscript, I.A.S. contributed to review and editing manuscript, M.U.Ö. contributed to material collection from herbaria, and D.D.S. contributed to manuscript writing and editing. All authors have read and agreed to the published version of the manuscript. Biodiversity and Ecosystem Research, Bulgarian Academy of Sciences) for providing material for the study, and two anonymous reviewers for their helpful comments.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
Taxa, sample code, locality, voucher information (herbarium code) for the material included in the phylogenetic analyses. GenBank accession numbers are given for the three markers sequenced, rps16 intron, trnL-F, ITS (new sequences indicated by an asterisk); Herbarium codes according to Index Herbariorum. An n-dash denotes a missing marker (GenBank accession numbers of new sequences will be inserted in the final version of the article).