Taxonomic and molecular characterization of 15 wild-growing tulip species of Greece using the internal transcribed spacer (ITS) nuclear marker in combination with the psbA-trnH and trnL/trnF plastid markers

Abstract

Due to the absence or scarcity of DNA sequence data in GenBank regarding Mediterranean, Balkan or Greek native tulip species, we performed molecular characterization of the 15 wild-growing Greek tulip species through analysis of the internal transcribed spacer nrDNA and the psbA/trnH and trnL/trnF cpDNA regions. To facilitate conservation and sustainable utilization efforts, the aim was to identify and complement the genetic distinctiveness of six Greek endemic and five Balkan or Aegean endemic and/or subendemic wild-growing Greek native tulips and two naturalized ones of Asiatic origin, among which seven are threatened with extinction. Molecular phylogenetic analysis of 148 DNA sequences was conducted, and 85 Tulipa taxa (species and subspecies) were aligned and compared to 40 new DNA sequences generated herein for the 15 wild-growing Greek tulips. These three molecular markers resulted in 354 single-nucleotide polymorphisms (SNPs) and 20 distinct indels for all Tulipa taxa, of which 100 SNPs and 12 indels were related to the Greek tulips. Phylogenetic dendrograms were constructed for all the above-mentioned molecular markers, and the DNA sequences were analyzed separately and jointly, leading to three distinct clusters for the 15 wild-growing Greek tulips. Group 1 included Greek members of subgenus Tulipa, that is, T. scardica, T. undulatifolia, T. rhodopea, T. agenensis and T. raddii; Group 2 included only T. clusiana (subgenus Clusianae) and Group 3 included Greek members of subgenus Eriostemones, that is, *T. bakeri, T. saxatilis, T. australis, *T. cretica, *T. goulimyi, *T. orphanidea T. bithynica, *T. hageri and *T. doerfleri (asterisks indicate single-country endemics).

Keywords:

• Barcoding
• PCR
• botanical collections
• tulip bulbs
• SNPs
• indel
• identification
• phylogenetic analysis

Introduction

All tulips (wild germplasm, botanical tulips and cultivated varieties) are elegant flowering plants that hold significant cultural and economic importance. Cultivated for centuries since Ottoman times, tulips have long been admired for their vibrant colours and distinctive shapes, rendering them worldwide popular ornamental flowers to date [1,2]. The tulip’s cultural and economic significance can be traced back to the eighteenth century and the ‘Tulip Era’ through the Dutch Golden Age and the ‘Tulip Mania’ [1,3], and still continues to date [2]. Despite being commonly cultivated across the world, the members of genus Tulipa L. (Liliaceae) naturally occur in the southern parts of Europe, in northern Africa, in the Middle East and across many of the mountain ranges of Central Asia from Transcaucasia to China [4]. This genus consists of numerous species; according to the checklist of the Plants of the World Online, there are 94 species of Tulipa with accepted names [5], while the World Flora Online (WFO) Plant List refers to 614 scientific names of Tulipa members, among which, only 92 species are taxonomically accepted to date [6]. The latter is related to the imperfect current knowledge regarding the circumscription of extant biological entities in the genus Tulipa and mirrors the complex taxonomic issues involved.

Apart from and/or in aid of conventional taxonomy, DNA barcodes have been commonly used to date to reveal cryptic species, identify a vast range of environmental samples, compare the phylogenetic diversity of closely related biological entities, resolve taxonomic issues and detect the genetic variations among different organisms or taxonomic units. Moreover, such approaches are also used as applied taxonomic-conservation tools for identification of plant species and subspecies as well as for conservation prioritization [7]. DNA barcoding is a useful tool in modern plant phylogenetic studies, enabling researchers to disentangle the evolutionary relationships among different taxa [8], using a significant number of acceptable DNA barcodes [9] with universal oligonucleotide primers [10]. The Internal Transcribed Spacer (ITS) region, the psbA/trnH and the trnL/trnF intergenic spacer are widely used DNA markers to date, allowing researchers to analyze plant relationships thoroughly and accurately, thus providing valuable guidance for targeted conservation efforts [11]. Tulip-wise, molecular phylogenetic analyses have been also performed using the ITS, psbA/trnH and trnL/trnF molecular markers, albeit to a limited extent compared to the importance of tulips in the global ornamental market. Currently, there are only 64 submissions available in GenBank for the molecular marker trnL/trnF, covering 49 Tulipa taxa (species and subspecies), and 52 accession numbers for the molecular marker psbA/trnH corresponding to 20 Tulipa taxa [12]. The ITS region in Tulipa species has been more extensively studied, with over 500 accession numbers corresponding to 84 Tulipa members (species, subspecies, varieties and hybrid cultivars) [12]. To gain genetic insights into important tulip groups, whole plastid genomes have been completed for nine Tulipa species of Asiatic origin, that is, T. altaica Pall. ex Spreng., T. alberti Regel, T. biflora Pall. (including T. buhseana Boiss.), T. fosteriana W. Irving, T. gesneriana L., T. iliensis Regel (including T. thianschanica Regel), T. sinkiangensis Z. M. Mao, T. suaveolens Roth (including T. schrenkii Regel), or of Eurasiatic origin such as T. sylvestris L. (including. T. patens C. Agardh) [13–15] (GenBank: NC063594; OL350836; OP749977; NC070244; OR458821). The latter indicates that all previous phylogenetic analysis of members of the genus Tulipa have primarily focused on Asiatic species, rather than on species occurring in other parts of the natural range of tulips, such as southern Europe and the Eastern Mediterranean. Despite this limitation, these studies have also provided alongside some genetic data regarding a tulip species found exclusively in Greece (i.e. T. cretica Boiss. & Heldr.); a subendemic species centred in Aegean area (T. saxatilis Sieber ex Spreng.); few species with (sub)Balkan distribution range such as T. orphanidea Boiss. ex Heldr. sensu lato (including T. hageri Heldr.), T. scardica Bornm. and T. undulatifolia Boiss.; or a couple of mainly Asiatic species with marginal populations in parts of the eastern Mediterranean (i.e. T. clusiana Redouté, T. agenensis Redouté) or the southwestern Mediterranean region (T. sylvestris) [1,16–20] (GenBank: MH555154-6; MH555199; MH555200-2; MH555212-4; MH555181-2; specimens originating from Greece, Turkey, Kosovo, Iran or Tunisia, respectively). Despite being valuable, these reference DNA sequences and available genetic information in GenBank clearly outline the scarcity of data for several tulip species.

Nested in the fringes of tulips’ natural distribution in Europe, Greece hosts from 9 up to 15 wild-growing tulip species, all protected by the Greek Presidential Decree 67/1981 [2]. However, different approaches in species’ circumscriptions and descriptions such as POWO [5] and Euro + Med [21] yield to date varied counts of tulip taxa in Greece (n = 9, though not consistent between the two). Two primary reasons contribute to this variability: (1) the circumscription of a ‘hyper-species’ namely T. orphanidea Heldr. sensu lato ranging from Bulgaria and Greece to Turkey [5], which includes as synonyms or cryptic species several tulips such as T. bithynica Griseb. ex Baker, T. doerfleri Gand., T. goulimyi Sealy & Turril, T. hageri Heldr. and T. orphanidea Heldr. sensu stricto (all considered local endemics of Greece except T. bithynica); as well as (2) the absence of established taxonomic distinctions between closely related taxa such as T. saxatilis and T. bakeri A. D. Hall or between T. agenensis and T. raddii. However, the Flora of Greece web (2023) adopts a different approach, suggesting the maximum number of Tulipa taxa in Greece (n = 15) by embracing a narrow species concept and clear taxonomic distinction of all these Greek endemic tulips (T. bakeri, T. saxatilis and T. hageri, T. goulimyi, T. doerfleri, T. orphanidea sensu stricto), or between the Asiatic tulips T. agenensis and T. raddii as well as by upgrading T. bithynica as a distinct Aegean subendemic species. Integrating the above-mentioned taxonomic discrepancies, most of these biological entities seem to refer to range-restricted plant species, many of which have been assessed as threatened according to the criteria of the International Union for the Conservation of Nature [22] such as the Critically Endangered T. bakeri and T. doerfleri; the Endangered T. cretica, T. hageri and T. orphanidea; as well as the Vulnerable T. doerfleri and T. goulimyi [2,23–25].

Due to lack of solid taxonomic basis and absence of integrated genetic knowledge about these biological entities, and in combination with (1) the national legal status of the Greek tulips (all species protected by the Greek Presidential Decree 67/1981), (2) the specialized habitat preferences of the Greek tulips [26], (3) their designated range-restriction as of Flora of Greece web (2023) and (4) the extinction risk associated with many of them [2,23–25], the present study aimed to generate missing molecular data regarding all wild-growing Greek tulips. Obtaining evidence on the genetic identity of these biological entities is expected to fill in or bridge previous knowledge gaps; to provide insight into respective taxonomic issues; to facilitate targeted and prioritized conservation efforts and to enable sustainable utilization strategies in the future, including artificial selection for breeding purposes.

Materials and methods

Plant material

The samples of 14 out of the 15 Greek tulip species were originally collected from wild-growing populations in different geographical areas of Greece [2,26], specimens of which have been taxonomically identified using standard floras [27]. Living materials of Greek tulips are currently conserved ex-situ at the premises of the Institute of Plant Breeding and Genetic Resources (IPBGR), Hellenic Agricultural Organization Demeter (ELGO-Dimitra) in Thermi (metropolitan Thessaloniki), Greece, each with an exclusive IPEN (International Plant Exchange Network) identifier (encoded accession number) as follows: T. agenensis GR-1-BBGK-22,2; T. australis GR-1-BBGK-20,148; T. bakeri GR-1-BBGK-19,346; T. bithynica GR-1-BBGK-20,515; T. clusiana GR-1-BBGK-22,46; T. cretica GR-1-BBGK-20,484; T. doerfleri GR-1-BBGK-19,354; T. goulimyi GR-1-BBGK-20,149; T. hageri GR-1-BBGK-21,46; T. orphanidea GR-1-BBGK-20,502; T. raddii GR-1-BBGK-20,512; T. rhodopea GR-1-BBGK-20,138; T. saxatilis GR-1-BBGK-20,315; T. scardica GR-1-BBGK-22,130 and T. undulatifolia GR-1-BBGK-20,144.

The bulbs of Greek tulips were maintained in 3.6 L plastic pots using the ‘Premium Soil’ (Anagnostou S. A., Athens, Greece) as suitable substrate for all species, and all bulbs were replanted in a new quantity of this substrate yearly. Plant material for DNA analysis was collected in spring from young leaves of ex-situ cultivated plants and was maintained in a freezer at −20 °C for all specimens examined for DNA barcoding studies. The specimen of T. scardica GR-1-BBGK-22,130 was sourced from well-documented originally wild material (Location: ca. 10 km from Pepeliste in Negotino-Lakavica Konecka planina, North Macedonia; Code: S & Z88-060A, B, C; Collectors: Strindberg & Zetterlund; Date: 21 April 1988) which was ex-situ cultivated at the premises of the Gothenburg Botanic Garden, Sweden.

Molecular analysis

Genomic DNA was isolated from frozen leaf tissue (ca. 250 mg) following a modified version of the Cetyl-Trimethylammonium bromide DNA extraction protocol [28]. The visual examination of the genomic DNA was conducted to evaluate both its qualitative and quantitative aspects. This was achieved through the electrophoresis of the DNA on a 0.7% agarose gel for ca. 1 h at 123–130 V.

The PCR procedure, the reagents, the sequencing procedure of PCR products and the respective analyses of DNA sequences, followed the methodology and use of technical tools as described by Tsoktouridis et al. [29]. In brief, each PCR reaction had a total volume of 50 μL and comprised of 30 μL of MasterMix (1.8 mM MgCl₂, 1 pmol/L of each oligonucleotide primer, 5 μL of 10× reaction buffer, 0.3 μL of dNTPs mix [dNTP set of 100 mmol/L, Cat: 110-012-4; GeneON GmpH, Ludwigshafen, Germany] and 25 ng of DNA template). The amplification was conducted using 3 units/reaction Maximo Taq DNA Polymerase (CAT: S101; GeneON GmpH, Ludwigshafen, Germany) and the procedure followed the manufacturer’s instructions. The oligonucleotide primers (5′ to 3′) synthesized for the amplification of the molecular markers (Eurofins Genomics, Ebersberg, Germany) were: CCTTATCATTTAGAGGAAGGAG and TCCTCCGCTTATTGATATGC for ITS (720 bp; [30]); CGCGCATGGTGGATTCACAATCC and GTTATGCATGAACGTAATGCTC for trnH/psbA (411 bp; [31]) and ATTTGAACTGGTGACACGAG and CGAAATCGGTAGACGCTACG for trnL/trnF (879 bp; [32]). The amplification protocol consisted of 4 min of denaturation at 96 °C followed by 35 cycles of 15 s at 95 °C (denaturation), 10 s (annealing), 30 s at 72 °C (extension), with a final 10-min extension at 72 °C. The annealing temperature for each pair of primers ranged from 58 to 62 °C, with the specific temperature determined based on the set of primers in each case.

The PCR products were purified using the Nucleospin Gel and PCR clean-up purification kit (CAT: 740609.50; Macherey-Nagel^TM) and the quantification of the purified PCR products was achieved using a 1.2% agarose gel electrophoresis by comparing the purified PCR products to known amounts of DNA (100 bp DNA ladder, Cat: 306-005; GeneON GmpH, Ludwigshafen, Germany). The sequencing for each amplified PCR product was performed directly by VBC-BIOTECH Service GmbH (Vienna, Austria). The PCR product was subjected to bidirectional sequencing using the amplification primers, ensuring comprehensive coverage of both DNA strands.

Phylogenetic analysis

A total of 40 new DNA sequences for the 15 wild-growing Greek tulip species were deposited in the GenBank (https://www.ncbi.nlm.nih.gov/GenBank/, accessed on 8 February 2024), each with assigned unique accession number corresponding to one of the respective molecular markers mentioned above. DNA sequences regarding the above-mentioned molecular markers for other Tulipa taxa were also retrieved from the GenBank. Pairwise genetic distance and Neighbor-Joining trees were constructed using the Maximum Composite Likelihood substitution model within the MEGA 11.0 software environment (Molecular Evolutionary Genetics Analysis; [33]).

Results

A dataset comprising 40 DNA sequences from 15 wild-growing Tulipa species of Greece (Figure 1) was genetically analyzed using the psbA-trnH, trnL-trnF plastid regions and/or in combination with the ITS (18S rDNA gene [partial] ITS1, complete 5.8S rDNA gene, ITS2 and a small part of 26S rDNA gene [partial]) of nrDNA region. This dataset was compared to all the Tulipa DNA sequences available in GenBank. Phylogenetic analysis of the psbA-trnH, trnL-trnF and ITS molecular markers was conducted separately for each molecular marker referring to 47, 58 and 83 Tulipa taxa (species and subspecies), respectively, as well as jointly for all examined molecular markers (psbA-trnH + trnL-trnF + ITS) including 32 Tulipa taxa (Table 1). Single-Nucleotide Polymorphisms (SNPs) and insertions/deletions (indels) for all taxa were initially defined (Supplemental Tables S1–S4) and then respective phylogenetic trees were constructed in each case (Figure 2, Supplemental Figures S1–S3), focusing, primarily, on the Greek tulip species. Indels and targeted SNPs variations were analyzed separately for the psbA-trnH and trnL-trnF molecular markers (Table 1), further extending the analysis to all NCBI-retrieved specimens related to the 15 wild-growing Tulipa species of Greece (Table 2).

Figure 1. Wild-growing Tulipa species documented in Greece.

Photos of T. agenensis, T. bithynica, T. goulimyi, T. hageri, T. orphanidea and T. raddii were taken by E. Dariotis; photos of T. bakeri, T. cretica, T. saxatilis were taken by M. Avramakis and the photo of T. doerfleri was taken by F. Samaritakis (all reproduced with permission). For T. scardica, see isotypus (https://iiif.bgbm.org/?manifest=https://herbarium.bgbm.org/object/B100365937/manifest.json; accessed 22 September 2023).

Figure 2. Phylogenetic tree of 34 Tulipa taxa with DNA sequences generated from wild-growing Greek tulip species (in green letters) and 22 retrieved from GenBank (in black letters). Sequence analysis was based on the ITS + psbA-trnH + trnL-trnF molecular markers and phylogenetic differences were estimated by using the Neighbor-Joining method. The evolutionary distances were computed using the Maximum Composite Likelihood method and are in the units of the number of base substitutions per site. Coding of samples and nomenclature of taxa as submitted in the GenBank database.

Table 1. Overview of nucleotide variations based on DNA sequence analysis of the psbA-trnH and trnL-trnF cpDNA molecular markers after comparison of 25 DNA sequences generated for the 15 wild-growing greek tulip species with 80 DNA sequences retrieved from GenBank (nomenclature of taxa as submitted in GenBank).

Molecular marker	Code of detected variation	Type of nucleotide variation	DNA sequences (in number of taxa detected)	Location on the T. gesneriana chloroplast genome (NC_063831.1)	Specimens of Tulipa spp. containing the nucleotide variations
psbA-trnH	INDEL-1	Insertion	CTAA (1)	151955–6	T. altaica (NC044780.1)
psbA-trnH	INDEL-2	Insertion	TAATAGTAACTTGGT (4)	151955 –6	T. australis (OR142668), T. bithynica (OR167001), T. doerfleri (OR167003), T. hageri (OR167005), T. altaica (NC044780.1), T. sylvestris subsp. australis (MZ147064.1) and T. sylvestris subsp. sylvestris (MZ147063.1)
psbA-trnH	INDEL-3	Insertion	AAACTT (1)	3 –4	T. clusiana (OR167013) and T. clusiana (EU939290.1)
psbA-trnH	INDEL-4	Insertion	TGG (17)	19–21	T. agenensis (OR167011), T. raddii (OR167007), T. rhodopea (OR167008), T. scardica (OR167012), T. undulatifolia (OR167010), T. affinis (MW847269.1), T. albanica (MZ147043.1), T. alberti (ON983980.1), T. altaica (NC044780.1), T. buhseana (NC052014.1), T. gesneriana (NC063831.1), T. iliensis (NC052697.1), T. kaufmanniana (ON423207.1), T. kosovarica (MZ147050.1), T. lehmanniana (EU939291.1), T. luanica (MZ147052.1), T. scardica (MZ147057.1), T. serbica (MZ147060.1) and T. vvedenskyi (MW847271.1)
psbA-trnH	INDEL-5	Insertion	TTACAAAATGTTAATG (1)	22–23	T. buhseana (NC052014.1)
psbA-trnH	INDEL-6	Insertion	ACTAG (1)	36–37	T. buhseana (NC052014.1)
psbA-trnH	INDEL-7	Insertion	ACTAAG (22)	117–122	T. raddii (OR167007), T. rhodopea (OR167008), T. undulatifolia (OR167010), T. agenensis (OR167011), T. scardica (OR167012), T. clusiana (OR167013), T. clusiana (EU939290.1), T. uniflora (EU939292.1), T. sinkiangensis NC063595.1, T. fosteriana (NC070244.1), T. gesneriana (NC063831.1), T. altaica (NC044780.1), T. schrenkii (NC063594.1), T. thianschanica (MW077738.1), T. iliensis (NC052697.1), T. alberti (ON983980.1), T. greigii (ON423208.1), T. affinis (MW847269.1), T. korolkowii (MW847268.1), T. lehmanniana (EU939291.1), T. vvedenskyi (MW847271.1), T. albanica (MZ147043.1), T. kaufmanniana (ON423207.1), T. kosovarica (MZ147050.1), T. luanica (MZ147052.1), T. scardica (MZ147057.1) and T. serbica (MZ147060.1), T. nalatensis OL339419
psbA-trnH	INDEL-8	Insertion	CGAAT (3)	125–129	T. clusiana (OR167013), T. clusiana (EU939290.1), T. uniflora (EU939292.1) and T. sinkiangensis (NC063595.1)
psbA-trnH	INDEL-9	Deletion	TATGTTATGGTGTTA TGTTACAAATGAATC AAACAAAAAAA (1)	12–53	T. greigii (ON423208.1)
psbA-trnH	SNPs	Variation	CAAA (17)	27–30	T. australis (OR142668), T. bakeri (OR142669), T. bithynica (OR167001), T. clusiana (OR167013), T. cretica (OR167002), T. doerfleri (OR167003), T. goulimyi (OR167004), T. hageri (OR167005), T. orphanidea (OR167006), T. saxatilis (OR167009), T. bifloriformis (ON423210.1), T. clusiana (EU939290.1), T. dubia (ON983982.1), T. korolkowii (MW847268.1), T. patens (NC061194.1), T. schrenkii (NC063594.1), T. sinkiagensis (NC063595.1), T. sylvestris (MT261172.1), T. sylvestris subsp. australis (MZ147064.1), T. sylvestris subsp. sylvestris (MZ147063.1), T. thianschanica (MW077738.1), T. turkestanica (ON423209.1) and T. uniflora (EU939292.1), T. nalatensis (OL339419), T. cretica (AM049257.1), T. sprengeri (AM085141.1)
psbA-trnH	SNPs	Variation	TGTT (17)	27–30	T. agenensis (OR167011), T. raddii (OR167007), T. rhodopea (OR167008), T. scardica (OR167012), T. undulatifolia (OR167010), T. buhseana (NC052014.1), T. affinis (MW847269.1), T. albanica (MZ147043.1), T. alberti (ON983980.1), T. altaica (NC044780.1), T. gesneriana (NC063831.1), T. iliensis (NC052697.1), T. kaufmanniana (ON423207.1), T. kosovarica (MZ147050.1), T. lehmanniana (EU939291.1), T. luanica (MZ147052.1), T. scardica (MZ147057.1), T. fosteriana (NC070244.1), T. serbica (MZ147060.1) and T. vvedenskyi (MW847271.1)
trnL-trnF	INDEL-1	Deletion	Poly-T (56)	45716–45721	Supplementary material SL - Figure S2
trnL-trnF	INDEL-2	Insertion	CAAAAAAGGA (1)	45726–45735	T. kolpakowskiana (HF952990.1)
trnL-trnF	INDEL-3	Insertion	CTATA (2)	45919–45920	T. vvedenskyi (HF953005.1) and Tulipa × tschimganica (HF953001.1)
trnL-trnF	INDEL-4	Insertion	CGTGT (4)	46148–46149	T. schrenkii (NC063594.1), T. thianschanica (MW077738.1), T. nalatensis (OL339419.1), T. intermedia (MW828756.1) and T. korolkowii (HF952991.1)
trnL-trnF	INDEL-5	Insertion	TTTATT (23)	46257–46263	T. agenensis (JQ280384.1), T. rhodopea (OR359613), T. undulatifolia (OR359615), T. raddii (OR359612), T. gesneriana (NC063831.1), T. iliensis (NC052697.1), T. altaica (NC044780.1), T. fosteriana (NC070244.1), T. agenensis (JQ280384.1), T. albanica (MN446897.1), T. aleppensis (JQ280378.1), T. armena (JQ280379.1), T. ferganica (MW828757.1), T. julia (HF952989.1), T. kolpakowskiana (HF952990.1), T. kosovarica (MN446900.1), T. lehmanniana (EU912324.1), T. luanica (MN446905.1), T. scardica (MN446910.1), T. serbica (MN446911.1), T. suaveolens (HF952998.1), T. systola (JQ280381.1), T. ulophylla (HF953003.1), T. undulatifolia (JQ280386.1), T. -----vvedenskyi (HF953005.1) and Tulipa × tschimganica (HF953001.1)
trnL-trnF	INDEL-6	Insertion	GACTACTTTTT (1)	46290–46291	T. sinkiangensis (NC063595.1)
trnL-trnF	INDEL-7	Deletion	TTTTTTTTATCTTTT (4)	46290–46305	T. iliensis (NC052697.1), T. altaica (NC 044780.1), T. ferganica (MW828757.1) and T. kolpakowskiana (HF952990.1)
trnL-trnF	INDEL-8	Insertion	GATCTT (13)	46290–46291	T. raddii (OR359612), T. saxatilis (OR359614), T. australis (OR359607), T. orphanidea (OR359611), T. bakeri (OR359608), T. goulimyi (OR359610), T. cretica (OR359609), T. sylvestris (MT261172.1), T. buhseana (NC052014.1), T. patens (NC061194.1), T. aucheriana (HF952988.1), T. biflora (JQ280377.1), T. cretica (HF952987.1), T. humilis (JQ280385.1), T. orphanidea (JQ280382.1), T. regelii (HF952995.1), T. saxatilis (JQ280387.1), T. saxatilis (HF952996.1), T. sylvestris (JQ280380.1), T. sylvestris subsp. australis (MN446915.1), T. sylvestris subsp. sylvestris (HF952999.1), T. sylvestris subsp. sylvestris (MN446916.1), T. theophrasti (HF952994.1) and T. turkestanica (HF953002.1)
trnL-trnF	INDEL-9	Insertion	TATCTT (9)	46290–46291	T. sinkiangensis (NC063595.1), T. schrenkii (NC063594.1), T. nalatensis (OL339419.1), T. thianschanica (MW077738.1), T. aitchisonii (HF952986.1), T. clusiana (EU912323.1), T. intermedia (MW828756.1), T. montana (HF952993.1), T. linifolia (HF952992.1) and T. korolkowii (HF952991.1)
trnL-trnF	INDEL-10	Insertion	----TT (18)	46290–46291	T. rhodopea (OR359613), T. undulatifolia (OR359615), T. agenensis (OR359616), T. gesneriana (NC063831.1), T. fosteriana (NC070244.1), T. agenensis (JQ280384.1), T. albanica (MN446897.1), T. aleppensis (JQ280378.1), T. armena (JQ280379.1), T. julia (HF952989.1), T. kosovarica (MN446900.1), T. lehmanniana (EU912324.1), T. luanica (MN446905.1), T. scardica (MN446910.1), T. serbica (MN446911.1), T. systola (JQ280381.1), T. ulophylla (HF953003.1), T. undulatifolia (JQ280386.1), T. uniflora (HF953004.1) and T. suaveolens (HF952998.1)
trnL-trnF	INDEL-11	Gap	------ (6)	46290–46291	T. kolpakowskiana (HF952990.1), T. ferganica (MW828757.1), T. iliensis (NC052697.1), T. altaica (NC044780.1), T. vvedenskyi (HF953005.1) and Tulipa × tschimganica (HF953001.1)
trnL-trnF	INDEL-12	Deletion	Poly-T (56)	46301–46313	Supplementary material SL – Figure S2
trnL-trnF	INDEL-13	Insertion	GATATA (5)	46315–46316	T. sinkiangensis (NC063595.1), T. aitchisonii (HF952986.1), T. clusiana (EU912323.1), T. linifolia (HF952992.1) and T. montana (HF952993.1)

Table 2. Overview of nucleotide variations in specific loci based on DNA sequence analysis of the psbA-trnH and trnL-trnF cpDNA molecular markers after comparison of the studied specimens of the 15 wild-growing greek tulip species with DNA sequences retrieved from GenBank (YES/NO indicates the presence/absence of the nucleotide variation, and dash the absence of data for comparison).

Tulipa specimens	psbA-trnH INDEL-2	psbA-trnH INDEL-3	psbA-trnH INDEL-4	psbA-trnH INDEL-7	psbA-trnH INDEL-8	psbA-trnH SNPs CAAA	psbA-trnH SNPs TGTT	trnL-trnF INDEL-1 (number of Ts deleted)	trnL-trnF INDEL-5	trnL-trnF INDEL-8	trnL-trnF INDEL-9	trnL-trnF INDEL-10	trnL-trnF INDEL-12 (number of Ts deleted)
T. agenensis (Greece)	NO	NO	YES	YES	NO	NO	YES	YES (6)	NO	NO	NO	YES	YES (13)
T. agenensis (Turkey)	–	–	–	–	–	–	–	YES (6)	YES	NO	NO	YES	YES (12)
T. australis (Greece)	YES	NO	NO	NO	NO	YES	NO	YES (6)	NO	YES	NO	NO	YES (11)
T. sylvestris (Turkey)	–	–	–	–	–	–	–	YES (4)	NO	YES	NO	NO	YES (11)
T. sylvestris (S. Korea)	NO	NO	NO	NO	NO	YES	NO	YES (5)	NO	YES	NO	NO	YES (11)
T. sylvestris subsp. australis (Kosovo)	YES	NO	NO	NO	NO	YES	NO	YES (6)	NO	YES	NO	NO	YES (11)
T. sylvestris subsp. sylvestris (Kew, U.K.)	–	–	–	–	–	–	–	YES (4)	NO	YES	NO	NO	YES (11)
T. sylvestris subsp. sylvestris (Kosovo)	YES	NO	NO	NO	NO	YES	NO	YES (6)	NO	YES	NO	NO	YES (11)
T. sylvestris (Halle, Germany)	NO	NO	NO	NO	NO	YES	NO	–	–	–	–	–	–
T. sylvestris (Italy)	YES	NO	NO	NO	NO	YES	NO	–	–	–	–	–	–
T. bakeri (Greece)	NO	NO	NO	NO	NO	YES	NO	YES (5)	NO	YES	NO	NO	YES (13)
T. bithynica (Greece)	YES	NO	NO	NO	NO	YES	NO	–	–	–	–	–	–
T. theophrasti (Kew, U.K.)	–	–	–	–	–	–	–	YES (3)	NO	YES	NO	NO	YES (11)
T. doerfleri (Greece)	YES	NO	NO	NO	NO	YES	NO	–	–	–	–	–	–
T. clusiana (Greece)	NO	YES	NO	YES	YES	YES	NO	–	–	–	–	–	–
T. clusiana (Kew, U.K.)	NO	YES	NO	YES	YES	YES	NO	YES (2)	NO	NO	YES	NO	YES (9)
T. cretica (Greece)	NO	NO	NO	NO	NO	YES	NO	YES (5)	NO	YES	NO	NO	YES (13)
T. cretica (Kew, U.K.)	–	–	–	–	–	–	–	YES (4)	NO	YES	NO	NO	YES (13)
T. cretica (Halle, Germany)	NO	NO	NO	NO	NO	YES	NO	–	–	–	NO	NO	YES (13)
T. goulimyi (Greece)	NO	NO	NO	NO	NO	YES	NO	YES (5)	NO	YES	NO	NO	YES (13)
T. hageri (Greece)	YES	NO	NO	NO	NO	YES	NO	–	–	–	–	–	–
T. orphanidea (Turkey)	–	–	–	–	–	–	–	YES (4)	NO	YES	NO	NO	YES (12)
T. orphanidea (Greece)	NO	NO	NO	NO	NO	YES	NO	YES (4)	NO	YES	NO	NO	YES (12)
T. raddii (Greece)	NO	NO	YES	YES	NO	NO	YES	YES (6)	YES	YES	NO	NO	YES (12)
T. rhodopea (Greece)	NO	NO	YES	YES	NO	NO	YES	YES (6)	YES	NO	NO	YES	YES (10)
T. saxatilis (Kew, U.K.)	–	–	–	–	–	–	–	YES (7)	NO	YES	NO	NO	YES (13)
T. saxatilis (Turkey)	–	–	–	–	–	–	–	YES (6)	NO	YES	NO	NO	YES (13)
T. saxatilis (Greece)	NO	NO	NO	NO	NO	YES	NO	YES (7)	NO	YES	NO	NO	YES (13)
T. scardica (repatriated from Sweden)	NO	NO	YES	YES	NO	NO	YES	–	–	–	–	–	–
T. scardica (Kosovo)	NO	NO	YES	YES	NO	NO	YES	YES (5)	YES	NO	NO	YES	YES (6)
T. undulatifolia (Greece)	NO	NO	YES	YES	NO	NO	YES	YES (3)	YES	NO	NO	YES	YES (12)
T. undulatifolia (Turkey)	–	–	–	–	–	–	–	YES (7)	YES	NO	NO	YES	YES (13)

Phylogenetic analysis on the basis of the psbA-trnH molecular marker

The PCR products isolated for the psbA-trnH plastid region exhibited lengths ranging from 402 to 417 bp (locus 151721 − 174 of T. gesneriana chloroplast genome) for the 15 wild-growing Tulipa species of Greece. The DNA sequences were analyzed and deposited in GenBank, obtaining the accession numbers OR142668, OR142669 and OR167001 − OR167013. These 15 DNA sequences were then compared to 32 psbA-trnH DNA sequences of Tulipa spp. retrieved from GenBank (Supplemental Table S1) and a phylogenetic dendrogram was constructed (Supplemental Figure S1) based on 47 sequences of different Tulipa taxa. A total of 58 SNPs was detected for the psbA-trnH region, thereof 13 for the Greek tulips; in addition, 9 indels were discerned, including 5 found in the Greek tulips. These nine indels comprised eight insertions (3–16 bp) and one deletion of 41 bp for T. greigii ON423208.1. In addition, a 4-bp locus separated 26 Tulipa taxa based on the CAAA sequence, another 20 Tulipa taxa with TGTT, and T. greigii ON423208.1 containing a gap at this locus (at 27–30 position of T. gesneriana NC_063831.1 complete chloroplast genome) (Table 1).

The psbA-trnH marker revealed nucleotide differences for six of the Greek tulip species, as referred to the Flora of Greece web (2023), that is, T. australis Link OR142668, T. bakeri OR142669, T. raddii OR167007, T. undulatifolia OR167010, T. agenensis OR167011 and T. clusiana OR167013. No differences were detected in a group consisting of T. rhodopea OR167008 and T. scardica OR167012, or in a second group including T. cretica OR167002, T. goulimyi OR167004, T. orphanidea sensu sticto OR167006 and T. saxatilis OR167009, or in a third group with T. bithynica OR167001, T. doerfleri OR167003 and T. hageri OR167005. The DNA sequences found in GenBank related to tulip species in Greece included T. clusiana (EU939290.1, AM085140.1), T. scardica (MZ147057-9) and T. cretica (AM049257). All these DNA sequences were identical to the results of the present study except for the sequence of T. clusiana AM085140.1, which showed one nucleotide difference compared to the herein examined wild-growing T. clusiana OR167013 from Chios Island, Greece.

Phylogenetic analysis on the basis of the trnL-trnF plastid molecular marker

The PCR products isolated for the cpDNA trnL-trnF molecular marker ranged from 823 to 835 bp in length (locus 45622 − 46455 of T. gesneriana chloroplast genome) for wild-growing Greek tulips. The sequences generated for the herein specimens of T. bithynica, T. clusiana, T. doerfleri, T. hageri and T. scardica were of poor quality; therefore, these taxa were excluded from the study of the specific molecular marker and only 10 out of the 15 wild-growing Greek tulips were studied for the trnL-trnF plastid molecular marker. The DNA sequences were analyzed and deposited in GenBank obtaining the accession numbers OR359607 − OR359616. These 10 DNA sequences were then compared to 48 trnL-trnF DNA sequences of Tulipa spp. retrieved from GenBank (Supplemental Table S2) and a phylogenetic dendrogram was constructed (Supplemental Figure S2) based on 58 sequences from various Tulipa taxa. A total of 80 SNPs was detected for the trnL-trnF region, including 11 for the Greek tulips. Additionally, 13 indels were discerned, 7 of which were found in the Greek tulips (Table 1). These indels included two poly-Ts regions, the first one exhibiting 2–7 (Ts) deletions and the second having 6–13 (Ts) deletions, both detected in the Greek tulips.

The trnL-trnF marker resulted to distinct nucleotide differences for 10 out of the 15 Greek tulip species as referred to the Flora of Greece web (2023) (T. australis OR359607, T. cretica OR359609, T. orphanidea OR359611, T. rhodopea OR359613, T. saxatilis OR359614, T. raddii OR359612, T. undulatifolia OR359615, T. agenensis OR359616, T. clusiana EU912323 and T. scardica MN446910), while one group consisting of T. bakeri OR359608 and T. goulimyi OR359610 showed no differences.

The DNA sequences obtained from GenBank related to tulip species present in Greece belonged to T. agenensis (JQ280384.1), T. cretica (HF952987.1 and AJ810118), T. orphanidea (JQ280382.1, JQ327138 and HF952996), T. saxatilis (HF952996.1 and HF952996), T. clusiana (EU912323.1) and T. scardica (MN446910.1). However, the results of the current study revealed many differences between the wild-growing Greek T. agenensis OR359616 and the Turkish T. agenensis JQ280384 (6 SNPs, 7 gaps and a 7 bp insertion). The comparison of the wild-growing Greek endemic T. cretica OR359609 with other specimens deposited in GenBank was technically problematic since the T. cretica HF952987 accession studied in the UK had poor sequence and origin data, and the T. cretica AJ810118 accession studied in Germany was identical to the specimen studied herein but its sequence was only 260 bp in length. The Greek wild-growing T. saxatilis OR359614 differed from T. saxatilis JQ260387 (2 SNPs and 2 gaps) and T. saxatilis HF952996 (3 SNPs and 4 gaps). The herein examined Greek wild-growing T. orphanidea OR359611 – either referring to a sensu lato circumscribed species or reflecting different aggregated cryptic species or considered as a species complex (‘hyper-species’) [5] – was also different from the Turkish T. orphanidea JQ280382 (2 gaps) and to T. orphanidea HF952996 studied in the UK (1 SNPs and 3 gaps). The herein examined wild-growing Greek T. undulatifolia OR167010 was also found different compared to the Turkish specimen of T. undulatifolia JQ280386 (1 SNPs and 3 gaps).

Phylogenetic analysis on the basis of the ITS molecular marker

The PCR products isolated for the ITS (18S rDNA gene [partial] ITS1, complete 5.8S rDNA gene, ITS2 and a small part of 26S rDNA gene [partial]) ranged from 711 to 714 bp in length and were sequenced for the 15 wild-growing Tulipa species of Greece. The DNA sequences were analyzed and deposited in GenBank obtaining the accession numbers OR192962, OR192963 and OR206552 to OR206564. These 15 DNA sequences were then compared to 68 ITS Tulipa DNA sequences retrieved from GenBank (Supplemental Table S3) and a phylogenetic dendrogram was constructed (Supplemental Figure S3) based on 83 sequences from different Tulipa taxa. A total of 216 SNPs was detected for the ITS region, including 76 for the Greek tulips.

The ITS marker revealed distinct nucleotide differences for 12 of the Greek tulip species as referred to the Flora of Greece web (2023), that is, T. clusiana OR206564, T. australis OR206562, T. undulatifolia OR206561, T. rhodopea OR206559, T. raddii OR206558, T. agenensis OR206562, T. goulimyi OR206555, T. scardica OR206563, T. cretica OR206553, T. saxatilis OR206560, T. bakeri OR206563 and T. orphanidea OR206557. However, a group consisting of T. bithynica OR206552, T. doerfleri OR206557 and T. hageri OR206556 presented no nucleotide differences, and their DNA sequences were identical. Species-wise, we detected several differences (shown in parenthesis) with the DNA sequences retrieved from the GenBank regarding the herein studied wild Greek tulip species. These referred to: the Turkish specimen of T. agenensis JQ776495 that showed eight SNPs and one gap compared to the Greek T. agenensis OR206562 of the current study; the specimens of T. cretica AM180461 (5 SNPs) examined in Germany as well as T. cretica MH555156 (21 SNPs and 2 gaps), MH555154 (23 SNPs and 2 gaps), MH555155 (23 SNPs and 2 gaps) and HF952961 (2 SNPs) examined in the UK, compared to the wild-growing Greek endemic T. cretica OR206553; the specimens of T. hageri MH555200 (15 SNPs and 3 gaps), MH555199 (20 SNPs and 3 gaps), MH555202 (21 SNPs and 3 gaps) examined in the UK compared to T. hageri OR206556 of the current study; the Turkish specimens of T. orphanidea JQ776499 (15 SNPs and 1 gaps) and JQ776493 (22 SNPs and 2 gaps) as well as T. orphanidea MH555214 (16 SNPs and 1 gaps), MH555201 (22 SNPs and 3 gaps), MH555213 (22 SNPs) and MH555212 (23 SNPs) examined in the UK compared to the wild-growing T. orphanidea OR206557; the Turkish specimen of T. saxatilis JQ776498 (16 SNPs and 2 gaps) as well as T. saxatilis HF952971 (2 SNPs), MH555182 (20 SNPs and 3 gaps) and T. saxatilis MH555181 (24 SNPs and 2 gaps) examined in the UK compared to the wild-growing T. saxatilis OR206560; and finally, the Kosovar specimen of T. scardica MN336211-3.1 (3 SNPs and 1 gap) compared to T. scardica OR206563.1 of the current study. The origin and details of all these specimens are presented in Tables 1 and 2.

Phylogenetic analysis using the combination of psbA-trnH + trnL-trnF + ITS molecular markers

The DNA sequences were concatenated for the psbA-trnH and trnL-trnF plastid regions and the ITS region (18S rDNA gene [partial] ITS1, complete 5.8S rDNA gene, ITS2 and a small part of 26S rDNA gene [partial]), resulting in sequences ranging from 1936 to 1966 bp in length for the 12 wild-growing Greek tulips as referred to the Flora of Greece web (2023). These 12 DNA sequences (psbA-trnH + trnL-trnF + ITS) were then aligned and compared to those of another 22 Tulipa taxa retrieved from the GenBank (Supplemental Table S4). A phylogenetic dendrogram was constructed (Figure 2) based on these 34 DNA sequences from different Tulipa taxa. The analysis of these sequences distinguished the wild-growing Greek tulip species into three groups: the first one (Group 1) included Greek members of subgenus Tulipa namely T. scardica, T. undulatifolia, T. rhodopea, T. agenensis and T. raddii; the second (Group 2) comprised only T. clusiana (subgenus Clusianae) and the third one (Group 3) included the Greek members of subgenus Eriostemones namely T. bakeri, T. saxatilis, T. australis, T. cretica, T. goulimyi and T. orphanidea (Figure 2). Tulipa bithynica, T. hageri and T. doerfleri were not included in this classification due to the lack of trnL-trnF DNA sequence data both in the current study and in GenBank.

Discussion

The whole-genome assembly of the genus Tulipa has long been anticipated, although it is estimated to 34 Gb, that is, about 10 times larger than the human genome with highly repetitive content [34]. Until such information becomes available, traditional plastid markers and the ITS nrDNA region could serve as appropriate DNA markers in attempts to decipher the phylogeny, the genetic characterization and circumscription of different Tulipa species and subspecies. Previous phylogenetic studies have offered insight into the complete cpDNA genomes of eight Asiatic Tulipa species such as T. altaica and T. illiensis [14,35,36], T. gesneriana [15], T. sinkiangensis [37], T. biflora [35], T. schrenkii (NC063594; OL350836), T. fosteriana (OP749977; NC070244), T. alberti (OR458821) and one of Eurasiatic origin namely T. sylvestris [13,14]. Undoubtedly, the same research line is needed to disentangle the taxonomic status of other tulip species occurring in Europe, especially in those centred in the biodiversity hotspot of Greece and the Eastern Mediterranean.

Comparison of DNA sequences from different individuals, species and subspecies has grown at a fast pace over the last decade, providing immensely valuable insights into the association between genomic diversity and phenotypic variance as well as into the evolution of population DNA sequences belonging to closely related species [38]. In this study, three genetic markers (ITS, psbA/trnH and trnL/trnF) were employed to enable a molecular phylogenetic analysis for the 15 wild-growing Greek tulip species. The data presented in this investigation highlighted the first molecular study on selected Tulipa genotypes of 9–15 distinct Tulipa species naturally thriving in diverse habitats of Greece. The ITS, psbA/trnH and trnL/trnF molecular markers were selected due to their involvement in previous Tulipa DNA studies, furnishing sequences deposited in GenBank, which constitute valuable material for comparative phylogenetic analysis combined with our data. Additionally, some of these DNA sequences in the GenBank referred to Greek tulips or tulip species extant in Greece, thus enabling comparisons to be made. Although more molecular markers and genetic data are needed to resolve specific taxonomic issues and problems, this study presents an initial overview of how the 9–15 tulip species of Greece are genetically related to each other or to other Tulipa taxa as described in the literature. In total, 148 DNA sequences representing 56 Tulipa species and subspecies (associated with 85 scientific names of tulip species and subspecies) were aligned and compared to new 40 new DNA sequences generated for the 15 wild-growing Greek tulip species as referred to the Flora of Greece web (2023); in particular, 48 Tulipa DNA sequences were retrieved from GenBank for the psbA-trnH, 32 for the trnL-trnF, 68 for the ITS molecular marker and 22 for the combination of psbA-trnH + trnL-trnF + ITS. The constructed phylogenetic dendrograms for all the above molecular markers and DNA sequence data enabled their analysis accordingly.

Naturally arising indel markers appear in loci as insertions or/and deletions in DNA sequences caused by cellular repair of DNA double-strand breaks, which may be generated by DNA-damaging agents such as UV light, ionizing irradiation or metabolic by-products [39]. Indels of different sizes contribute polymorphisms among plant genomes, exerting a greater impact than SNPs on nucleotide differences between individuals of the same or different species [40] and can be used as genetic markers in natural populations, especially in phylogenetic studies [41]. Indel markers have proven successful for genetic diversity studies in rice [42], chickpea [43], wheat [44] and Arabidopsis thaliana (L.) Heynh. [45], currently serving as a valuable tool in developing species-specific DNA signatures and ensuring the purity and authentication of each cultivar against unintended contamination, further endorsing the high value of products in the agricultural industry [46].

A total of 20 distinct indels were identified herein, out of which 12 refferred to the Greek tulips. Focusing on the 15 Greek tulip species as referred to the Flora of Greece web (2023), it is worth analyzing and grouping these species in the first place, and moreover, comparing them with other specimens of extant tulips of Greece as found in GenBank. These distinct indels and nucleotide variable loci were identified only for the psbA/trnH and trnL/trnF molecular markers and not for the ITS nrDNA region (Table 2). These indels were used to initially cluster the Greek tulip species into four major groups as follows: (1) T. agenensis, T. scardica, T. undulatifolia, T. raddii and T. rhodopea (all in subgenus Tulipa), (2) T. cretica, T. orphanidea, T. goulimyi, T. saxatilis and T. bakeri (all in subgenus Eriostemones), (3) T. australis, T. bithynica, T. hageri and T. doerfleri (all in subgenus Eriostemones) and (4) T. clusiana (subgenus Clusianae). Undoubtedly, these indels may indicate problematic taxonomic identification of the tulip samples analyzed but may also unveil technical mistakes or wrong labelling regarding the DNA sequences retrieved from GenBank (see, e.g. similar discrepancies in Kloukina et al. (2020) for Sideritis syriaca L. subsp. syriaca). Such problematic cases may stem from misinterpreted nomenclatural issues confusing scientists among closely related and resembling tulip species but may also indicate absence of taxonomic identification of collected tulip samples in the first place as well as uninformed taxonomic identification or sample misidentification performed by unskilled scientists and/or amateur botanists (see also Ref. [47]). Consequently, such issues may lead to inconsistent or controversial results as compared to the data generated herein, especially provided that the latter are based on taxonomic identification of wild-growing population specimens prior to genetic analysis. Therefore, the re-examination of wild-growing or cultivated status of studied material, informed taxonomic identification of samples and/or data verification associated with the above-mentioned tulip sequences is strongly suggested (Tables 1 and 2).

The phylogenetic tree based on the psbA-trnH molecular marker used for the 15 Tulipa Greek species and compared to another 32 Tulipa species and subspecies (3 Greek tulips) from GenBank (Supplemental Table S1) revealed a distinct separation into major clades. One clade comprised 10 Greek tulip species, while the other consisted of 5 Greek species (Supplemental Figure S1). These major clades presented a 4-bp locus variation (CAAA) for T. cretica, T. orphanidea, T. goulimyi, T. saxatilis, T. bakeri, T. australis, T. bithynica, T. hageri, T. doerfleri and T. clusiana, and the second major clade showed a 4-bp locus variation (TGTT) for T. agenensis, T. scardica, T. undulatifolia, T. raddii and T. rhodopea. Specimens of T. clusiana [18,20], T. scardica [16] and T. cretica [17] have been previously investigated for psbA-trnH and were herein proved to be identical based on the datasets generated in the present study. A total of 13 SNPs and 5 indels were detected for the Greek tulips and nucleotide variations effectively distinguished T. bakeri, T. raddii, T. undulatifolia, T. agenensis and T. clusiana from other species. T. clusiana’s DNA sequence of the psbA-trnH region was notably different from other tulip species, featuring two insertions (AAACTT and CGAAT) with the first one being exclusive for T. clusiana (subgenus Clusianae) and the second one shared with T. uniflora and T. sinkiangensis (both in subgenus Orithyia), thus possibly indicating a genetic relation to tulip species of Asiatic origin and range.

The trnL-trnF marker (Supplemental Figure S2) resulted in distinct nucleotide differences for T. australis, T. cretica, T. orphanidea, T. rhodopea, T. saxatilis, T. raddii, T. undulatifolia, T. clusiana, T. scardica and T. agenensis, while a separate group consisted of T. bakeri and T. goulimyi, both with identical DNA sequences. Greek tulip species and tulip specimens from GenBank showed many differences compared to the DNA sequences generated in the present study. For the Turkish specimen of T. agenensis, a 7 bp insertion (TATTCTT) was observed [19]. Regarding T. cretica, one specimen had poor data for comparison [1] and another specimen exhibited identical sequences but only 260 bp in length [18], thus of compromised value. Nucleotide differences were also detected for T. saxatilis, T. orphanidea and T. undulatifolia in specimens originating from Turkey [19] and Greece [1] but with no further information regarding the original collection locality or the taxonomic identification of the specimens, leading to an unclear interpretation. Additionally, DNA sequences for T. clusiana from a specimen deposited at Kew Gardens, UK [20] and for T. scardica from a specimen at Kosovo [16] were also included in the present study (Tables 1 and 2).

The phylogenetic analysis regarding the ITS marker (Supplemental Figure S3) resulted in distinct nucleotide differences for T. clusiana, T. australis, T. undulatifolia, T. rhodopea, T. raddii, T. agenensis, T. goulimyi, T. scardica, T. cretica, T. saxatilis, T. bakeri and T. orphanidea. However, the ITS DNA sequences of T. bithynica, T. doerfleri and T. hageri were identical. A total of 19 Greek tulip specimens for ITS were retrieved and compared to the DNA sequences generated herein; these belonged to T. agenensis [19], T. cretica [1,17] (MH555154-6), T. hageri (MH555199; MH555200-2), T. orphanidea (MH555212-4; [19]), T. saxatilis [1,19] (MH555212-4) and T. scardica [16]. All 19 specimens were found to be different compared to the DNA sequences of the present study. Therefore, a re-examination and/or verification of these specimens is recommended for the same reasons as previously mentioned for the differences observed in indels (Supplemental Table S3).

Due to the lack of trnL-trnF DNA sequence data from the present study and in GenBank as well, the Greek tulips T. bithynica, T. hageri and T. doerfleri were not included in the jointly aligned and analyzed DNA sequences from the above-mentioned three molecular markers. This combination contrasted 12 wild-growing tulip species of Greece with another 22 Tulipa taxa retrieved for this combination of molecular markers from the GenBank (Supplemental Table S4). The phylogenetic tree obtained from the combined psbA-trnH, trnL-trnF and ITS DNA sequences (Figure 2) provided the most strongly supported tree structure for the analyzed tulip specimens. The phylogenetic tree was divided into three main clades corresponding to three of the subgenera of genus Tulipa. The first clade (Group 1) comprised distinct entities of T. scardica, T. undulatifolia, T. rhodopea, T. agenensis and T. raddii (Subgenus Tulipa), the second clade (Group 2) discerned only T. clusiana (Subgenus Clusianae) and the third clade (Group 3) included the taxa T. bakeri, T. saxatilis, T. australis, T. cretica, T. goulimyi and T. orphanidea as distinct entities of Subgenus Eriostemones (Figure 2). The molecular characterization of T. bithynica, T. hageri and T. doerfleri as referred to the Flora of Greece web (2023) was complemented through the combination of ITS and psbA-trnH molecular markers (Supplemental Figures S1 and S3), indicating that these belong to Group 3.

Despite the economic significance of tulips in horticulture or the importance of biodiversity conservation in the age of climate crisis, at present, the molecular characterization is primarily focused on Asiatic species, notably lacking for many wild tulip species of Balkan and East Mediterranean ranges, and particularly for threatened Greek endemic ones. This investigation has provided specific genetic indications regarding the distinctiveness of 13 out of the 15 Greek tulip species, thus facilitating their taxonomic circumscription and enabling future conservation efforts and sustainable utilization strategies. Certainly, more molecular markers and more data are needed to clarify their taxonomic status before prioritizing species and key areas for conservation strategies in Greece and beyond, or prior to their sustainable utilization.

Conclusions

The study presented herein, represents the first molecular characterization of all the wild-growing tulip species of Greece (n = 15) on the basis of the psbA-trnH, trnL-trnF and ITS molecular markers offering new insight into their genetic distinctiveness and allowing to make steps towards future conservation and sustainable utilization efforts. Undoubtedly, additional research is required, with more and new species-specific molecular markers to accurately circumscribe and characterize the genetic relatedness, taxonomic distinctiveness, and phylogenesis of such wild and valuable tulip germplasm.

Authors’ contribution

Conceptualization: G.T., N.K.; methodology: G.T., N.K., M.S., I.S.; validation: G.T., N.K., I.S., A.T., M.S., I.T., P.M.; formal analysis: G.T., N.K., I.S., A.T., M.S.; investigation: G.T., N.K., I.S., A.T., P.M.; resources: G.T., N.K.; data curation: G.T., N.K., A.T.; writing – original draft preparation: G.T., N.K., I.S., A.T.; writing – review and editing: G.T., N.K., I.S., A.T.; visualization: G.T., N.K., I.S.; supervision: G.T., N.K., I.T.; project administration: G.T., N.K.; funding acquisition: G.T., N.K.; All authors read and approved the final submitted manuscript.

Ethical standards

The authors declare that this study complies with the current laws of E.U.

Acknowledgements

The authors are thankful to V. Papiomytoglou and M. Avramakis (Natural History Museum of Crete) for assistance in field work and the permission to reproduce photographs of Tulipa cretica and T. saxatilis, respectively. We are also thankful to the staff of the Balkan Botanic Garden of Kroussia and the Institute of Plant Breeding and Genetic Resources, Hellenic Agricultural Organization Demeter for hosting the ex-situ conservation of the tulip accessions of TULIPS.GR. We would also like to thank Stina Weststrand and Mats Havström of Gothenburg Botanical Garden (Sweden) for the Tulipa scardica samples analyzed in this investigation.

Disclosure statement

The authors declare that they have no conflict interests.

Data availability statement

The data supporting the findings of this study are available from the corresponding author upon reasonable request.

Additional information

Funding

This research has been co-financed by the European Union and Greek national funds through the Operational Program Competitiveness, Entrepreneurship and Innovation, under the call RESEARCH-CREATE-INNOVATE (project code: T2EDK-05115; acronym: TULIPS.GR), entitled ‘Value chain for Greek native tulips: development of documented propagation material and integrated conservation for sustainable exploitation’.