For each dataset (each marker) a nucleotide substitution model was selected with MrModeltest 2.2 (Nylander, 2004). For the mitochondrial datasets the model was GTR+I+G, for H3 it was HKY+I+G. Bayesian analyses were done in MrBayes 3.2.1. (Ronquist and Huelsenbeck, 2003) hosted on the CIPRES Science Gateway (Miller, 2010). For each marker the analysis consisted of two simultaneous, four chain, MCMC runs (10 M generations). Trees were sampled every 1000 generations, the first 2500 trees were discarded as burnin (relburnin = yes, burninfrac = 0.25). Examination of the .p output files in Tracer v.1.5 (Rambaut and Drummond, 2007) showed stationarity was reached with proper effective sample sizes for all parameters (ESS > 200). Sumtrees (Sukumaran and Holder, 2010) was used to calculate 25% majority rule consensus trees. Subsequently the datasets for the individual markers were combined into two concatenated datasets (from hereon referred to as): the ‘stringent’ and ‘relaxed’ datasets. The stringent dataset (89 taxa) consisted of only protein coding genes (i.e. H3, COI and CytB) and had no missing data. The relaxed dataset (103 taxa) consisted of all markers (H3, COI, CytB and 16S); taxa for which only one marker was missing were also included. A partitioned analysis was set up in MrBayes (same version) for both datasets; for each partition the GTR+I+G model was selected using the above described procedure.