A phylogenetic analysis of the wild Tulipa species (Liliaceae) of Kosovo based on plastid and nuclear DNA sequence

In Kosovo, the genus Tulipa is represented by eight taxa, most of which form a species complex surrounding Tulipa scardica. To investigate the phylogenetic relationship of these Tulipa species a Bayesian analysis was undertaken using the ITS nuclear marker and trnL‐trnF, rbcL and psbA‐trnH plastid markers. The resulting phylogenetic trees show that Kosovarian Tulipa species consistently group into two main clades, the subgenera Eriostemones and Tulipa. Furthermore, our analyses provide some evidence that the subspecies of Tulipa sylvestris are genetically distinguishable, however not significantly enough to support their reclassification as species. In contrast, the markers provide some novel information to reassess the species concepts of the T. scardica complex. Our data provide support for the synonymisation of Tulipa luanica and Tulipa kosovarica under the species Tulipa serbica. Resolution and sampling limitations hinder any concrete conclusion about whether Tulipa albanica and T. scardica are true species, yet our data do provide some support that these are unique taxa and therefore should continue to be treated as such until further clarification. Overall, our work shows that genetic data will be important in determining species concepts in this genus, however, even with a molecular perspective pulling apart closely related taxa can be extremely challenging.


| INTRODUCTION
Species of the genus Tulipa L. (Liliaceae) have great economic, horticultural, and ecological value 1 while also being culturally significant in many areas of the world. 2 They are bulbous monocots characterized by a diverse range of variable vegetative and floral traits, which were traditionally used to define species concepts in this genus. Furthermore, the vegetative and floral traits often show a high degree of plasticity, sometimes, even within populations of a species. [2][3][4] Due to this and the long horticultural history of tulips, creating a stable taxonomic framework for the genus has been extremely difficult, despite the existence of a large body of literature, [3][4][5][6]

and so classifications of
Tulipa have been revised several times, 7 The total number of extant Tulipa species varies between publications, although generally ranges from 40 to 150 species. 5,8 In the World Checklist of Selected Plant Families, 9 516 names are listed for Tulipa, but only 102 taxa have been accepted, while in the Plant List 10 499 names are listed for Tulipa and 120 taxa have been accepted. According to the most complete evaluation of the genus to date, 2 only 76 species are accepted, but since this work, a number of new species have been described. [11][12][13][14] The number of Tulipa species native to the Balkan Peninsula is only a small proportion of the global diversity, varying from 15 15 to 22 9 species. In Kosovo, the genus Tulipa is represented by eight taxa (six species and two subspecies), belonging to the subgenera Eriostemones and Tulipa.
In general, researchers working on these species have used different morphological traits to define the taxonomic relationship between them. The subgenus Eriostemones, is generally represented by Tulipa sylvestris and at the lower taxonomic level by two subspecies, 16 T. sylvestris subsp. sylvestris only accepted by the World Checklist of Selected Plant Families 9 and Tulipa sylvestris subsp. australis (Link) Pamp (accepted subsp.). While the subgenus Tulipa is represented by several species, Tulipa gesneriana L. 16 is sometimes treated as a wild species, although it is not thought to grow in a truly wild state 2 and is believed to be a complex hybrid derived from T. agenensis DC, T. armena Boiss, Tulipa suaveolens Roth and so on. 9 It is here not treated as a true species in line with previous research, but is included in a range of analyses. 2 Tulipa scardica Bornm. which has a distribution that encompasses Southern Kosovo and North Macedonia. 17 In Kosovo, it occurs near the village Krivenik, close to the border of North Macedonia. It is synonymised as T. gesneriana L., 2,6,10 accepted as a species by the World Checklist of Selected Plant Families, 9 but not accepted by Flora Europea. 18 Tulipa serbica Tatic & Krivošej occurs on serpentine soil in the South of Serbia (community Knjaževac: Mt. Rogozna near Donja Kamenica) and Northern Kosovo (Beli Laz hill, near Ibar river). 19 Tulipa kosovarica Kit Tan, Shuka & Krasniqi is endemic to Kosovo, in the serpentine area of Mirusha region at the foot of Mt. Koznik, between Mrasori and Llapçevë villages, 20 as well as in the localities Guriç, Llapushnik, Qafë Prush and Devë. 16 Tulipa luanica Millaku is also endemic to Kosovo found growing on limestone substrate on Mt. Pashtrik, located in the district of Prizren, Southern Kosovo, near the border with Albania. 11 Tulipa albanica Kit Tan & Shuka was originally described as a new species from a locality in Albania 21 (Kukësi district: from Kolshi to Surroj village, on serpentine slopes), but has been recently found growing in the Kosovar village of Deva. 16 T. scardica, T. serbica, T. albanica, T. kosovarica, and T. luanica are all morphologically similar [19][20][21] and form the species complex known as the T. scardica complex (scardica complex), named after the oldest species in the group. 2 Due to the similarities between these species, they have sometimes been treated as synonyms, and are often erroneously identified and misclassified.
Studies focused on defining species concepts within the scardica complex have primarily used morphological characteristics and geographical distributions. However, in addition, karyological analyses have been undertaken for T. albanica, 21 23 Molecular phylogenetic analysis using sequences from nuclear ribosomal DNA (nrDNA) and chloroplast DNA (cpDNA) have previously been successfully used to determine relationships between species within the genus Tulipa. Thus, we decided to use Tulipa DNA sequences from the ITS region, 2,7,24,25 trnL-trnF region, 26 psbA-trnH region, and rbcL region 27 to undertake a phylogenetic analysis of Kosovarian tulip diversity. This work aimed to improve understanding of species concepts across the wild-growing Tulipa species of Kosovo, especially the scardica complex, with a view to inform tulip conservation, evolutionary understanding, and the broader taxonomic positioning of Kosovarian tulip species.

| RESULTS
The ITS sequences (ITS1, complete 5.8S rDNA gene, ITS2 and a small part of 26S rDNA gene) of Tulipa species in the dataset ranged from 616 to 655 bp. The in-group alignment included 66 ambiguous positions. Sixty-seven positions were potentially informative, 33 potentially informative indels, and 60.0% G + C content ( The trnL-trnF sequences of Tulipa species in the dataset ranged from 765 to 788 bp in length. The in-group alignment had 46 ambiguous positions. Analyzed sequences showed eight potentially informative characters, 16 potentially informative indels, and 31.2% G + C content ( Table 1). The trnL-trnF region was made up of trnL 631 to 692 bp, trnF 57-64 bp and IGS 25 bp for each sequence, respectively.
The rbcL sequence length in the dataset ranged from 488 to 597 bp.
In-group alignment includes three ambiguous positions. Analyzed sequences showed three potentially informative characters, five potentially informative indels, and 44.0% G + C content ( Table 1). The psbA-trnH sequences in the dataset ranged from 488 to 597 bp in length. The in-group alignment had 35 ambiguous positions. Analyzed sequences showed 15 potentially informative characters, 93 potentially informative indels, and 32.6% G + C content ( Table 2). The combined ITS + trnL-trnF + psbA-trnH + rbcL sequences for species ranged from 2405 to 2469 bp in length. The alignment showed 134 ambiguous positions, 2272 conserved sites, 113 potentially informative characters, and 125 potentially informative indels, and an average 42.1% G + C content (Table 1).

| Phylogenetic analysis
In total, 106 sequences were used in the phylogenetic analysis. The phylogenetic trees for all datasets (separate ITS, trnL-trnF, psbA-trnH, and rbcL trees as well as the combined ITS + trnL-trnF + psbA-trnH + rbcL datasets) were generated through a Bayesian analyses. Resolution was relatively weak for all trees produced from single markers while the best resolution was obtained from the phylogenetic tree created using the combined ITS + trnL-trnF+ psbA-trnH + rbcL dataset.

| ITS region
The phylogenetic analysis based on 31 ITS sequences is shown in

| trnL-trnF region
The phylogenetic tree obtained from 28 trnL-trnF sequences was again divided into two major clades representing the two subgenera sampled ( Figure 2). The first clade, Eriostemones (T. sylvestris including both subspecies), was strongly supported by Bayesian analyses (BPP = 1), although the structure of the tree slightly differed from the phylogenetic tree created using ITS data; the T. sylvestris subsp. australis specimen (T22) was more closely related to a T. sylvestris subsp.

| psbA-trnH region
Twenty-six psbA-trnH sequences were used to construct a phylogenetic tree, which again resulted in the clear division of the sequences into two major clades ( Figure 3). Again, the first strongly supported clade (BPP = 1) consists of members of the subgenus Eriostemones, with, similarly to the tree generated based on ITS sequences, wild T. sylvestris subsp. sylvestris separated from T. sylvestris subsp. australis. Surprisingly, the outgroup specimen A. erythronioides fell within the Eriostemones clade suggesting this may be a poor marker for taxonomic understanding. The second strongly supported clade (BPP = 1) consists of members of the subgenus Tulipa. The specimens of subgenus Tulipa appear to divided into two groups; the first consisting of all T. albanica, specimens, while the second group encompasses specimens of T. scardica, T. kosovarica, T. luanica, and T. serbica, although none of these were monophyletic with the tree lacking structure (BPP = <0.5).

| rbcL region
The phylogenetic tree obtained from 29 rbcL sequences was the least informative single marker ( Figure 4) with generally low posterior probability scores. Even so, the marker was able to distinguish between T A B L E 1 Data set and parsimony-based tree characteristics for ITS and trnL-trnF, rbcL, and psbA-trnH analyses  Here, we also note that Erythronium japonicum appears more closely related to Tulipa specimens than to Amana specimens, which contradicts the expected relationship of these genera providing some evidence that this is not a taxonomically informative marker.
2.6 | Combined ITS, trnL-trnF, psbA-trnH, and rbcL dataset The phylogenetic tree obtained from the combined ITS, trnL-trnF, psbA-trnH, and rbcL sequences provided the most strongly supported tree structure for the specimens analyzed ( Figure 5). The phylogenetic tree is divided into two main clades, the subgenus Eriostemones and the subgenus Tulipa with strong support for this separation (BPP = 1).
Within the Eriostemones subgenus, the specimens of T. sylvestris

| DISCUSSION
In this study, we use the genetic markers ITS, trnL-trnF, psbA-trnH, and rbcL to undertake a molecular phylogenetic analysis of Kosovarian tulip diversity. Our data highlight the informativeness and limitations of the ITS nuclear marker 2,7,24,25 and plastid markers trnL-trnF, 26 rbcL, [28][29][30] and psbA-trnH 30,31 in investigating evolutionary relationships between species of wild Tulipa. In general, we found that subgenera can be reliable separated by a range of single genetic markers; however, that separating more closely related species requires a combination of markers. Our most informative tree provides evidence that the scardica complex has been over split and specifically that T. luanica In general, phylogenetic trees generated using ITS sequence data had better resolution than those generated from single plastid markers, including the trnL-trnF marker which is in line with previous research. 7,32 The rbcL tree was the least informative as it had extremely weak resolution across the analyzed taxa, which supports previous reports of the marker performing poorly. [28][29][30] The psbA-trnH marker provided somewhat better resolution than rbcL, 30 The use of sections within the genus Tulipa was actively discouraged 2 until further in-depth genetic studies could be undertaken. Yet, we wanted to briefly explore how our ITS tree fits into the taxonomic framework developed by Zonneveld. 6 The phylogenetic tree based on the ITS marker we generated had monophyletic groups that represented the Eriostemones section Sylvestres and the Tulipa sections Tulipa but do note that these may not all hold as more genetic data become available.
Unsurprisingly, our most informative tree was generated using the combined dataset that included all the markers (ITS, trnL-trnF, psbA-trnH, and rbcL). This, like the single marker trees, separated the  Further cytotaxonomic studies will therefore be needed to investigate the chromosome numbers of the specimens located in Kosovo to confirm their taxonomic identity, while extensive in-depth molecular work will be needed to unentangle this widespread, complicated taxon. In Today, there are five species recognized as part of this complex.
T. scardica was the first species described from this complex, 33 and individuals of this species show significant variation in several morphological characters, such as leaf form, flower color, length of filaments, and anthers in different distribution areas. 17 Tulipa serbica, the second species named in this complex was described from Mt Rogozna 34 and was originally thought to be a population of T. scardica, before being described as a new species. 20 Both species are thought to be closely related with T. serbica only morphologically differing from T. scardica in its paler, unspotted periapt segments, pale (not blackish) staminal filaments, dull violet (not yellowish), and acute anthers. 19 T. albanica was recorded as a new species in Northeast Albania in 2010; it has recently been found growing in Kosovo as well. 16  with radiation showed that blotch margin and flower color can easily be influenced. 35 Flower color is therefore not regarded as a suitable trait from which to make taxonomic decisions. 2 Apart from flower morphological features, the characteristics of the bulb tunic have often been used to differentiate between Tulipa species and has generally been found to be a reliable character. 36  F I G U R E 5 Phylogenetic trees based on a combined ITS+trnL-trnF+ rbcL + psbA-trnH sequence set including, including posterior probabilities (BPPs) (>0.5) provided above each node especially given that differences in genome sizes within species could be correlated with differences in habitat, 37 plant phenotype, 38 or caused by technical artifacts. 39 In addition, the DNA content of T. serbica has not been measured so this cannot be linked to other species in the scardica complex. Overall, this means that our DNA sequence data are likely the best assessment of this species complex to date and should be used as a guide on how to classify these taxa into species over and above current cytogenetic data.

| Plant material
Eight taxa (six species and two subspecies) of the genus Tulipa were collected from wild populations between the months of April and May across 2017, 2018, and 2019. All Tulipa species were collected in Kosovo, except T. albanica, which was collected in Albania ( Figure 6).
One unidentified plant specimen of Tulipa sp. (sample T9,

| Phylogenetic analyses
In total, 106 sequences obtained from 14 taxa were analyzed, 87 of them were newly generated sequences generated from eight Tulipa taxa (six species and two subspecies) collected from wild populations in Kosovo and 19 sequences were obtained from GenBank ( Table 2).
ITS, trnL-trnF, psbA-trnH, and rbcL sequences of most of the taxa were amplified and then sequenced from three specimens for each species, while the T. kosovarica (locality Goriç) and T. luanica (locality Qafë Prush) were amplified and sequenced successfully from two specimens per species. Due to the amplification failure of some specimens (ITS T8 and T10; trnL-trnF T9, T16, and T18), some species were represented by only one or two sequences.
The sequences were aligned using MEGA X software. 43 For ITS analyses, in total 31 sequences were aligned to determine sequence statistics, 21 of them were newly generated, and 10 were obtained from GenBank, for trnL-trnF statistical analyses included 28 sequences (20 newly generated and eight obtained from GenBank) ( Table 1). For rbcL analyses of 29 sequences were used, of them 23 were newly generated and six of them were obtained from Genebank, while for psbA-trnH 26 sequences were used for analyses of them 23 newly generated and three obtained from gene bank (outgroup species).
Bayesian analyses were conducted through a Markov Chain Monte Carlo (MCMC) approach using BEAST v1.10.4 with the help of BEA-GLE v3.1.0 library. The input files for BEAST were prepared in the corresponding BEAUti program and maximum clade credibility trees generated and annotated in TreeAnnotator. 44 The MCMC was run for 10 000 000 generations, with resulting phylogenetic trees sampled every 1000. A burn in period of 1 000 000 was used. All trees were visualized using Figtree (V. 1.4.4) and Mega X software.

| CONCLUSIONS
Our phylogenetic analyses show that Kosovarian tulips can easily be distinguished as either in the subgenera Eriostemones or Tulipa. Yet, within these subgenera, we found limited resolution to determine clear species relationships using the markers we selected. Nonetheless, we note that there was some genetic distinguishability between the subspecies of Tulipa sylvestris (australis and sylvestris) and that these should therefore continue to be classified as different subspecies but our work does not suggest that they should be raised to species level. In contrast, our data suggest that within the Tulipa subgenus, there has been over splitting of species within the scardica complex. With our novel genetic perspective, we suggest that T. luanica and T. kosovarica can be synonymised under T. serbica, while both T. albanica and T. scardica were genetically distinct enough to continue to be treated as species. Further analyses with more extensive sampling and additional genetic markers will be necessary for a better understanding of the natural variability within the taxa of the scardica complex, but for now our study provides the most comprehensive genetic understanding of the complex diversity of tulips growing in and around Kosovo. This understanding will not only be crucial for taxonomic stability and future research, but also for identifying conservation priorities, especially given that threats to wild tulips are likely to increase in the near future. 45