Laboratory procedure and haplotyping
We collected tissue and blood samples from the Russian Far East forL. c. cristatus , from Hokkaido and Nagano prefectures forL. c. superciliosus and from the Ryukyus for L. c. lucionensis . A total of 95 genetic samples was obtained. We also used database sequences of mitochondrial genes from GenBank and BOLD (Ratnasingham & Hebert, 2007) for a coverage of its wider distribution range (Figure 1c). Since the Ryukyus and the Philippines are occupied only by non-breeding (i.e. migrating and wintering) L. c. lucionensis (Lefranc & Worfolk, 1997; Severinghaus, 1996), they were used as representatives of L. c. lucionensis . See Table S1.1 for details about samples and sequences.
Total genomic DNA was extracted from blood and tissue samples using DNeasy Blood & Tissue Kit (Qiagen). Two mitochondrial genes (complete cytochrome b (cytb ), and partial cytochrome oxidase c subunit I (COI)) and two autosomal introns (myoglobin intron-2 (MB), transforming growth factor beta 2 intron-5 (TGFb2)) were amplified by polymerase chain reaction (PCR). The autosomal introns were sequenced only for a subset of the samples (Table S1.1). Automated sequencing was run on either an ABI7130 or ABI3130 Genetic Analyzer (ABI). Sequences were aligned using ProSeq v. 3.5 software (Filatov, 2009). To determine haplotypes of nuclear introns with multiple heterozygous sites, we used PHASE v. 2.1.1 (Stephens et al., 2001; Stephens & Donnelly, 2003). Haplotypes inferred with support values below 0.5 were not used in the following analyses. See Appendix S1 for details about primers, PCR reactions, and PHASE analysis.
Phylogenetic analyses
We constructed a median-joining (MJ) haplotype network (using NETWORK 5.0.0.3; Bandelt, Forster, & Rohl, 1999) using partial COI (521 bp) from 108 samples. We also constructed a multi-locus network of four loci (cytb , COI, MB and TGFb2) under the NeighborNet algorithm (using Splits Tree 4.14.4; Huson & Bryant, 2006). We only included individuals with sequences of all four loci and calculated uncorrected pairwise (p ) distances among sequences in each locus using MEGA v. 7 (Kumar, Stecher, & Tamura, 2016). Then, we combined the fourp -distance matrices to produce a standardized genetic distance matrix among individual samples using POFAD v. 1.03 (Joly & Bruneau, 2006), which was used as an input file for the NeighborNet analysis.
We reconstructed mitochondrial phylogeny by using the two mitochondrial sequences. A dataset consisting of unique haplotypes of concatenated sequences of cytb and COI was prepared for tree reconstruction by the Bayesian inference (BI) method using BEAST v. 2.4.8 (Bouckaert et al., 2014). In the BEAST analysis, two genes were independently assigned to different partitions with the best substitution models determined by the lowest Bayesian Information Criterion on MEGA v. 7; HKY + G was selected for each gene (α = 0.10 for cytb and 0.19 for COI). Application of the standard avian clock rate for cytb , 0.0105/lineage/million years was justified since shrikes are passerines (Nabholz, Lanfear, & Fuchs, 2016; Weir & Schluter, 2008). Markov chain Monte Carlo (MCMC) sampling proceeded for 50 million states with 100 thousand pre-burn-in chains. Convergence for each parameter was checked through Tracer v. 1.7.1 software (Rambaut, Drummond, Xie, Baele, & Suchard, 2018). We followed Aoki et al. (2018) for other settings.