In total, pods of 74 individuals including 62 individuals from eight populations of M. sinskiae, five individuals from a single population of M. constricta and seven individuals from a single population of M. rigidula were collected from western and southwestern Iran in July 2017 (Table 1). Considering the limited occurrence of M. sinskiae, the restricted population size, and the criterion of individuals of each population being sampled at 20 m minimum intervals, no more effective number of individuals could be gathered. We could not find any material from Turkmenistan. All collected samples were identified and labeled based on authoritative identification keys (Heyn, 1984; Small & Jomphe, 1989; Mehregan et al., 2002; Small, 2011). Six to eight seeds of each individual were cultivated in separate pots on a research farm in southwestern Iran with similar ecological conditions to their natural habitats. Samples were taken at different growth stages up until fully ripened pods were developed. Total DNA was extracted from young leaves. Morphological examinations were performed on the fully grown plants and pods. Voucher specimens were deposited at IAUH (Islamic Azad University Herbarium).

Total genomic DNA was extracted from young leaves dried in silica gel using CTAB (cetyltrimethylammonium bromide) method of Doyle & Doyle (1987) employing Nucleospin© Plants kits (Machery-Nagel, Germany) after manufacturer’s instructions. Success of DNA extraction was initially checked on 1% agarose gel. Density and purity of extracted DNA were examined on a NanoDrop™ 2000 (Thermo Scientific).

Method of Yang et al. (2015) was used to identify and develop SSR markers. 100 ng of the genomic DNA of a single sample was used to generate an Illumina DNA library. After DNA was fragmented, repaired at the ends, “A-tailed”, and ligated to the TruSeq adapters, the library was amplified in eight cycles. The average size of the library was 670 bp, corresponding to an average integral length of 500 bp. Sequencing of the “library” was carried out in an Illumina-MiSeq system (Illumina, San Diego, CA) with 300 bp each in the “paired-end” mode. In order to remove residues of adapter sequences, the overlapping “paired-end-reads” were first trimmed at the ends. The quality score was set to at least 20. FLASH software (Magoč & Salzberg, 2011) was used to assemble the reads. The resulting sequences were bioinformatically analyzed for existing microsatellites (Faircloth, 2008). Examination of 817,510 potential loci with any tandem repeat were performed with following criteria: mismatch = 0, motive-length = 3–5 and length of repeats = 40–185.

Primer3 software (Koressaar et al., 2018) was used to design primers for 513 out of 2410 eligible loci suitable for primer designing with the following criteria: PCR-product size = 150–350 bp, primer length = 18–22 bp, and TM-value = 58–62°C. Primer pairs were synthesized for 124 loci and tested on 12 samples for suitability. Those primer pairs which yielded unique PCR fragments were tested on a batch of four individuals. In the best case those primers should result in different unique fragments in each individual. PCR products from individuals were sequenced to confirm the tandem repeat pattern. Out of the 124 primer pairs only six polymorphic loci could be detected. Primers with different labeling were synthesized (Table 2). In order to avoid the possible problems with null alleles, we tried to improve our methods by using labelled primers, finding the best annealing temperature, and sequencing the PCR products for estimating the product size (Dakin & Avise, 2004).

Polymerase chain reactions (PCR) were carried out in 25 μl reaction volumes containing 5 μl of 5× PCR-buffer, 2 μl of MgCl₂ (25mM), 2 μl of dNTPs (2.5 mM), 0.1 μl of Taq-Polymerase (1.25 U), 1 μl of forward primer (10 pmol/μl), 1 μl of reverse primer (10 pmol/μl), 1 μl of genomic DNA (~15 ng), and 12.9 μl of ddH₂O. The PCR reactions were performed on a Labcycler Gradient (SensoQuest GmbH, Göttingen, Germany) under the following conditions: initial denaturation at 95°C for 180 s, 34 cycles of denaturation at 95°C for 30 s, annealing at primer specific temperature for 30 s and extension at 72°C for 30 s, and a final extension at 72°C for 300 s. PCR products of MED-01/MED-02/MED-03 and MED-04/MED-05/MED-06 primer pairs were pooled separately. Two μl of each pool was mixed with 7.75 μl of HiDi formamide (Applied Biosystems) and 0.25 μl ROX-500 internal size standard (Applied Biosystems) and then injected to an 3730xl Applied Biosystems capillary sequencer. Raw data were visualized with GeneMarker v4.0 (Applied Biosystems, Foster City, CA, USA). Output files were aligned with the ROX-500 size standard using GeneMarker v2.4.2 (GeneMarker, SoftGenetics, State College, PA, USA). Each peak with a signal intensity of more than 1000 was scored as present. The binary matrices (1: presence, 0: absence) of each two primer pools were combined and prepared for further analyses. Some analyses need data to be entered as co-dominant. To do so, allele sizes for each locus were entered in GenAlEx software after publisher’s tutorials.

UPGMA (Unweighted Pair Group Method with Arithmetic mean) algorithm of clustering with Dice similarity index as well as PCoA (Principal Coordinate Analysis) analyses of dataset were performed with PAST3 software package (Hammer et al., 2001). After circumscribing M. sinskiae with UPGMA and PCoA, further analysis was performed on M. sinskiae populations only. In POPTREEW software genetic distances were measured and phylogenetic tree was constructed.

Genetic structure of M. sinkiae populations was estimated using a Bayesian Markov Chain Monte Carlo model (MCMC) implemented in Structure v2.3.4 (Pritchard et al., 2000). The true number of subpopulations (K) was calculated using Evanno et al. (2005) method summarized in CLUMPP_Windows v1.1.2 software (Jakobsson & Rosenberg, 2007) on the Structure Harvester website (Earl & vonHoldt, 2012). It was tested for K = 1 to K = 8 with 20 independent simulations at 60,000 samplings with a burn-in period of 10,000 first iterations. The final analysis with the resulting K (= 3) was conducted by 1,250,000 repetitions after a burn-in of first 500,000 replications. Individuals with at least 80% probability of membership in a cluster were considered to belong to that cluster. Individuals with probabilities of membership below 80% were interpreted as a hybrid genotype.

Different parameters of M. sinskiae populations including the haploid number of migrants (Nm), number of different alleles (N _a), number of effective alleles (N _e), number of private alleles (N _p), expected heterozygosity (H _e), unbiased expected heterozygosity (uH _e), and Shannon information Index (I) for each population were calculated using GenAlEx v6.503 (Peakall & Smouse, 2012) and POPGENE (Yeh et al., 1999). F-statistics employ inbreeding coefficients to describe the partitioning of genetic variation within and among populations and can be calculated at three different levels (F _IS, F _ST, F _IT). GenAlEx was also used to calculate the mean of haploid number of migrants or gene flow (Nm), genetic differentiation between subpopulations (F _ST), inbreeding coefficient of an individual relative to the subpopulation (F _IS), and inbreeding coefficient of an individual relative to the total population (F _IT and G _ST) for each primer pair (locus). G _ST is assumed to be an analogue of F _ST, and G _ST is equivalent to F _ST when there are only two alleles per locus and is the weighted average of F _ST for all alleles in case of multiple alleles per locus (Freeland, 2020). F _IS was estimated by dividing (H _e – H _o)/H _e (Pagnotta, 2018). The selfing rates (S) were calculated by dividing 2F _IS/ (1+ F _IS) in our populations (Burkil et al., 2017).

The polymorphism information content (PIC) was also calculated manually in spreadsheet as PIC = 1 – Σⁿ _i = 1p _i ², where p _i (i = 1, 2, 3, … I) are the frequencies of i ^th alleles for the given locus and “i” is the number of distinct alleles at a locus (Luo et al., 2019). Genetic distances and genetic diversity with pairwise test and Nei’s Genetic Distance (GD) were calculated using POPGENE and GenAlEx. Data were analyzed with POPTREE2 (Takezaki et al., 2009). A neighbor-joining (NJ) method based on Nei’s genetic distances using POPTREE2 with 1000 replicates of bootstrapping was used to illustrate the relationships between populations. The SplitsTree v4.15.1 software (Huson & Bryant, 2006) was used to calculate the genetic distances for a split neighbor net. The genetic difference within and among the studied populations (Table 1) as well as regions (Region 1, western Iran: populations ABD, KHW, KHR, SPC, SPD, and PLS; Region 2, southwestern Iran: populations FTH, BSN) was tested by AMOVA (Analysis of Molecular Variance) with 1000 permutations using GenAlEx. In order to evaluate the impact of genetic distance and geographical distance on population differentiation, a Mantel test (Mantel, 1967) was performed to correlate two matrixes of genetic distance based on F _ST and Nei’s genetic distance using GenAlEx.

The UPGMA dendrogram of SSR data with Dice similarity index showed that the 74 individuals sampled fell into three distinct groups (Fig. 2). Each of these groups represented the individuals of a different species. Similar results were also observed in the PCoA analysis of the same dataset (Fig. 3). In both analyses, all individuals of M. sinskiae were circumscribed as a single distinct group and therefore is in accordance with the proposal of Small & Brookes (1991) and the results of Zareei et al. (2020). Further analyses were performed on the 62 individuals of M. sinskiae (see Figs. 4–10 and Tables 4–6).

A total of 19 alleles were detected from the analysis of six loci (Table 3). The number of alleles per locus ranged from 2 to 6. The mean of Shannon Index (I) for eight populations was 0.458 (Table 4). The Shannon Index (I) values were higher than expected heterozygosity (H _e) values. Population ABD showed the maximum values of N _a (2.667), N _e (1.870), I (0.701), H _e (0.418) and uH _e (0.440), and population KHR showed the minimum values of N _a (1.333), N _e (1.308), I (0.224), H _o (0), H _e (0.160) and uH _e (0.178). The mean of observed heterozygosity (H _o), expected heterozygosity (H _e), unbiased expected heterozygosity (uH _e), number of different alleles (N _a), number of effective alleles (N _e) and fixation Index (F) for all eight populations were 0.040, 0.296, 0.317, 1.917, 1.535, and 0.849 respectively (Table 4). The F _IS value ranged from 0.582 to 1 in eight populations and population BSN showed the lowest level of F _IS (0.582), compared to the highest levels of F _IS (1) observed in populations PLS and KHR. Selfing rate (S) ranged from 0.735 to 1 with a mean value of 0.928 (Table 4). Nei’s genetic distance (GD) results (Table 5) showed the populations with more similar alleles having smaller genetic distances. Populations KHW and SPC both from Lurestan province were the most similar (Nei’s GD = 0.013). The largest amount of genetic distance was observed between populations BSN and ABD (Nei’s GD = 0.223).

Mean of polymorphism percentage was 77.08% (Table 4). Although populations ABD, KHW, SPC and FTH showed to be 100% polymorphic, population KHR showed the lowest polymorphism (33.3%) as well as lowest values of genetic parameters (Table 4). Our SSR markers displayed a high level of polymorphism, and this species seems to be a polymorphic plant. Both PIC and H _e (= gene diversity) values are measures of genetic diversity, although PIC values are not useful in linkage analyses when determining the inheritance between offspring and parental genotypes, and expected heterozygosity (H _e) is useful for haploid markers (Luo et al., 2019). The mean of following parameters was observed for six loci: haploid number of migrants or gene flow (Nm) = 1.999, coefficient of genetic differentiation among populations (F _ST) = 0.157, inbreeding coefficient of an individual relative to the subpopulation (F _IS) = 0.820, inbreeding coefficient of an individual relative to the total population (F _IT) = 0.847, and among-population genetic differentiation (G _ST) = 0.402. Locus MED-02 showed the highest value of F _{IS ,} G _ST and Nm and the highest value of F _ST and F _IT were observed in locus MED-04 (Table 3). The SSR markers showed the PIC values ranging from 0.196 (locus MED-01) to 0.999 (locus MED-05) (Table 3). Among the primers, locus MED-06 with six alleles had the highest number of polymorphic bands and locus MED-01, MED-02, MED-03 with two alleles had the lowest number of polymorphic bands (Table 3). The PIC value was used to measure the informativeness of primers. Having the highest PIC value (0.999), locus MED-05 showed higher polymorphism and had more impact in differentiation of individuals. The minimum amount of PIC value (0.196) was observed in the monomorphic locus MED-01, which was uniform in all individuals. All other primers were polymorphic (Table 3).

The true number of subpopulations (K = 3, Fig. 4) was obtained using the method of Evanno et al. (2005). As seen in Fig. 4, none of the populations seems to be uniformly consisting of a single cluster. All populations had individuals from different genetic clusters. Results of UPGMA clustering analysis showed that there was no major cluster formed by individuals solely from a single population or individuals with geographical proximity (Fig. 5). Different clusters included individuals from different populations with no geographical proximity. For example, individuals of population SPC were present in all four major clusters A1a, A1b, A2 and B. As seen in Figs. 6 and 7, individuals of different populations were scattered all over the PCoA plot and Neighbor-net network. In accordance with the results of Structure analysis (Fig. 4), UPGMA clustering dendrogram (Fig. 5), the PCoA plot (Fig. 6) and Neighbor-net network (Fig. 7), none of populations studied were unmixed. Despite being clearly circumscribed when analyzed alongside with M. rigidula and M. constricta, populations of M. sinskiae showed no grouping based on geographical proximity. In POPTREE software genetic distances measured for constructing phylogenetic trees of eight populations showed no relationship based on geographical proximity (Fig. 8). The AMOVA test was performed (Table 6) to study population differentiation and to estimate the percentage of intrapopulation and interpopulation genetic variation. Most of the genetic variation occurred among individuals (95%). The calculated genetic variation among populations was 5%, and the number for variation among regions was 0% (Table 6). The Mantel test indicated no meaningful correlation between genetic distance and geographical distribution (R = -0.019, P = 0.659).