Tulipa suaveolens Roth is a decorative bulbous, polycarpic species in the Liliaceae family. Currently, the species populations are declining due to the degradation of natural steppes. In the present article, we evaluated the genetic diversity of 216 specimens from 22 natural T. suaveolens localities in the European part of the species range using 10 inter simple sequence repeat (ISSR) markers. The polymerase chain reaction (PCR) yielded 250 unambiguous and reproducible polymorphic bands with a mean polymorphism information content (PIC) value of 0.27. Neither the principal component analysis nor the clustering method split the dataset. However, the Bayesian model-based STRUCTURE analysis detected two genetic clusters of T. suaveolens. The analysis of the biogeographical distributional pattern of the revealed genetic groups confirmed our hypothesis that the evolution and distribution of T. suaveolens were likely affected by the Early Khvalynian transgression of the Caspian Sea and the following Enotayevsk regression. Furthermore, this hypothesis was supported by the results of the NEW HYBRIDS analysis identifying the specimens from the populations located outside the Early Khvalynian flooding limit as pure parents, while the specimens of the populations within the flooding limit were classified as backcrosses and second generation hybrids.

The samples were collected from 22 natural T. suaveolens localities in the Astrakhan, Volgograd, Orenburg, Rostov, Samara and Saratov provinces, the Krasnodar region, the Republics of Kalmykia and Dagestan as well as from the Crimea and Western Kazakhstan (Table 1). Ten plants were selected from each locality except the population no. 22. Leaves were dried in silica gel and put in plastic zip-bags.

DNA was extracted using the NucleoSpin® Plant II kit (MACHEREY-NAGEL, Germany) according to the manufacturer’s protocol. DNA purity (absorbance ratio A260/A280) and quantity (absorbance at 260 nm) were measured by the Qubit fluorometer (Thermo Fisher Scientific, Carlsbad, CA). The samples were normalized at 10 ng/μl for the ISSR analysis following the Qubit concentrations.

The polymerase chain reaction (PCR) was carried out in the Mastercycler gradient amplifier (Eppendorf, Germany) with 10 ISSR primers from the “Sintol” company (Moscow, Russian Federation) that were pre-selected for T. suaveolens (Kashin et al. 2016; Kritskaya et al. 2019). The PCR was performed in 20 µl reaction volumes. The reaction mixture contained 4 µl ready-to-use PCR mix MaGMix (200 µM of each dNTP, 1.5 mM MgCl₂, 1.5 U SmarTaqDNA-polymerase and buffer; Dialat Ltd., Moscow, Russian Federation), 15 µl of deionized water, 3.4 pmol of each primer and 1 µl of template DNA. The PCR amplification regime was as follows: preliminary denaturation during 5 minutes at 95° C, then 35 cycles of 30 seconds each at 95° C, 30 seconds at 44° C and 2 minutes at 72° C, then final elongation during 10 minutes at 72° C. DNA bands were visualized on 2% agarose gels (10 cm long and 15 cm wide) ran at 5V/cm for 2 h. Gel images were documented using the gel documentation system (Doc-print VX2, Germany). One Kb and 100 bp ladder were used as reference.

For each specimen, two technical replicates with each ISSR primer were conducted (Gel replicability score ~ 95%).

All DNA samples were amplified and produced clear, reproducible bands which were scored as present (1) or absent (0), using the CorelDRAW X8 software. A binary matrix (1/0) containing 250 polymorphic loci was generated and used for the subsequent data analysis for diversity and population structure.

The informativeness of the ISSR markers was evaluated using polymorphism information content (PIC), resolving power (RP), mean resolving power (MRP), the marker index (MI) and the Shannon index (H’). PIC is the probability in detecting polymorphism by a primer: PIC = 1−∑(P_i)², where P_i is the frequency of the i^th allele (Sehgal et al. 2009). RP is the ability of each primer to detect the level of variation between individuals and was calculated according to Prevost and Wilkinson (1999): RP = ∑Ib, where Ib (band informativeness) takes the values of 1–[2|0.5–p|], where p is the proportion of individuals containing the band. MI for each primer was calculated as a product of polymorphism information content and the effective multiplex ratio (MI = PIC × EMR (Varshney et al. 2007)). The Shannon index (H’) was calculated by the formula H’ = −_i=1^Rpilnpi (Shannon and Weaver 1949).

Genetic diversity within each population was estimated through the percentage of polymorphic loci (P), the mean effective number of alleles (Ne), the mean Shannon Information Index (I), Nei’s gene diversity (h'), the overall Genetic differentiation (G_ST), the Gene flow (N_m) and the number of rare alleles/loci per population using the POPGene v 1.32 (Yeh 1997) software.

The analysis of the obtained matrix was performed by the Neighbour Net method in the SPLITS TREE 4.6 software (Huson 1998). The principal component analysis (PCA) was carried out using R environment (R Core Team 2020) with RStudio (Version 1.2.5033) by running the scriptfile (see Appendix 1).

The population structure analysis was conducted by the Bayes method in the STRUCTURE 2.3 software (Pritchard et al. 2000, Evanno et al. 2005, Jakobsson and Rosenberg 2007). The analysis was performed twice based on the genetic admixture model. The burn-in took 500,000 iterations; the subsequent building of the Markov chain took one million iterations for K (the hypothetical number of populations) from 1 to 10 with five-time repetition for each value of K. The obtained data were interpreted in the STRUCTURE HARVESTER software (Earl and VonHoldt 2012).

Since we supposed that T. suaveolens populations were long isolated from one another due to the geological processes of the Late Pleistocene, we performed an additional NEW HYBRIDS 3.1.1 analysis (Anderson and Thompson 2002). The method evaluates the posterior probability of each specimen belonging to one of the probable hybrid classes and requires no preliminary knowledge about the parent classes (Anderson and Thompson 2002). NEW HYBRIDS is also used to reveal the restrictions of the gene flow between populations of the same species (Misiewicz and Fine 2014, Palm et al. 2019, Reid et al. 2019, Ortíz-Gamino et al. 2020) resulting from the physical barriers separating populations during the species distribution. All data were analyzed ten times from overdispersed starting values, over one million sweeps after a burn-in of 500,000 sweeps, following the recommendations of Anderson and Thompson (2002).

The ArcGIS® 10.6 software was used to analyze the geographical distribution of the revealed genetic groups.