Phylogenetic analyses and population structure
DNA sequences were aligned using Clustal W (Thompson, Higgins, & Gibson, 1994) with default options in Geneious Pro v8.1.9 (Biomatters, Auckland, New Zealand). Haplotypes were identified with a Bayesian coalescent-based program performed in DnaSP v5.10.1 (Librado & Rozas, 2009). All mtDNA cox1 and ITS1 haplotypes were deposited in GenBank (cox1 : MN728370-MN728538; ITS1: MN728333-MN728363). Phylogenetic relationships among haplotypes were inferred using maximum likelihood (ML) analysis and Bayesian inference (BI) implemented in RAxML v8.2.9 (Stamatakis, 2014) and MrBayes v3.2.6 (Huelsenbeck & Ronquist, 2001) using nucleotide sequences of M. keenae (GenBank Accession nos. for cox1 : MN728366 – MN728369; GenBank Accession nos. for ITS1: MN728364 – MN728365) as an outgroup. The best fit models for the two genes were determined using the Akaike information criterion (AIC) in jModelTest 2.1.10 (number of substitution schemes = 11; Darriba, Taboada, Doallo, & Posada, 2012): GTR+I for cox1 and TVM+G (0.47) for ITS1. The ML analysis was performed using 1,000 bootstrap replicates. The BI analysis was conducted with two independent runs, each with four simultaneous Monte Carlo Markov chain (MCMC) runs for 1 x 106 generations. Trees were sampled in every 100 generations and posterior probability was estimated after removing an initial 25% of generations (reaching a stationarity of the chains) as burn-in. Haplotype networks were constructed using PopArt v1.7 (Leigh & Bryant, 2015) software via the TCS method (Clement, Posada, & Crandall, 2000).
Molecular diversity indices such as haplotype diversity (h ), nucleotide diversity (π ), and mean number of pairwise differences (k ) of each population were estimated with ARLEQUIN v3.5 (Excoffier & Lischer, 2010). Pairwise genetic divergence between populations was evaluated using the fixation index (F st) (Excoffier, Smouse, & Quattro, 1992). A hierarchical analysis of molecular variance (AMOVA) was conducted to estimate population structure among geographical groups, which were determined by a spatial analysis of molecular variance performed with SAMOVA v2 (Dupanloup, Schneider, & Excoffier, 2002). These analyses were carried out with 1,000 permutations after sequential Bonferroni adjustments using ARLEQUIN.