Genetic diversity analyses using iPBS data

All bands that could be reliably read across all individuals were scored as either present (1) or absent (0) across genotypes and treated as single dominant loci. Based on the binary matrix obtained (Table S1), the following genetic parameters were estimated using GenAlEx 6.5 (Peakall & Smouse, 2006, 2012): total number of bands per population (N_B), percentage of polymorphic bands (P), Shannon’s information index (I) and expected heterozygosity (H_E).

Two methods were used to investigate the genetic structure of the samples. The first approach was the Bayesian model-based clustering method implemented in STRUCTURE ver. 2.3.4. (Pritchard, Stephens & Donnelly, 2000). The model assigns individual multilocus genotypes probabilistically to a user-defined number of clusters (K), achieving linkage equilibrium within clusters (Pritchard, Stephens & Donnelly, 2000). We conducted 10 replicate runs for each K, ranging from 1 to 10 (Fig. 2A). Each run consisted of a burn-in of 500,000 iterations, followed by data collection of over 2,000,000 iterations. The analysis using admixture model was conducted without any prior information on the original population. To determine the optimal number of clusters, an ad hoc statistic ΔK was used (Evanno, Regnaut & Goudet, 2005). The ΔK was evaluated in Structure Harvester ver. 0.6.94 (Earl & Vonholdt, 2012). The second method was a principal coordinates analysis (PCoA), based on the matrix of Euclidean distances between individuals from all analyzed populations, performed in PAST software (Hammer, Harper & Ryan, 2001).

(A) The values of the second-order rate of change of L(K), ΔK, of data between successive K values. (B) The population structure bar plots generated at K = 2.

Analysis of molecular variance (AMOVA) was performed with Arlequin 3.5. For this analysis, the iPBS data was treated as haplotypic, comprising of a combination of alleles at one or several loci (Excoffier, Laval & Schneider, 2005). The significance of the fixation indices was tested using a non-parametric permutation approach, the method implemented in Arlequin 3.5 (Excoffier, Smouse & Quattro, 1992; Excoffier, Laval & Schneider, 2005). Moreover, Tajima’s D, Fu’s F_S neutrality test, and the mismatch distribution and demographic processes affecting populations were estimated using the same software. Bottleneck ver. 1.2.02 (Cornuet & Luikart, 1996) software was used to investigate recent effective population size reductions based on allele data frequencies (Cornuet & Luikart, 1996; Piry, Luikart & Cornuet, 1999) for each population. In populations that have experienced a recent reduction in their effective population size, the H_E becomes larger than the heterozygosity expected at mutation-drift equilibrium. In order to study such effect using dominant markers, the infinite allele model (IAM) was used to test the mutation-drift vs bottleneck hypothesis (Tero et al., 2003). The significance of potential bottleneck was estimated using a sign test, standardized differences test and one-tailed Wilcoxon sign rank test for heterozygosity excess.