DNA was extracted from small pieces of foot muscle from up to four specimens per lot by use of a QIAGEN DNA extraction kit for animal tissue following the standard procedure of the manual. Fragments of the mitochondrial 16S rRNA (16S) and cytochrome c oxidase subunit 1 (COI) genes were amplified by PCR using the primer pairs 16Scs1 (Chiba, 1999) and 16Sbd1 (Sutcharit et al., 2007) and L1490 and H2198 (Folmer et al., 1994), respectively. Reactions were performed with annealing temperatures / elongation times of 55°C / 90 s for 16S and 50°C / 60 s for COI, respectively. Both strands of PCR fragments were purified and cycle sequenced by use of the PCR primers. Electropherograms were corrected for misreads and forward and reverse strands were merged into one sequence file using CodonCode Aligner v. 3.6.1 (CodonCode Corporation, Dedham, MA). All sequences have been deposited in GenBank (Table S1). Sequence alignments were generated using MUSCLE as implemented in MEGA5 (Tamura et al., 2011). Sequence saturation was assessed for each mtDNA fragment by comparing values of the entropy-based index of substitution saturation (Iss) with its critical value (Iss.c). The test is implemented in DAMBE (Xia et al., 2003). Uncorrected pair-wise genetic distances were calculated using MEGA5 under the option ‘pair-wise deletion of gaps’. For phylogenetic analyses, 16S and COI sequence datasets were examined separately and concatenated into one partitioned dataset. Prior to the model-based phylogenetic analyses, the best-fit model of nucleotide substitution was identified for each gene partition separately by means of the Bayesian Information Criterion calculated with MrModeltest (Nylander, 2002). Partitioned models were applied in the Bayesian Inference (BI) analyses with parameters estimated from the data set. A Maximum Likelihood (ML) analysis was performed using RaxMLGui by applying the GTRGAMMA model (Silvestro and Michalak, 2010). One-hundred thorough ML bootstrap replicates, each with 10 runs, were performed to assess the branch support of the ML tree. Bayesian posterior probabilities of phylogenetic trees were estimated by running a 10,000,000 generations Metropolis-coupled Markov chain Monte Carlo (2 runs each with 4 chains, one of which was heated) as implemented by MrBayes vs. 3.2.1 (Ronquist and Huelsenbeck, 2003). A data partition was applied that allowed parameters to be estimated separately for each gene fragment and for each codon position of the COI gene. Sampling rate of the trees was 1,000 generations. Generations sampled before the chain reached stationary were discarded as burn-in. Stationarity was reached when the average standard deviation of split frequencies shown in MrBayes was less than 0.01 and the log likelihood of sampled trees reached a stationary distribution (Ronquist and Huelsenbeck, 2003).