Migration‐tracking integrated phylogeography supports long‐distance dispersal‐driven divergence for a migratory bird species in the Japanese archipelago

Abstract Long‐distance dispersal (LDD) outside a species' breeding range contributes to genetic divergence. Previous phylogeographic studies of migratory bird species have not discriminated LDD from vicariant speciation in their diversification process. We conducted an integrative phylogeographic approach to test the LDD hypothesis, which predicts that a Japanese migratory bird subspecies diverged from a population in the coastal region of the East China Sea (CRECS) via LDD over the East China Sea (ECS). Haplotype networks of both mitochondrial and nuclear genes of its three subspecies were reconstructed to examine whether the Japanese subspecies of the Brown Shrike (Lanius cristatus superciliosus) diverged from an ancestral CRECS population. A species distribution model (SDM) for the Japanese subspecies was constructed using bioclimatic variables under the maximum entropy algorithm. It was projected backwards to the climate of the last glacial maximum (LGM) to infer the candidate source area of colonization. A migratory route of L. c. superciliosus, which possibly reflects a candidate past colonization route, was tracked by light‐level geolocators. Molecular phylogenetic networks suggest that the Japanese subspecies diverged from a population in the CRECS and maintained anciently diverged haplotypes. The SDM inferred that the emerged continental shelf of the ECS and the present CRECS were suitable breeding areas for the Japanese subspecies during the LGM. A major migratory route for L. c. superciliosus was inferred between the CRECS and the Japanese archipelago across the ECS. Our integrative approach supported the LDD hypothesis for divergence of the Japanese subspecies of the Brown Shrike. Shrinkage of the ECS may have been responsible for successful population establishment, due to a sufficient number of migrants overshooting to the Japanese archipelago from the CRECS. Our framework provides a new phylogeographic scenario for this region. Discriminating LDD and vicariance models helps improve our understanding of the phylogeographic histories of migratory species.

Haplotype networks of both mitochondrial and nuclear genes of its three subspecies were reconstructed to examine whether the Japanese subspecies of the Brown Shrike (Lanius cristatus superciliosus) diverged from an ancestral CRECS population. A species distribution model (SDM) for the Japanese subspecies was constructed using bioclimatic variables under the maximum entropy algorithm. It was projected backwards to the climate of the last glacial maximum (LGM) to infer the candidate source area of colonization. A migratory route of L. c. superciliosus, which possibly reflects a candidate past colonization route, was tracked by light-level geolocators. Molecular phylogenetic networks suggest that the Japanese subspecies diverged from a population in the CRECS and maintained anciently diverged haplotypes. The SDM inferred that the emerged continental shelf of the ECS and the present CRECS were suitable breeding areas for the Japanese subspecies during the LGM. A major migratory route for L. c. superciliosus was inferred between the CRECS and the Japanese archipelago across the ECS. Our integrative approach supported the LDD hypothesis for divergence of the Japanese subspecies of the Brown Shrike. Shrinkage of the ECS may have been responsible for successful population establishment, due to a sufficient number of migrants overshooting to the Japanese archipelago from the CRECS. Our framework provides a new phylogeographic scenario for this region. Discriminating LDD and vicariance models helps improve our understanding of the phylogeographic histories of migratory species.

| INTRODUC TI ON
How animal migratory behavior contributes to species distribution and speciation is a poorly studied area of biogeography (Greenberg & Marra, 2005;Winger et al., 2019). Evolutionary studies have recently detected greater opportunities for both population establishment in distant areas and genetic divergence in migratory than nonmigratory lineages (Lees & Gilroy, 2014;Rolland et al., 2014).
Vagrancy, long-distance dispersal (LDD) well outside a species' known range, is a candidate process for this evolutionary pattern.
Vagrancy occurs when organisms in movement are influenced by wind or other meteorological conditions, and it derives mostly from migratory populations that have, for instance in the case of birds, drifted past or overshot their anticipated destinations during their seasonal migration (Newton, 2008). Therefore, vagrancy allows successful establishment of an allopatric breeding population (O'Connor, 1986; but see Lees & Gilroy, 2014) and its subsequent genetic divergence. Meanwhile, because migration-related traits, such as a migration program forcing vagrants to continue migrating, may restrict establishment of a new population (Bensch, 1999), LDDdriven divergence in these species could be condition-dependent from the phylogeographical perspective. So far, for the divergence process of migratory species, many phylogeographic studies have assumed vicariance (e.g., Weir & Schluter, 2004;Zink et al., 2006), fragmentation of a once continuous breeding population by the appearance of a barrier to dispersal (Kropf et al., 2006). An alternative scenario, a disjunct breeding distribution, and divergence occur by LDD across a pre-existing barrier, has never been tested. This is because migratory bird species are continental and widespread in their distribution (Somveille et al., 2018), thus the phylogeographic processes of their speciation cannot be empirically assigned to process (Albert et al., 2017). To understand the diversification process of migratory species, it is crucial to propose a phylogeographic scenario that would explain how and when migratory species diverge via LDD.
Owing to its insular nature, well-documented paleogeography, and relatively high endemism of breeding migratory passerine lineages (Saitoh et al., 2015), the Japanese archipelago provides a suitable system within which to identify a biogeographic mode of speciation of a migratory lineage. During glacial periods of the Quaternary (approximately 2.7 million years ago to the present), due to lower sea levels, some parts of the Japanese archipelago were connected to the continent via land bridges, while other parts remained separated (Gallagher et al., 2015;Matsuzaki et al., 2019;Ohshima, 1990;Ota, 1998). Most terrestrial animal species currently occurring in Japan colonized the archipelago from the East Asian continent during glacial periods, and they subsequently diverged as Japanese endemics (McKay, 2012;Motokawa, 2017). Moreover, source continental populations could be either on the coastal region of the East China Sea (CRECS) or the northeast Asian continent (the Korean Peninsula, and northern China to Far East Russia), where distinct species or populations occupy their breeding ranges (e.g., Dong et al., 2015;Päckert et al., 2012;Saitoh et al., 2015;Zhao et al., 2017).
This system, predicting how the Japanese archipelago was colonized (from which regional population and along which route), will help indicate which of the two following hypotheses are more plausible. (a) Colonization across more than 500 km wide East China Sea (ECS) from the CRECS beyond the normal dispersal distance of many migratory bird species (Paradis et al., 1998) Although molecular phylogenetics is one promising solution for distinguishing between these two hypotheses, it is often the case that postdivergence gene flow biases inference of species relationships and past demography (Leaché et al., 2014). The Japanese archipelago would provide a rare opportunity to discriminate the LDD and vicariance hypotheses, yet phylogenetic analyses only may not be sufficient for scenario testing.
We integrated a species distribution model (SDM) and migration tracking using light-level geolocators with phylogeography.
SDMs can infer suitable breeding areas of species during the last glacial maximum (LGM), ca. 20,000 years ago (20 kya) from present species distributions (Elith et al., 2011). Although predicted distributions during the LGM do not directly reflect the ancestral area of the Japanese population at the time divergence occurred, any continental range that was also suitable during one of several glacial periods may be interpreted as an area from which colonization occurred more likely than others (Zink & Gardner, 2017). An avian migratory route over a sea barrier between the Japanese archipelago and the Asian continent may be interpreted as a candidate for the past colonization route to the Japanese archipelago for the following reasons. First, an existing migration route over a large geographic barrier may reflect LDD in the past when the barrier was smaller (Newton, 2008;Winger et al., 2019). Second, the migratory route also partly retraces past shifts in the breeding range (Alvarado et al., 2014;Newton, 2008;Ruegg et al., 2006;Winger et al., 2019). Therefore, in this system, both modes of past colonization by a breeding population may be reflected in the migratory route over a present sea barrier.

K E Y W O R D S
Brown Shrike, Japanese archipelago, light-level geolocator, long-distance dispersal, migratory route, species distribution modeling We studied the Brown Shrike (Lanius cristatus), a long-distance migratory passerine, as our first attempt to test these hypotheses because of the breeding distribution of its three subspecies. Lanius cristatus superciliosus breeds in the Japanese archipelago ( Figure 1), Lanius cristatus lucionensis breeds in the CRECS, and Lanius cristatus cristatus breeds in the northeast Asian continent and most of Sakhalin Island (Figure 2a; Lefranc & Worfolk, 1997). The ancestral breeding area of this species has been inferred to be in the Oriental region (Fuchs et al., 2019), including the CRECS region. We sequenced two mitochondrial genes and two nuclear introns. We constructed their individual haplotype networks and a dated mitochondrial phylogenetic tree to detect genetic relationships among the three subspecies and genetic signatures related with past population demography. Results of the three components can be differently predicted from the two hypotheses in this framework. (a) The LDD hypothesis predicts that the archipelagic population has been derived either from the ancestral Brown Shrike population or L. c. lucionensis in the CRECS. Signatures of old genetic bottlenecks and genetic divergence may be expected because LDD and colonization should occur with a relatively small population. For the archipelagic population, the SDM predicts that CRECS was suitable during the glacial period, and a migratory route across the ECS between the Japanese archipelago and the CRECS is considered to have been likely ( Figure 2b). (b) The vicariance hypothesis predicts that the archipelagic population genetically nests within, or is not well divergent from, the population on the northeast Asian continent, because the land bridge between the Korean Peninsula and the Japanese archipelago formed and disappeared repeatedly. For the archipelagic population, the SDM predicts that either current northern China or the Korean Peninsula were suitable during the glacial period, and a migratory route is expected to have existed over the Korean Peninsula from and to the Japanese archipelago ( Figure 2c). Although each method individually involves limitations and provides weak support, in combination together they improve confidence in scenario selection for the colonization process of the Japanese migratory population. From the more plausible scenario based on our approach, we discuss the conditions under that LDD could have contributed to the genetic divergence of this migratory species.

| Study species and field procedure
To track the migratory routes of L. c. superciliosus, we captured 25 adult male shrikes in Hokkaido and marked them with leg rings. We attached MK5740 geolocator tags (BioTrack) with leg-loop harness strings (Rappole & Tipton, 1991). We ensured that the total mass of the system (<1.1 g) did not exceed 5% of the birds' mass (Bridge et al., 2013). Eight birds (32%) returned to the breeding study site in 2018. Deployment of tags did not significantly change the return rate when compared with a previous report by Takagi (2003) (chisquare test for independence; χ 2 = 0.11, df = 1, p = 0.74). Four males lost their tags before recapture. One male marked with leg rings disappeared before breeding in 2018 and did not return to the site in 2019, so it remained uncaptured (the presence of this bird's geolocator was not confirmed).

| Laboratory procedure and haplotyping
We collected tissue and blood samples of L. c. cristatus from the Russian Far East; L. c. superciliosus from Hokkaido, and from Nagano Prefecture (Honshu Island); and L. c. lucionensis from the Ryukyu Islands. A total of 94 genetic samples were obtained. We also used database sequences of mitochondrial genes from GenBank and BOLD (Ratnasingham & Hebert, 2007) for coverage of its wider distribution ( Figure 2a). Since the Ryukyu Islands and the Philippines are occupied only by nonbreeding (i.e., migrating and wintering) L. c. lucionensis (Lefranc & Worfolk, 1997;Severinghaus, 1996), and as this subspecies is easily separable from the other two subspecies based on morphological characteristics (Lefranc & Worfolk, 1997), samples from those locations were used as representatives of L. c. lucionensis.
See Table S1.1 for details of 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 & Donnelly, 2003;Stephens et al., 2001). Haplotypes inferred with support values below 0.5 were not used in the following analyses.
See Appendix S1 for details about primers, and PCR reactions.

| Phylogenetic analyses
We constructed parsimony haplotype networks based on the median-joining algorithm with postprocessing (using NETWORK 5.0.0.3; Bandelt et al., 1999;Polzin & Daneschmand, 2003) for each gene region with outgroup species (the Gray-backed Shrike Lanius tephronotus and the Bull-headed Shrike Lanius bucephalus). For the nuclear networks, we also included sequences from other Asian Lanius species available on the database because genetic structure was less apparent in these regions than in mitochondrial DNA (Table S1.1).
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 (Kumar et al., 2016); 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 et al., 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 et al., 2018). We followed Aoki et al. (2018) for other settings.

| Species distribution modeling
We modeled the breeding range of L. c. superciliosus and projected it into the LGM. Data on occurrence in Japan were obtained from a summary of breeding localities reviewed by Kitazawa et al. (2020) (hereafter referred to as "review data"). Additionally, on 12 May 2019, we downloaded data from Global Biodiversity Information Facility (GBIF; https://www.gbif.org/, http://doi.org/10.15468/ F I G U R E 2 (a) DNA sampling localities. Circles and triangles represent samples from breeding and nonbreeding localities, respectively, and their colors indicate haplogroups inferred in the median-joining haplotype network of COI ( Figure 3a). Land colors indicate breeding ranges of presumed subspecies following Worfolk (2000). The striped region indicates a presumed hybrid zone between subspecies. Predicted colonization route (black arrows) and expected results from integrative approaches in (b) long-distance dispersal (LDD) and (c) vicariance hypotheses for the Brown Shrike. Predicted migratory routes (lines with filled circles), and suitable glacial regions (open circles on the continent) under the two hypotheses are shown. Regions indicated in gray represent terrestrial areas during the Early to the Middle Pleistocene, and dotted lines indicate the present coastlines dl.igj5uw). A polygon bounding the presumed distribution of L. c. superciliosus according to Lefranc and Worfolk (1997) was used to clip the data. Records were limited only to the breeding season (May-August) to avoid integration of migrating birds. After a few qualitychecking processes, final coordinates were thinned to include only one occurrence per 20 km square covering the Japanese archipelago ( Figure S1.1).
We constructed a model from the final coordinate data. Our purpose was to infer the continental area that might have been suitable for L. c. superciliosus as its breeding habitat during the LGM. Therefore, any ecological differences between the current archipelago and the continental area need to be integrated in the model (Elith et al., 2011 We used the R package "ENMeval" v. 0.3.0 (Muscarella et al., 2014) to execute the maximum entropy (MaxEnt) algorithm (Phillips et al., 2006) across a range of several settings.
Combinations of three feature classes (linear, a combination of linear and quadratic, and hinge) and regularization multipliers ranging from 0.5 to 3 (increments of 0.5) were tested. The "maxent.jar" option was used as the model algorithm. From the geographical modeling domain defined above, 10,000 background points were sampled randomly. The "block" method was used to partition data for cross-validation since it assesses model transferability better than other methods (Muscarella et al., 2014). The hinge feature, with a regularization multiplier of 1.5, was selected by the lowest AICc. The present distribution was predicted, and then it was projected to the LGM environmental data using cloglog trans- environmental similarity surface (MESS) was also calculated and extrapolation was restricted to regions indicated by negative values from this analysis (see Appendix S1 for the detailed procedures on the SDM analysis).
We determined sunset and sunrise transitions for specific dates by light intensity thresholds of 0.5. We incorporated a statistical modeling approach, called the "group model" in a Bayesian framework implemented in "SGAT." This model refines the estimated locations by means of three components: the twilight model, a movement model, and spatial masks. The twilight model incorporates the deviation between determined and actual twilight times in calibration periods of each tag. We used twilight data for 26 days after geolocator deployment as the calibration period.
We fitted gamma distribution as an error distribution of detection of twilight. To define a movement model and spatial masks separately for stationary and movement periods, the ChangeLight function in "GeoLight" was applied to discriminate between these periods from the twilight data. A probability threshold for discrimination was set to 0.9. For the movement model, we set a flight speed prior with a gamma distribution for periods of movement (shape = 13, rate = 0.3), determined allowing for the average air speed of 12.9 m/s for the Red-backed Shrike Lanius collurio (Bruderer & Boldt, 2001). We conducted a sensitivity analysis for the gamma distribution and found that priors did not greatly influence our interpretation of the results (Table S1.2). In the group model, twilight times grouped for a stationary period are treated as records at a single location, thus no movement model was defined. For a probability mask, birds must cross the sea between the Japanese archipelago and the continent and must stopover on land. Therefore, a spatial mask for stationary periods was designed only to include terrestrial land while both sea and land areas were allowed for movement periods. We modified our analysis slightly for bird #V5604-025 because its geolocator stopped recording during the autumn equinox period (see Appendix S1 for details).
We interpolated latitudinal estimations around the equinox periods by setting a "tol" value between 0.08 and 0.1. As burn-in and fine-tuning, we ran 120,000 chains with three iterations, and the final paths were determined by the last 15,000 samples with three iterations. Regarding large uncertainties involved in migration estimations by geolocators (Lisovski et al., 2019), we standardized the posterior probabilities of location estimates on each raster cell, calculated their means among the three individuals, then presented them as a probability distribution in a region of interest where we expected birds to retrace the route of past colonization. We also inferred median paths of the three migratory tracks. Because birds spend more time at stopover sites than at an instantaneous point along a migratory path, the former are indicated with relatively higher probabilities than the latter.

| Phylogenetic analyses
An MJ haplotype network of COI was constructed using 107 sequences from the species' range ( Figure 3a). This network inferred three genetic groups, which were supported as distinct clades in the BI phylogenetic tree based on the concatenated mitochondrial sequences ( Figure S1.2). The archipelagic, northern, and southern clades are thought to reflect the three extant subspecies (L. c. superciliosus, L. c. cristatus, and L. c. lucionensis, respectively; Figure 2a). In this network, a split between the southern and archipelagic clades was observed with nine mutational gaps. The archipelagic clade was inferred to be more ancestral than the northern clade. The major COI haplotype of the northern clade (H1) was shared among samples from a wide range of localities in Far East Russia (Table S1) The BI tree indicated that the basal lineage of the Brown Shrike was the southern clade ( Figure S1.2). The first split occurred be-

| Species distribution modeling
The present breeding distribution of L. c. superciliosus was modeled based on presence-only data and seven noncorrelated bioclimatic variables using the MaxEnt algorithm ( Figure 4). Minimum training presence omission rates for test localities were 1.9% and the model was evaluated to estimate suitable habitat without problematic overfitting. The present prediction of suitable areas (Figure 4a) reflected the taxon's known breeding range according to Lefranc and Worfolk (1997).
Projection of the distribution model to three climatic models of the LGM (Figure 4b-

| Migratory route analysis
Under the Bayesian framework, autumn and spring migration routes were inferred for three individual L. c. superciliosus (Figures 5, S1.5).
Full annual data were retrieved from two tags, and similarity was observed between their autumn and spring routes. One tag stopped recording on the continent due to cell failure 3 months after deployment. Despite uncertainties involved with location estimates, higher mean probabilities for migratory locations were observed among the three individuals both in the western part of the Japanese archipelago and along the CRECS to northern Taiwan. The probability map and the median migratory tracks suggest that birds migrated across the ECS between these two regions and did not detour over the Korean Peninsula. The data from the two successfully retrieved geolocators revealed that both shrikes had wintered in Indonesia, in the Sunda Islands, specifically on Java and Flores.

| D ISCUSS I ON
The archipelagic population of the Brown Shrike was inferred as a genetically divergent clade by the mitochondrial tree and networks.
In the COI network, a phylogenetic split between the southern and archipelagic clades was inferred, while the northern clade was inferred to be a later-derived clade (Figure 3a). Moreover, the archipelagic population preserved single haplotypes of the two nuclear introns, which were found outside the major Brown Shrike clusters and as divergent as those found from different Asian Lanius species (Figure 3c,d). These results agree with the LDD hypothesis that the Japanese population was once established from the ancestral population of this species and then isolated from the continental  Integration of the SDM and migration tracking were supportive of the LDD hypothesis. Suitable glacial breeding habitat for L. c. superciliosus was estimated to have existed around the CRECS during the LGM by the SDM (Figure 4), and its current migratory route was observed between the CRECS and the Japanese archipelago across the ECS ( Figure 5). Within our scenario-testing framework, our results indicate that the LDD-driven divergence of a Japanese migratory subspecies across the ECS was more likely than the vicariance hypothesis. Signatures of old genetic divergence of the archipelagic haplotypes in both the COI and nuclear networks seem to reflect the same phylogeographic event (Zink & Barrowclough, 2008), possibly the ancient LDD-driven colonization and divergence in the Japanese archipelago (Figure 6a,b). Therefore, less divergent nuclear haplo- during the LGM, and no migratory route was indicated along them, unlike the routes reported for other species (Koike et al., 2016;Shiu et al., 2006). Furthermore, since northern Taiwan was connected to the CRECS during the glacial period (Figure 6), the inferred migration aligns with our prediction.

| Hypothetical phylogeographic scenario of the Brown Shrike
The ancestral breeding area of the Brown Shrike is believed to be in the Oriental region (Fuchs et al., 2019) Ota, 1998). Our SDM analysis is in line with existing paleopalynological records indicating that the emergent continental shelf of the ECS was covered with grassland and sparse forest (Xu et al., 2010), which would have provided suitable breeding habitat for the Brown Shrike (Lefranc & Worfolk, 1997). Present vagrancy records of many migratory species in the Japanese archipelago during spring are probably due to overshooting (birds traveling too far because of wind and anticyclonic conditions; Newton, 2008) from the CRECS (Senzaki et al., 2019). A sufficient number of migrants toward the CRECS, therefore, may have accidentally traveled too far and arrived in the Japanese archipelago, resulting in a rare population establishment event by LDD (Figure 6a Meanwhile, the evolutionary origins of the two continental populations remain inconclusive. Genetically, the northern clade was found to be most distant from the southern clade but close to the archipelagic clade in our COI network, possibly suggesting divergence of the northern continental population from the archipelagic population by reverse colonization (Nishiumi & Kim, 2015).

| Justification for the use of a migratory route as a proxy for the past colonization pathway
Our data on the migratory route of the Brown Shrike across the ECS is far too limited for population-level inference, thus the universality of the route for the entire archipelagic subspecies remains unclear. Careful interpretation of our tracking data is required because individual variation in migratory routes is recognized (Stanley et al., 2012). Migratory routes inferred from the analysis also included large uncertainties ( Figure 5) because of the methodological limitations of light-level geolocators (Lisovski et al., 2019). However, Lefranc and Worfolk (1997) stated that, based on direct observation during migration, many L. c. superciliosus fly at least 700 km across the sea from Japan to the coast of China's Jiangsu Province or the Zhoushan Islands. Our geolocator data concur with their statement.
Although further deployment of geolocators is needed, their statement guarantees that the focal sea crossing route is at least one of the major routes for the archipelagic population.
We have taken the recent view that idiosyncratic routes of obligate long-distance migratory species are highly persistent traits throughout their evolutionary histories (Winger et al., 2019) in order to interpret part of our tracking data as the route of past colonization. However, the dominant view is to consider migratory routes as highly labile traits that are easily subjected to natural selection (Pulido, 2007;Winger et al., 2019). Thus, one could argue that the observed migratory route across the ECS is a result of optimization to the present environmental conditions and does not reflect past colonization (Alerstam, 2011;Berthold et al., 1992;Sutherland, 1998;Zink & Gardner, 2017).
For the following two reasons, we believe Winger et al.'s (2019) and our interpretation to be more plausible than the counterview. First, our system does not include any region that experienced severe glaciation. Earlier studies have examined species that inhabit regions previously covered by ice sheets during the LGM (Haché et al., 2017;Milá et al., 2006;Ruegg et al., 2006;Sokolovskis et al., 2018). Reverting from being migratory to becoming sedentary residents was probably a major response to glaciation, which was inferred by SDM analyses (Zink & Gardner, 2017). This reversion implies that migration is labile throughout a species' evolutionary history and migratory routes are likely to trace the postLGM expansion to the present breeding range. By contrast, we have shown that suitable breeding habitats for the archipelagic population remained in the Japanese archipelago even during the LGM (Figure 4b-d). This implies that the focal migratory route traces a distribution change before the LGM, supporting our view that idiosyncratic routes over a large barrier have been conserved. Second, similarity between autumn and spring migration across the ECS ( Figure 5) may support the conservatism of the focal route. If a migratory route is optimized for the present environment, then different selection pressures upon passages in different seasons should independently shape the route (Stanley et al., 2012;Tøttrup et al., 2012). Generally, loop migration, in which spring and autumn migration routes are divergent, is thought to have evolved by retaining one ancestral route retracing colonization while deriving the other to adjust to current environmental conditions (Newton, 2008). Wind is an important determinant (Alerstam, 2001). This loop migration phenomenon has been shown for a population of the Oriental Honey Buzzard Pernis ptilorhynchus breeding in Japan. Its autumn migration involves crossing the ECS, whereas its spring migration route detours around the Korean Peninsula to cross a narrow sea channel to enter the Japanese archipelago (Nourani et al., 2016). It has been shown that a route over the Tsushima Strait was specifically suitable for spring wind conditions when wind directions and strengths over the ECS were unstable (Nourani et al., 2016;Yamaguchi et al., 2012). These conditions would be unfavorable also for passerines such as shrikes (Alerstam, 2001;Tøttrup et al., 2017).

| Conclusion
Previous studies have illustrated a divergence scenario for migratory species that assumed vicariance in generating isolated populations (Hewitt, 2004;Weir & Schluter, 2004). While this scenario has provided suitable explanations for many cases, establishment of a new population from a distant area across a barrier as a result of LDD is likely to be their alternative explanation, given that migratory bird species have heightened opportunities for LDD and genetic divergence than nonmigratory lineages (Lees & Gilroy, 2014;Rolland et al., 2014). However, traditional frameworks based on the continental distribution of most migratory species have perhaps hindered the opportunities to test it (Warren et al., 2015). Our integrative approach using a Japanese migratory bird species has provided a suitable system within which to discriminate the mode of genetic divergence under a hypothesis-testing framework. We have successfully provided a hypothetical scenario for how LDD and paleoenvironment have both contributed to the divergence of a migratory species. We have also shown that discriminating LDD from vicariance in our framework benefited a more comprehensive understanding of diversification mechanisms of migratory species in this region. We believe that consideration of LDD is also important for continental regions where most migratory species have diversified. Attention has recently been paid to LDD and colonization as an important phenomenon for distribution and speciation of continental migratory species (Rushing et al., 2015;Winkler et al., 2017).
Further verification of the usability of our integrative approach among several species will facilitate our understanding of the relative importance of LDD to vicariance in the speciation of migratory birds and other animals. from the Japanese Society for the Promotion of Science (JSPS).

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
Genetic data is available at DDBJ (accession numbers are listed in