Phylogenetic tree reconstruction
For the mitochondrial COI dataset, the alignment of 323 sequences was used to construct a Bayesian phylogenetic tree using the software Beast2 (Bouckaert et al, 2014), with out-group sequences of N. purpurata , N. meanyi , N. vandykei , N. praedicta , and N. metallica downloaded from GenBank (accession numbers: MN346054, KU641246, KU641248, JQ711382, KU875540). The model HKY+G was selected as the best substitution model based on model tests conducted with the software MEGA v7 (Kumar, Stecher, & Tamura, 2016). A clock rate for the COI gene was set to 0.0113 per site per million years based on the best available COI mutation rate estimate for beetles, in this case related carabid beetles (Andújar, Serrano, & Gómez-Zurita, 2012). Ten million steps of a Markov chain Monte Carlo (MCMC) were run to estimate model parameters and trees, with samples drawn every 1000 steps. Two separate runs were conducted. The posterior samples were assessed for convergence using TRACER v1.7 (Rambaut, Drummond, Xie, Baele, & Suchard, 2018) and each parameter value was assessed to ensure a high effective sample size (>200). Additionally, a statistical parsimony cladistic analysis, or TCS network, was reconstructed using Popart-1.7 (Clement, Snell, Walker, Posada, & Crandall, 2002; Leigh & Bryant, 2015). Based on the divisions in the TCS network, haplotypes were assigned to northern or southern haplogroups.
A population tree based on genome-wide SNP data (12,498 SNPs,fully-filtered dataset ) was constructed using maximum likelihood in the program IQtree version 1.6.12 (Nguyen, Schmidt, Von Haeseler, & Minh, 2015). The substitution model GTR+ASC was selected and the branch support was estimated based on 1,000 bootstrap replicates. The output tree was then visualized and colored using FigTree v1.4.1 (Rambaut, 2012). Additionally, we inferred the species tree using the program SNAPP (Bryant, Bouckaert, Felsenstein, Rosenberg, & RoyChoudhury, 2012), implemented in BEAST2 version 2.6.2. Thefully-filtered dataset was used, but due to computational constraints, one random individual was selected from each of eight populations. Individuals were chosen from populations with relatively pure ancestry coefficients to avoid violating assumptions in the species tree approach, namely the effect of gene flow from other populations (see the result of sNMF analysis; Bryant et al., 2012, Stoltz et al., 2019). SNAPP was run for 30 million Markov chain Monte Carlo (MCMC) generations, with the custom program SNAPPER employed to run SNAPP at faster speed (Stoltz et al., 2019). Two independent runs were conducted and the output trees were visualized after 10% burnin using DensiTree version 2.01 (Bouckaert 2010).