Microsatellite and mtDNA data analysis
Allele size range and private alleles were determined using CONVERT version 1.31 (Glaubitz, 2004). Genetic diversity parameters such as number of alleles per locus (NA), number of private alleles (PA), observed (HO ) and expected (HE ) heterozygosities were calculated by using POPGENE version 1.31 (Yeh, Yang, Boyle, Ye, & Mao, 1997). Polymorphic information content (PIC) and inbreeding coefficient (FIS) were calculated with Microsatellite Tollkit software (Park, 2001) and FSTAT version 2.9.3.2 (Goudet, 1995). Genetic differentiation tests between populations (pairwise FST) and molecular variance analyses (AMOVA) were conducted using ARLEQUIN (Excoffier, Laval, & Schneider, 2005). The genetic structures among populations and individuals were investigated using STRUCTURE version 2.3.1 (Pritchard, Matthew, & Donnelly, 2000).We estimated the optimalK value using the ΔK statistic as described by Evanno, Regnaut, & Goudet (2005) in STRUCTURE HARVESTER (Earl & vonHoldt, 2012). In addition, GENETIX v 4.05 (Belkhir, Borsa, Chikhi, Raufaste, & Bonhomme, 2004) program was applied for Factorial Correspondence Analysis.
Sequences of mtDNA gene regions were aligned with CLUSTAL W2 (Larkin et al., 2007). mtDNA sequences were checked with reference sequences (Gen Bank accession number for COI and cyt b , respectively : GU085204.1, JQ820820.1) by using BLAST 2.2.20 (Zhang, Schwartz, Wagner, & Miller, 2000). DNA sequence alignment was performed by uploading all sequence information to the MEGA6 program (Tamura, Stecher, Peterson, Filipski, & Kumar, 2013) for both gene regions. Sequences for detection of SNP regions were uploaded to DnaSP v.5 (Librado & Rozas, 2009) software and each bee was then assigned to a mitochondrial haplotype. We compared mitochondrial haplotype frequencies in each population using ARLEQUIN 2.00 (Schneider, Roessli, & Excoffier, 2000). As a result of statistical analyzes intra and inter-group genetic variations (F statistics) and genetic distances of 14 B. terrestris populations were calculated. The genetic distance information of the populations was uploaded to SplitsTree 4.0 (Huson & Bryant, 2006) software, to conduct phylogenetic analysis and to construct Neigbour-Joining (NJ) tree. In order to determine the evolutionary ancestral origin of each individual in B. terrestris populations, the Maximum Likelihood (ML) method based on the Tamura-Nei model was used (Tamura & Nei, 1993). Trees were created in MEGA X program using 633 bp nucleotide sequences for 66 samples for COI gene region and 428 bp nucleotide sequences for 68 samples for cyt b region (Kumar, Stecher, Li, Knyaz, & Tamura, 2018). B. lucorum , used as an out group (OG) in the present study, were retrieved from Gen- Bank (Accession number: JQ843492.1) to construct both NJ trees.