Data scoring and analysis

All the samples were amplified and produced clear, reproducible bands (S2 Fig) which were scored as present (1) or absent (0), using the BioNumerics (v7.5; Applied Maths, Belgium) software. A binary matrix (1/0) was generated containing 195 loci which was used for subsequent data analysis for diversity and population structure (S4 Table).

The informativeness of ISSR markers was evaluated using the polymorphic information content (PIC) resolving power (RP), mean resolving power (MRP), marker index (MI) and Shannon’s index (H′). PIC is the probability in detecting polymorphism by a primer PIC  =  1-Σ(pi)², where Pi is the frequency of the ith allele [27]; RP is the ability of each primer to detect level of variation between individuals and was calculated according to Prevost and Wilkinson [28]: 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 the polymorphic information content and effective multiplex ratio (MI  =  PIC × EMR [29]. Shannon’s index (H′) was calculated by the formula H′ = -pilnpi [30].

Genetic diversity within each population was estimated through percentage of polymorphic loci (P), mean effective number of alleles (Ne), mean expected heterozygosity (H_E), mean Shannon’s Information Index (I), Nei’s gene diversity (h'), and Nei’s pairwise genetic distances by using the GenAlEx 6.5 software [31,32]. The POPGene v 1.32 [33] software was used to calculate the overall Genetic differentiation (G_ST = H_T-H_S/H_T). The G_ST was corrected according to Hedrick [34] and Jost [35] to obtain the estimates for G’_ST and D, respectively. We also performed an analysis of molecular variance (AMOVA) to examine the distribution of genetic variability within and among populations using the Arlequin software (version 3.5) [36] and estimated an overall F_ST (Fixation index). Based on the F_ST the Gene flow (N_M = 1(1/F_ST-1)/4 was also calculated [37]. The pairwise F_ST from the Arlequin software were used to derive the Slatkins relationship (F_ST/1-F_ST) to create an isolation by distance plot through the Mantel’s test employing 10,000 permutations [36].

Population structure was determined by the neighbour-net split decomposition network generated by SPLITS Tree v4.6 and bootstrapping runs of 1000 replicates [38], a PCoA analysis of pairwise F_ST between the populations [36] and Bayesian model-based clustering analysis using the software program STRUCTURE 2.3.3 [39]. The admixture model and correlated allele frequencies were used for each run with a burn-in period of 1,000 and 100,000 Markov chain Monte Carlo replications. The optimal K value, which indicates the number of genetically distinct clusters in the data, was determined from 10 replicate runs for each value of K [40]. The value of ΔK was based on the change in the log probability of the data between successive K values. Structure Harvester version 6.0 [41] was used to calculate parameters described by Evanno et al. [40]. The usepopinfo option of STRUCTURE was employed after K was determined by the Evanno’s method to estimate Q (proportion of membership of each pre-defined population in each of the clusters).