Mitochondrial phylogenetic tree and divergence dating
Complete mitochondrial genomes were annotated with MitoFinder (48) and aligned with MAFFT l-ins-i. We first investigated phylogenetic relationships among Glaucopsyche individuals by analyzing the entire mitochondrial genome. We used IQ-TREE2 (49) to select the best fitting substitution model for each partition and merge similar partitions and build a maximum likelihood tree and assessed support with 1000 ultrafast bootstrap replicates.
To infer a time-calibrated phylogeny, we selected one individual each of Xerces Blue and Silvery Blue in addition to 12 related Polyommatinae. We extracted the sequences for protein-coding genes, aligned each with the codon-aware aligner MACSE, and concatenated all. We used BEAST2 (50) with the bModelTest package to perform phylogenetic site model averaging for each of the merged partitions. Because there is no accepted molecular clock rate for butterflies and no fossils to apply in this part of the phylogeny, we used two strategies to apply time constraints to the analysis. First, we used two published molecular clock rates for the mitochondrial COX1 gene (1.5% divergence/Ma estimated for various invertebrates (51), and the ‘standard’ insect mitochondrial clock 2.3% divergence/Ma (51). We applied a strict clock with a normal prior set up to span the 1.5-2.3% range with the 95% HPD interval (mean=1.9%, sigma=0.00119). Second, we borrowed the age of the most recent common ancestor of our sampled taxa from fossil-calibrated analyses across butterflies, which has been estimated to ~33 Ma(7, 52). We fixed the root age to 33 Ma and allowed the remaining node ages to be estimated using a strict clock. Analyses were run twice from different starting seeds for 10 million MCMC generations and trees were sampled every 1000 generations. Runs were checked for convergence with Tracer and all effective sample size (ESS) values were >200. Runs were combined with LogCombiner, after removing the first 30% of topologies as burn-in, and a maximum credibility tree was generated with TreeAnnotator (50). Phylogenetic analyses were performed on the National Life Science Supercomputing Center - Computerome 2.0 (www.computerome.dk).
Xerces Blue and Silvery Blue population histories
We used the Pairwise Sequentially Markovian Coalescent (PSMC) model to explore the demographic history of both species (Xerces Blue and Silvery Blue). We obtained a consensus fastq sequence of the mappable fraction of the genome for each autosomal chromosome (total of 22 chromosomes ofG. alexis assembly). Only positions with a depth of coverage above 4X and below 15X were kept. Posteriorly a PSMC was built using the following parameters: -N25 -t15 -r5 -p ”28*2+3+5”. We used 1 year for the generation time and a mutation rate of 1.9x10-9, estimated in Heliconius melpomene (53). Considering that calling consensus sequences from low coverage samples (< 10x) can underestimate heterozygous sites (54), and given the different coverage between samples, we corrected by False Negative Rate the samples with coverage lower than the coverage of L005 (for Xerces Blue) and L013 (for Silvery Blue), as recommended by the developers of the software, so that all samples are comparable with each other. However, since in our dataset we do not reach a coverage > 20x, we acknowledge that we are not capturing the whole diversity and thus our PSMC might infer lower historical effective population sizes.