When Morphology and Biogeography Approximate Nuclear ITS but Conflict with Plastid Phylogeny: Phylogeography of the Lotus dorycnium Species Complex (Leguminosae)

Lotus dorycnium s.l. is a complex of taxa traditionally regarded as members of Dorycnium. It has a wide Mediterranean range, extending in the north to Central and Eastern Europe, and in the east to the Crimea, the Caucasus, and the Western Caspian region. Molecular phylogenetic data support placement of the L. dorycnium complex in the genus Lotus. The present study investigated the phylogeny, phylogeography and morphological variability of the L. dorycnium complex across its distribution range to reveal the main trends in genetic and morphological differentiation in this group. The results of the morphological analyses demonstrated some degree of differentiation, with L. d. ssp. herbaceus, ssp. gracilis, and ssp. anatolicus more or less well defined, whereas ssp. dorycnium, ssp. germanicus, and ssp. haussknechtii can be hardly distinguished from each other using morphology. Analyses of the L. dorycnium complex based on nrITS revealed a tendency towards a geographic differentiation into Western, Eastern, and Turkish groups. Phylogenetic and phylogeographic analyses of the same set of specimens using concatenated plastid markers trnL-F, rps16, and psbA-trnH demonstrated a low resolution between the L. dorycnium complex and L. hirsutus, as well as among the taxa within the L. dorycnium complex, which can be interpreted as evidence of an incomplete lineage sorting or hybridization. The evolutionary processes responsible for incongruence in phylogenetic signals between plastid and nuclear sequences of the morphologically well-defined species L. dorycnium and L. hirsutus were most likely localized in the Eastern Mediterranean. A possibility of rare gene exchange between the L. dorycnium complex and the group of L. graecus is revealed for the first time.


Introduction
Lotus is the largest and most taxonomically complicated genus of the tribe Loteae (Papilionoideae-Leguminosae). Dorycnium Mill. was traditionally accepted as a distinct genus by European botanists [1,2], but on the global scale it cannot be properly separated from Lotus in terms of morphology [3][4][5][6]. Molecular phylogenetic data showed that even at the European scale separation of Dorycnium is strongly problematic [7,8]. In a monograph of Dorycnium which was published 120 years ago but which still remains the latest detailed worldwide study of the group, Rikli [9] recognized three sections within the genus: Canaria Rikli, Bonjeanea Taubert, and Eudorycnium Boiss. (the valid name of the latter section is Dorycnium). Members of the sections Canaria and Bonjeanea combine morphological  (Figure 1). Distribution range: East Mediterranean, Balkan Peninsula, extending westwards to Italy and southeastern France, northwards to Germany, Poland, Czech Republic, Slovakia, and Transcarpathian Ukraine, eastwards to the Northern part of Asia Minor, the Crimea, the Caucasus, and Transcaucasia, and to the Western part of the Caspian region [11][12][13]. Dorycnium intermedium has been described from the Crimea as a species close to D. herbaceum, but differing from the latter mainly in a patent (not appressed) pubescence on the calyx [14]. Rikli [9] believed that D. intermedium and D. herbaceum do not differ either morphologically or geographically, but Steinberg [15] considered D. intermedium the easternmost race of D. herbaceum.

2.
Lotus dorycnium ssp. gracilis (Jord.) Kramina Figure 6). Restricted to Balearic Islands. According to morphological [17,18] and molecular data [7], it is a member of the L. dorycnium complex; we have therefore included it in the study.      Diagnostic morphological characters of the taxa listed above are presented in Table 1. Greuter et al. [11] also included the Turkish endemics Dorycnium amani Zohary and D. axilliflorum Huber-Morath in their concept of the group of D. pentaphyllum (=L. dorycnium s.l.). Dorycnium amani was not included in the present study (or any other molecular-based work), because it is known from the type collection only and we had no material for this species. Besides, such a character of D. amani as a pronounced leaf rachis makes the inclusion of this species into the L. dorycnium complex debatable. The previous studies of D. axilliflorum demonstrated that it belongs to the L. graecus L. species group [7,8], rather than to the L. dorycnium complex. Table 1. Morphological characteristics of subspecies of Lotus dorycnium s.l. [1,9,15,17,19,20]. As the basic taxa recognized here within the L. dorycnium complex are subspecies, one should not expect that every single specimen can be identified up to that level. Indeed, some specimens included in the present study shared morphological and/or molecular features with two subspecies. For simplicity, these are tentatively called here hybrids, though in some instances more data should be accumulated to demonstrate clear evidence of reticulate evolution. Previous studies demonstrated that the L. dorycnium complex is not genetically isolated from L. hirsutus of Lotus sect. Bonjeanea ( Figure 7) using plastid data, but monophyletic according to analyses of nuclear ribosomal markers ITS1-5.8S-ITS2 (nrITS) and 5 ETS, though sampling was relatively low [7,8,21]. Based on the nrITS data set, the L. dorycnium complex is younger than L. strictus, L. rectus, and L. hirsutus (estimated divergence times 2.52 Ma, 6.1 Ma, 4.94 Ma, and 4.16 Ma, respectively) [8]. Geographical distribution of the Lotus dorycnium complex and its subspecies based on specimens included in our molecular analyses is presented in Figure 8.
The main aims of the present study were: (1) to investigate genetic diversity and differentiation within the L. dorycnium complex across the whole distribution range using nrITS and a set of plastid DNA markers; (2) to compare degrees of genetic and morphological differentiation in this group; (3) to reveal geographic trends in genetic variability; (4) to refine data on interrelationships between the L. dorycnium complex and L. hirsutus using a more representative sampling of specimens and an enlarged set of plastid markers; (5) to test the hypothesis that nuclear markers allow precise separation of L. hirstus and the L. dorycnium complex.

Morphometric Analyses
Analysis of 89 individuals of the L. dorycnium complex using two characters (flower length and number of flowers per umbel) revealed a division of the dataset into two groups ( Figure 9A):

1.
Blue group: umbels with numerous (usually more than 13) small (usually shorter than 4.5 mm) flowers; this group includes ssp. herbaceus and ssp. gracilis.
germanicus took the position between the two groups. The type specimen of Dorycnium intermedium included in the study for a comparison with other samples of L. d. ssp. herbaceus was placed in the center of the L. d. ssp. herbaceus cluster. The linear correlation coefficient between the variables "flower length" and "number of flowers per umbel" is −0.685 (p = 1.3204 × 10 −12 ).
In DA, we set the group numbers for 'pure' subspecies only: ssp. anatolicus, ssp. haussknechtii, ssp. dorycnium, ssp. germanicus, ssp. herbaceus, and ssp. gracilis. The specimens of presumed hybrid origin were left without a group number. The analysis revealed three main clusters in the dataset ( Figure 9B): 1. ssp. herbaceus; 2. ssp. gracilis; 3. ssp. dorycnium, ssp. germanicus, ssp. anatolicus, and ssp. haussknechtii. In the third group, ssp. anatolicus differs slightly from others along the second axis. The specimens L. d. ssp. anatolicus × L. d. ssp. haussknechtii took the position between corresponding subspecies. The positions of other hybrid specimens and the type specimen of D. intermedium were similar to those obtained on a 2D scatterplot. Note that we were unable to sample type material of Dorycnium herbaceum or any other specimen from the region from which it was described (southeastern France).
All three methods of analysis demonstrated the absence of clear morphological differentiation within the group that includes L. dorycnium ssp. dorycnium, L. d. ssp. haussknechtii, and L. d. ssp. germanicus.

Phylogenetic Analysis of nrDNA ITS1-2 Dataset
The nrITS dataset included 125 accessions (Supplementary Materials, Dataset S1), 64 accessions from the Lotus dorycnium complex, 58 from other Lotus species, and three accessions from outgroups represented by Hammatolobium kremerianum, Cytisopsis pseudocytisus, and Tripodion tetraphyllum. For each of the samples L7 of L. corniculatus, 425 of L. cytisoides, and GRAC2 of L. dorycnium ssp. gracilis (Table 2), two different nrITS sequences were obtained through direct sequencing (without cloning). The total alignment length was 672 bp (601 bp after the exclusion of gap-rich and ambiguous positions). From 601 sites, 252 were variable and 175 parsimony informative. Table 2. Taxa, sample code, voucher information, and GenBank accession numbers of Lotus dorycnium complex specimens used in molecular and morphological analyses. Herbarium codes according to Index Herbariorum. New sequences indicated by an asterisk (GenBank accession numbers for rps16 will be added to the final version of the article.).    The Bayesian phylogenetic tree topology is presented in Figure 10 where posterior probabilities (PP) of nods are given together with bootstrap support (BS) obtained in a maximum likelihood (ML) analysis constructed using IQ-tree software. The results obtained in the ML analysis conducted using RAxML software generally correlate with those obtained in IQ-tree. ML trees based on the nrITS dataset are presented in Supplementary Materials (Figures S1 and S3). Both the Bayesian and ML analyses revealed a highly supported monophyly of the genus Lotus. All specimens of the Lotus dorycnium complex were clustered within a highly supported clade (PP = 1.00/BS = 100%), which was sister to a less supported L. hirsutus clade. The latter was subdivided into Eastern and Western subclades. Lotus graecus and related species (L. axilliflorus and L. sanguineus) were more distantly related to the clade of L. dorycnium and L. hirsutus.
Within the Lotus dorycnium complex clade, the following subclades can be distin-   Table 2 and Appendix A for voucher information.
The Bayesian phylogenetic tree constructed by the plastid DNA dataset is presented in Figure 11 and supplied by PP values and BS values of nods obtained in ML analysis conducted using IQ-tree, as in the analysis of the ITS dataset. ML trees based on the plastid dataset are presented in Supplementary Materials (Figures S2 and S4). The monophyly of the genus Lotus was supported by both Bayesian and ML methods (PP = 1.00/BS = 99%). Two main clades within Lotus, namely Lotus Northern clade and Lotus Southern clade, were well confirmed (1.00/99% and 1.00/97%, respectively). All samples of the Lotus dorycnium complex were placed in the Lotus Northern Clade; however, they did not form a separate subclade, but were combined with L. hirsutus within a highly supported common clade (1.00/100%). Within this common clade, examined specimens of L. hirsutus formed two subclades (Western subclade and Eastern 1 subclade) and several specimens were scattered among the specimens of the L. dorycnium complex (Eastern 2 group). None of the clades within the L. dorycnium complex revealed in ITS analyses were observed in analyses based on plastid markers. Only two subclades more or less corresponding to the clades in the ITS analyses were found. These were (1) a clade combining accessions of the ssp. anatolicus and ssp. haussknechtii and (2) the L. dorycnium Western clade, but their composition differed slightly from those obtained by ITS data. In accordance with its distribution in the Balearic Islands, Lotus d. ssp. fulgurans was revealed as a member of the L. dorycnium Western clade. The most surprising result from the analysis of the plastid DNA is the very well supported placement of two samples of L. d. ssp. haussknechtii (4980 and 9391) within the clade that includes L. graecus plus related taxa.

Haplotype Network Based on the Concatenated Plastid DNA Dataset
TCS analysis included 85 sequences of concatenated plastid DNA regions trnL-F, rps16 intron, and psbA-trnH, 60 sequences of the L. dorycnium complex, 21 sequences of L. hirsutus, and four sequences of the outgroup (L. corniculatus, L. rectus, L. strictus, L. graecus). Sixty-seven haplotypes have been revealed in the dataset, four in the outgroup, and 63 in L. hirsutus and L. dorycnium s.l. (Figure 12). The majority of the haplotypes, 58 of 67 (or 86,6%), are singletons. We revealed three haplotypes shared between subspecies of L. dorycnium and one haplotype shared between L. dorycnium and L. hirsutus. Three of the shared haplotypes were inner, and only one haplotype, shared between HERB4 and PENT02 specimens, was a tip haplotype.
To study genetic diversity within the L. dorycnium complex and L. hirsutus, we divided all specimens of these taxa according to their geographical origin ( Figure 8, [8]) and placement in the clades of the ITS phylogenetic tree ( Figure 10). Thus, we have three geographical groups of specimens of L. dorycnium s.l. (Western, Eastern and Turkish) and two geographical groups of specimens of L. hirsutus (Western and Eastern). We found fairly high levels of genetic diversity in these geographical groups (Table 3).
Within the Lotus dorycnium complex, maximal haplotype diversity was observed in the eastern group (0.986), but its nucleotide diversity is of middle value (0.0043). The western group, on the contrary, had maximal nucleotide diversity among the studied populations (0.0057), and middle Hd. The Turkish geographical group is characterized by minimal population diversity parameters, both Hd and Pi ( Table 3). The eastern group demonstrated a unimodal mismatch distribution, suggesting recent population expansion [22], while the other groups of the L. dorycnium complex showed multimodal distribution patterns indicating a stable population size over time.
Both geographical groups of L. hirsutus were characterized by high levels of haplotype diversity and average values of nucleotide diversity, which together with the multimodal nature of mismatch distribution suggests a long existence in the given territory. At the same time, the western group of the species has higher indicators of diversity than the eastern one.
The haplotype network does not allow for determining an exact ancestral haplotype of the Lotus dorycnium complex, because of many loops occurring in the network. The hypothetical haplotypes X and Y may be the divergence points between the studied clade (i.e., the L. dorycnium complex and L. hirsutus) and outgroups. The haplotype group #1 (L. hirsutus Eastern 1 clade) and group #3 (Lotus dorycnium Western clade) are the closest to the hypothetical haplotypes X and Y. Two L. hirsutus haplotype groups (Eastern group (#1) and Western group (#2)) and two L. dorycnium complex haplotype groups (Western group (#3) and Turkish group (#4)) have the most pronounced branching character without loops or with single loops. Other haplotypes, mainly representing L. dorycnium ssp. herbaceus, L. dorycnium ssp. germanicus, and L. hirsutus specimens, form a very complex and difficultto-interpret part of the network with several loops and three out of four shared haplotypes.  Table 2 and Appendix A for voucher information.  The networks constructed separately for each plastid DNA marker (not presented), despite differences in details, demonstrated the following common features: the presence of loops, many missing/hypothetical haplotypes (differences between haplotypes with two or more mutations), the presence of a haplotype shared by L. d. ssp. germanicus, L. d. ssp. herbaceus and L. hirsutus.

Taxonomic Identification of Specimens
For the taxonomic identification of the specimens, we used keys from several sources [1,9,15,19,20]; however, only the account of Rikli [9] covers the studied group worldwide, though it provides rather short keys and does not consider a lot of observations and collections made during last 120 years. Note that Rikli did not include in his account the peculiar Balearic endemic known at that time as Anthyllis fulgurans [23] and only much later identified as a member of Dorycnium [24]. The rest of the works used here for identification of samples are of a regional scale and do not include all taxa of the complex, which makes it difficult to compare taxa from different regions. The use of fruit characters was restricted because the majority of the specimens studied were in the flowering phase. Another problem was the ambiguity in the formulation of some quantitative features. For example, when the authors use the calyx tube length [1,9,19,20], they do not specify whether the tube length includes a hypanthium. In a few cases, for taxonomic identification, we considered the geographical location and phylogenetic position of a specimen.
The results of our morphometric analyses demonstrated that the Lotus dorycnium complex can be subdivided into two groups by flower length and the number of flowers per umbel, which agrees with the ideas of Rikli [9]. These characters make it possible to separate ssp. herbaceus and ssp. gracilis with small flowers and multi-flowered inflorescences from the rest of the complex. Between themselves, these two subspecies (ssp. herbaceus and ssp. gracilis) differ in the length of the calyx teeth and pubescence. The use of additional traits, especially those of pubescence, distinguishes ssp. anatolicus from other representatives of the group with large flowers and few-flowered inflorescences. However, even the use of a large number of morphological characters does not allow for precise separation of ssp. dorycnium, ssp. germanicus, and ssp. haussknechtii from each other. In addition, we found ten samples of intermediate morphology between subspecies or even between two groups of subspecies, considered here as intersubspecific hybrids. Ball [1] noted the difficulty of separating taxa within the L. dorycnium complex on a purely morphological basis and treated them as subspecies of Dorycnium pentaphyllum Scop. Demiriz [19] also accepted the subspecific rank of taxa within D. pentaphyllum Scop. Hybridization between the members of the Lotus dorycnium complex is well known [12,17,18]. An introgression zone between ssp. herbaceus and ssp. germanicus in southern Moravia and western Slovakia has been described [12]. Hybrids between ssp. dorycnium and ssp. fulgurans have been discovered in the island of Minorca among the Balearic Islands, and their hybrid nature has been confirmed by molecular methods [17,18]. The available evidence suggests the subspecies rank of the studied taxa of the L. dorycnium complex.

Phylogenetic Placement of Lotus dorycnium s.l. and Relationships among Its Subspecies
The phylogenetic analyses conducted using nrITS and plastid datasets support earlier results on the phylogenetic position of the Lotus dorycnium complex within the genus Lotus [7,8]. The ITS data confirmed sister relationships between the L. dorycnium complex and L. hirsutus.
The structure of the clade of the Lotus dorycnium complex in the nrITS phylogenetic trees reflects its geographical differentiation. Three of four large subclades (i.e., the Western, Eastern, and Turkish ones) represent three geographical groups of specimens. On the other hand, specimens of ssp. germanicus, distributed in an area that partially overlaps with that of ssp. herbaceus do not form a statistically supported clade. Strong or weak segregation of eastern and western evolutionary lineages within widely distributed Mediterranean species has been described for many groups of plants and animals [25,26]. However, clustering of ssp. fulgurans from the western Mediterranean (Balearic Islands) with the Turkish clade is inconsistent with geography. Unfortunately, our study did not include any samples of the L. dorycnium complex from the middle parts of the Mediterranean region, such as the Italian Peninsula and the large islands Corsica, Sardinia, and Sicily. Further investigation of the material from these regions may shed light on the evolutionary history of ssp. fulgurans. Interestingly, ssp. fulgurans resembles in habit another local endemic taxon of the Balearic Islands that also belongs to the tribe Loteae, Anthyllis hystrix (Willk. ex Barceló) Cardona, Contandr. et Sierra. Both plants are thorny shrublets (ssp. fulgurans is the only thorny member of Lotus). Anthyllis hystrix appears to share some aspects of its phylogeography with L. d. ssp. fulgurans. Like L. d. ssp. fulgurans, A. hystrix has all its potential relatives occurring to the east of its range. Plastid and nuclear data revealed a remarkable incongruence regarding the position of A. hystrix [27], and this result, together with the octoploid chromosome number, supports its possible hybrid origin [28]. However, the chromosome number of L. d. ssp. fulgurans, is the same as in other members of the L. dorycnium complex, 2n = 14 [29].
Within the highly supported Western subclade of the L. dorycnium complex, the specimens of ssp. pentaphyllum and ssp. gracilis are not separated, which implies the presence of a gene flow between the two subspecies. The same can be addressed to a pair of taxa from Turkey, ssp. anatolicus and ssp. haussknechtii.
The phylogenetic analyses by plastid DNA dataset demonstrated a clear separation of a clade consisting of the L. dorycnium complex and L. hirsutus in the phylogenetic trees. At the same time, they showed the impossibility of separating these two taxa based on plastid data, which confirms the results obtained earlier using a more restricted material [7,8,21]. As a whole, plastid markers are of limited value for defining taxonomic boundaries within the studied complex. Moreover, the geographic structure within the L. dorycnium complex revealed by nrITS is much less pronounced in the analysis of plastid DNA.

Phylogeography of the Lotus Dorycnium Complex
Many loops present in the plastid haplotype network suggest homoplastic mutations that are common in plastid sequences [30,31]. These loops and the general network pattern do not imply a single ancestral haplotype of the Lotus dorycnium complex. It is possible that the origin of the complex is associated with a number of hybridization events. Lotus hirsutus is obviously the closest species to the L. dorycnium complex. Their relationships are not completely resolved in analyses of plastid data, but rather are resolved in ITS analyses. This can be explained by the lower evolutionary rate of plastid sequences, which, given the recent evolution of the group [8], may not be sufficient to differentiate genetic lines. In contrast, the concert evolution of nrITS may play an important role in the evolution of low-level taxa, acting as a process analogous to lineage sorting [31]. Most of the shared haplotypes are internal, and the derived haplotypes are subspecies-specific, suggesting retention of ancestral variability and incomplete lineage sorting, e.g., [32,33] rather than recent hybridization. However, it is not so easy to distinguish between these two processes.
The haplotype network of the L. dorycnium complex is branched, with many missing/hypothetical haplotypes and a predominance of singletons. The network of another complex of species from the genus Lotus, the L. corniculatus complex, contains several widespread haplotypes, each of which has several derived haplotypes, mainly differing in one mutation [34]. Despite the differences in methodology (the trnL-F plastid DNA region in L. corniculatus and three plastid DNA regions in L. dorycnium), it is possible to outline several important differences between these networks, associated with the different histories of these complexes. It is assumed that the origin of both complexes is associated with the Mediterranean, and then the members of the complexes spread to more northern and eastern regions. Some representatives of the L. corniculatus complex (for example, L. krylovii Schischk. & Serg.) have migrated much further to the north and east and have undergone a recent expansion there, as evidenced by the presence of widespread haplotypes and a low number of derived haplotypes. In contrast, most subspecies of L. dorycnium apparently existed for a long time in the Mediterranean region, undergoing fluctuations in abundance, as evidenced by the presence of many missing haplotypes and the multimodal distribution of pairwise substitutions. L. d. ssp. herbaceus is the subspecies, the most advanced to the east. We hypothesize that it may have undergone a relatively recent expansion, as evidenced by the unimodal mismatch distribution.
Remarkably, the Western clade of L. hirsutus and the Western clade of L. dorycnium are both well-supported in the plastid phylogeny, whereas eastern accessions of both species are intermixed with each other. Nuclear data, again, show a well-supported western clade in L. dorycnium (though it does not include ssp. fulgurans). These data indicate that the evolutionary processes responsible for data incongruence between plastid and nuclear sequences of these two morphologically well-defined species were most likely localized in the eastern part of the Mediterranean region.
In the eastern Mediterranean, in southern Turkey, two other members of the L. dorycnium complex with incongruent position in phylogenetic trees constructed using plastid and nuclear data were discovered in this study. These two samples, 4980 and 9391, were identified as L. dorycnium ssp. haussknechtii, which agrees with their clusterization within the L. dorycnium Turkish clade in the ITS trees. However, according to plastid data, these specimens turned out to be close to L. graecus L., a species not included in the L. dorycnium complex. This result suggests gene exchange between the two taxa occurring in recent or older times.

Plant Material
The molecular study involved 122 specimens, including 61 specimens of Lotus dorycnium complex, 21 specimens of L. hirsutus, 37 specimens of other Lotus species representing all main sections of the genus, and 3 specimens of genera Cytisopsis, Hammatolobium, and Tripodion, closely related to Lotus. Samples for molecular studies were taken from herbarium specimens stored in herbaria ANK, GAZI, LE, MA, MHA, MW, P, and ZA. Voucher information and GenBank accession numbers are presented in Table 2 and Appendix A. Geographical distribution of specimens included in molecular analyses is presented on a map ( Figure 8) prepared using SimpleMappr [35].
The morphological study was conducted on herbarium specimens and involved 89 specimens belonging to the Lotus dorycnium complex stored in herbaria GAZI, ISTE, LE, MA, MHA, MW, P, and ZA. Voucher information of specimens included in the morphological analysis only is presented in Appendix B.

DNA Extraction, Amplification, and Sequencing
DNA was extracted from herbarium specimens (ca. 20 mg of leaf tissue) with Nucle-oSpin Plant II kit (Macherey-Nagel, Germany) according to the manufacturer's instructions or using the CTAB method [36]. The nrDNA ITS and plastid DNA psbA-trnH intergenic spacer, trnL-trnF intergenic spacer and trnL intron, and rps16 intron were selected for the analysis. The sequences of the nrITS were amplified with primers NNC-18S10, C26A [37], ITS2, and ITS3 [38]. The amplification of the psbA-trnH spacer was conducted using primers trnH2 [39] and psbAF [40]. The sequences of the trnL-trnF region of plastid DNA were amplified using standard primers 'c', 'd', 'e' and 'f' [41], and the sequences of rps16 intron using primers rpsF, rpsR2 [42], Lot-rps16-F and Lot-rps16-intR [8]. PCRs were performed in a 0.02 mL mixture containing 10-20 ng DNA, 3. PCR products were purified using the Cleanup Mini kit (Evrogen, Moscow, Russia) and then used as a template in sequencing reactions with the ABI Prism BigDye Terminator Cycle Sequencing Ready Reaction Kit v. 3.1. Sequencing was performed on the ABI PRISM 3100 genetic analyzer (Applied Biosystems, Foster City, CA, USA). 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 an electropherogram.
The sequences were aligned using MAFFT version 7.215 [43,44] and then adjusted manually in BioEdit version. 7.2.5 [45]. The matrices of psbA-trnH spacer, rps16 intron, and trnL-F plastid DNA regions were combined into a single matrix. Gap-rich and ambiguous positions were excluded from the analyses. The aligned data matrices are presented in the online Supplement (Datasets S1-S2).

Phylogenetic Analyses
Maximum likelihood analyses were performed in IQ-tree version 2.1.1 [46], internal branch support was assessed using the ultrafast bootstrap [47] with 10,000 re-samplings. The GTR + R model of nucleotide substitutions for plastid data and the SYM+ Γ model for nrITS were selected as the most appropriate by the Bayesian information criterion in the built-in ModelFinder utility [48]. In addition, ML analysis with 500 nonparametric bootstrap re-samplings was performed and a majority rule consensus tree was constructed for both data sets in RAxML version 8.2.10. The GTR + G model was used and each tree search procedure started with a random tree [49].
The Bayesian inference was performed using MrBayes v. 3.2.6 [50] considering the optimal model of nucleotide substitutions selected by AICc in PAUP version 4.0a [51] for each marker: SYM + Γ (symmetrical model with substitution rate heterogeneity) for nrITS, and GTR + Γ for plastid data. The Bayesian analysis used four independent runs of 25 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 [52]. 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 plastid DNA haplotypes were reconstructed using statistical parsimony analysis as implemented in TCS v1.2 [53]. Long indels were reduced to one character, then gaps were treated as fifth state. L. rectus, L. strictus, L. hirsutus, L. graecus, and L. corniculatus were used as outgroups. At the first stage, the haplotype networks were constructed separately for each plastid DNA marker (not presented), then for the concatenated set of three markers (trnL-F, rps16 intron, and psbA-trnH). Parameters of genetic variability were calculated using DnaSP 6 software [54].
To test the hypothesis about subdivision of the L. dorycnium complex into two main groups according to flower length and number of flowers in a head [9], we constructed a 2D scatterplot with these two characters. Then, to test the hypothesis about the separation of taxa within the L. dorycnium complex by a set of quantitative morphological characters, we conducted a discriminant analysis (DA) using Statistica v.7.0 for Windows [55]. Finally, to test the hypothesis about the separation of taxa within the studied complex by a set of quantitative and qualitative characters, we performed a principal coordinate analysis (PCoA) using PaSt 4.08 software [56]. For the PCoA method we used the Gower distance metric, which is suitable for a combination of quantitative and qualitative characters.

Conclusions
Our phylogenetic and phylogeographic study of the Lotus dorycnium L. (=Dorycnium pentaphyllum Scop.) complex revealed a tendency towards a geographical differentiation into Western, Eastern (more precisely, north-eastern) and Turkish groups supported by nrITS data. The analysis of the same set of specimens using plastid markers demonstrated a low resolution both between the L. dorycnium complex and L. hirsutus and among the taxa of the L. dorycnium complex. Interestingly, our plastid phylogeny also revealed some geographical differentiation, but in a different way. Namely, our plastid tree has a wellsupported clade that combines accessions of L. dorycnium s.l. and L. hirsutus from the eastern parts of their ranges (including Turkey), whereas western accessions of the two species form separate clades.
The discordance between our plastid and nuclear data can be interpreted as evidence of an incomplete lineage sorting and/or hybridization. The only potential way of further improving plastid phylogenetic data is though the generation of numerous complete plastid genomes to see whether current features of plastid phylogeny can be partly explained by the still inadequate informativeness of the markers used so far. However, experience from other plant groups suggests that phylogenies inferred from complete plastid genomes may still be incongruent with morphology and taxonomy. A recent study of Ophrys (Orchidaceae) in Europe and the Mediterranean revealed that plastomes represent geographic location more strongly than taxonomic assignment and correlate poorly with morphology, suggesting widespread plastid capture and possibly post-glacial expansion from multiple southern refugia [57].
We propose to treat the taxa of the complex as subspecies of Lotus dorycnium L. Presumed recent and more ancient hybridization apparently played an important role in the formation of the up-to-date pattern of genetic variability of the L. dorycnium complex and made it difficult to establish the ancestor (or ancestors) of this group.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/plants11030410/s1, Figure S1: Maximum Likelihood phylogenetic tree of the genus Lotus based on nrITS dataset (constructed using IQ-tree software). Figure S2: Maximum likelihood phylogenetic tree of the genus Lotus based on the combined plastid dataset (constructed using IQ-tree software). Figure S3: Maximum likelihood phylogenetic tree of the genus Lotus based on nrITS dataset (constructed using RAxML software). Figure S4: Maximum likelihood phylogenetic tree of the genus Lotus based on the combined plastid dataset (constructed using RAxML software). Dataset S1: nrITS dataset. Dataset S2: combined plastid dataset (trnL-F, rps16 and psbA-trnH).