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.