Response of an Afro-Palearctic bird migrant to glaciation cycles.

Migration allows animals to exploit spatially separated and seasonally available resources at a continental to global scale. However, responding to global climatic changes might prove challenging, especially for long-distance intercontinental migrants. During glacial periods, when conditions became too harsh for breeding in the north, avian migrants have been hypothesized to retract their distribution to reside within small refugial areas. Here, we present data showing that an Afro-Palearctic migrant continued seasonal migration, largely within Africa, during previous glacial-interglacial cycles with no obvious impact on population size. Using individual migratory track data to hindcast monthly bioclimatic habitat availability maps through the last 120,000 y, we show altered seasonal use of suitable areas through time. Independently derived effective population sizes indicate a growing population through the last 40,000 y. We conclude that the migratory lifestyle enabled adaptation to shifting climate conditions. This indicates that populations of resource-tracking, long-distance migratory species could expand successfully during warming periods in the past, which could also be the case under future climate scenarios.

Seasonality resulting from the tilted Earth orbiting the Sun is an important factor shaping the distribution of life on Earth (1). As seasonality increases toward the poles, year-round availability of resources declines (1, 2). Nonmobile species need to sustain survival throughout the season with the lowest resource availability, limiting population size outside of the peak of resources in the most resource-rich months during which breeding normally occur (3, 4). Migration allows animals to exploit the seasonal variation in resources (5–8). Every year, as the northern latitude winter approaches, billions of birds leave their breeding grounds, traveling thousands of kilometers to warmer, lower latitudes (5, 9). Comparable migrations are found in a wide range of taxa from small insects (10) to the largest of the whales (11), but seasonal migrations have been lost in much of the terrestrial megafauna, probably owing to human impacts (12).

Migratory behavior is common across birds, although not all species migrate (13). The spatiotemporal schedules can be surprisingly fine tuned so that migratory species adjust the timing and migratory routes so to optimize the harvest of resources as they become available at continental scales (7). Because of the complex spatiotemporal schedules, migratory species are generally considered vulnerable (14). However, given the ability to track changes in habitat availability across the globe, migration likely represents an advantageous adaptation to exploit the changing distribution of seasonal resources over time, which would potentially be reflected in population size (15). The trait is often considered as plastic and changeable over short, decadal time scales (16, 17). Current migratory patterns derive from successful adjustments dating back to glacial periods when higher-latitude breeding grounds were unsuitable for breeding (18). During these glacial stages, populations of migrants retracted to the south and then expanded back to the north during interglacial periods (13, 15, 19, 20). Zink and Gardner (21) found that most currently long-distance migrants lacked suitable breeding area in North America during glacial maxima and suggested that a “migratory switch” would turn migrations off. That is, migratory populations would become sedentary to areas that were wintering grounds during warmer interglacial periods as indicated by the phylogeography of sedentary and migratory populations of chipping sparrows Spizella passerina (20). In contrast, Ponti et al. (22) found no support for a switch in long-distance migrants between Eurasia and Africa. While previous research (21–23) has explored the maintenance of migrations at breeding–wintering resolution, our study assesses migratory behavior at a higher temporal resolution (monthly estimates) and combines it with corresponding estimates of past population size inferred from whole-genome data, an independent line of evidence not combined with spatial prediction in previous studies.

Inferring past migration patterns and population dynamics offer insights into future responses of migrant birds to accelerating climate change. Here, we combine independent lines of evidence from spatial models and genomics to explore the past migrations and shed light on this long-lasting debate on migratory behavior through shifting climates. We used direct tracking data of individual red-backed shrikes Lanius collurio (24), an Afro-Palearctic songbird with a complex, spatiotemporal, annual schedule that allows individuals to optimally exploit seasonal surplus of vegetation greenness across continents (7). Schedules consisting of alternating stationary and traveling periods are normally not captured in Species Distribution Models (SDMs), particularly not those inferring distributions of migrants in the past (15, 25, 26) or the future (27, 28). To capture these complex schedules, we used dynamic seasonal modeling, which takes into account the actual conditions experienced throughout the annual cycle (29) and paleoclimatic and vegetation conditions (30) at a fine spatiotemporal resolution to hindcast the potential seasonal distributions back to the previous interglacial, the Eemian. The habitat availability was based on modeled climate suitability and potential build-up of vegetation. Our paleoclimate data consisted of simulated monthly climatic conditions and photosynthetic activity spanning 120,000 y B.P. Moreover, we use genomics data to estimate effective population size back in time (31) to infer whether past species distribution and migration patterns have left imprints on population size (15). Together, estimated past spatiotemporal schedules and population dynamics shed light on previously unknown past migration patterns and population consequences for a long-distance migratory species.

Results

We find that a complex migration system in red-backed shrikes has likely been present since before the Last Glacial Maxima (LGM) 21,000 B.P. (Fig. 1), at least throughout the last 120,000 B.P. (Fig. 2). Our results suggest that breeding took place in Europe as today during warm periods but in northern Africa during glacial periods, including the LGM (Fig. 1). However, refugial populations might have existed during the main breeding season (June and July) in the south-central region of Europe during the Last Glacial Period (LGP). Although the migration distance between the breeding and nonbreeding grounds was reduced in the LGP, migration likely persisted in tropical regions, given that we find no overlap of the estimated latitude of all months from cells with maximum projected suitability at any time step (Fig. 1).

Fig. 1.

Maximum predicted monthly suitability (top 5% cell occupancy) for migratory red-backed shrikes from the LGM (21,000 B.P.) until present. (A) Annual migration of red-backed shrikes breeding in South Scandinavia as revealed by geolocator tracking; positions during stationary periods are indicated (color represent months). (B) Distance between breeding and wintering latitudes of maximum predicted suitability from the LGM to present; wintering latitude was in Africa throughout, and breeding latitude during LGM was in Africa (blue) but in Europe (red) from 14,000 B.P. (C) Predicted area of occupancy during breeding (July, orange) and wintering (January, green) at the LGM (Left) and present (Right). Red-backed shrikes do not currently breed in Africa, and cells with predicted area of presence in Africa during breeding are crossed. (D) Predicted area of occupancy during stationary periods across the annual cycle from LGM (Left) to present (Right); lines show the estimated mean latitude of maximum suitability cells for each month.

Fig. 2.

From the Eemian (120,000 B.P.) to the present (0 B.P.), temporal alignment of (A) changes in global surface temperature during the last 1 million y (in Celsius compared to the present indicated by dotted, vertical line), (B) estimated latitudinal migration distance (kilometers), (C) predicted area of suitable breeding habitat (number of cells with suitability above monthly threshold; dotted line indicates the present condition), and (D) effective population size estimates (BSP) for Lanius clades 1 and 2.

None of the regions across the study area were found to be suitable all year around any time step (Fig. 1). Breeding (June) and wintering (January) areas were well separated throughout (Fig. 1 ; overlap occurring only for cells with suitability less than 0.4). Estimates of mean longitude across the annual cycle suggest that migrations may have occurred across a broader longitudinal span of Africa in the past, shifting toward the eastern part of Africa during the Holocene but also during the previous interglacial period. Alternatively, suitable breeding areas were exclusively in the southern Palearctic during cold periods as shown by breeding-only hindcasting (). These areas are small and fragmented during cold periods (). In combination with the lack of concurrent low-population size (Fig. 2), this indicates that during cold periods, the areas in Africa predicted by the full annual model () were used. Breeding and wintering areas also were well separated when breeding and wintering were modeled separately ().

The seasonal migration dynamics were overall similar in thrush nightingales () with separated breeding and wintering () throughout. Breeding areas of thrush nightingales showed less-pronounced shifts between Palearctic and Africa with main breeding in Africa indicated only during the coldest periods ().

We find overall consilience between the estimated changes in suitable habitat and the population sizes during the warm periods, since the Eemian (Fig. 2 and ). Generally, the area of predicted presence was larger during the warm periods (the Eemian and the Holocene) than during the LGP (Fig. 3). Since the LGM, suitable conditions expanded gradually in Europe during the breeding season (Fig. 3). Genomic data indicate that the effective population size of red-backed shrikes has been consistently large during the recent warm period when most habitat was available primarily in Europe (Figs. 2 and 4 and ), with a growing trend in population size dating back 40,000 y when less habitat was available.

Fig. 3.

Summed area of all cells with predicted presence from binary projection (threshold maximizing sensitivity and specificity for each month), presented as mean of all months (Left) and quarterly means (Right)  from 120,000 B.P. to present (0 B.P.). Shaded area indicates 95% CI.

Fig. 4.

Effective population size estimates back in time. BSPs of the two Lanius clades based on complete mitochondrial genomes. The solid lines represent the median estimates of the effective population size, with the 95% high posterior density interval denoted by the dashed lines (note the logarithmic scale on the y-axis).

Discussion

Our modeling imply that populations of red-backed shrikes performed seasonal long-distance migration in variable forms through glacial cycles. Using a temporally explicit SDM based on tracks of individual shrikes, hindcasted habitat suitability reveals a seasonal migratory pattern similar to the one currently observed but with breeding and nonbreeding stopover areas displaced during glaciation when the birds remained in Africa throughout the year, moving between what is now the Sahel to Southern and Central Africa and back again. Relying on seasonal movements, suitable habitat was available to the birds year-round throughout the climate cycles of the past 120,000 y with larger areas of habitat becoming available during the warm periods. Independent genomic data on effective population size back in time indicate that both cold periods with shorter migration and warm periods with the large areas of suitable habitat requiring a longer migration sustained a large population size.

Migration tracks of several songbirds have recently revealed complex, seasonal migration schedules with little variation among individuals (7). Our temporally explicit SDMs reproduced the observed spatiotemporal progression across the annual cycle—essential to justify the use in hindcasting suitable habitat availability. The fact that suitable habitat was available in the past and showed a qualitative relationship to the estimated effective population size back in time confirms the general assumption underlying hindcasting (25, 32) that the niche has been stable over time as a reasonable assumption.

Given that the patterns found here are representative of glaciation cycles during the Pleistocene, this suggests that the complex migration dynamics observed today are old, potentially preceding the glaciation cycles of the Pleistocene. In general, trends for effective population size back in time varies considerably both quantitatively and qualitatively among populations (33). However, the population structure of the species studied here is relatively well known with two distinct clades that largely mix across the entire distribution (34). Our estimated effective population sizes reflected this pattern indicating that we have captured the relevant population dynamics. Furthermore, the lack of bottlenecks is largely consistent with the amount of hindcasted available habitat.

We find support for migration persisting through time (22, 23) with no switch from resident to migratory populations (21). During the LGM, migration seems to have occurred within Africa in the same seasonal spatiotemporal framework as currently observed, likely as a loop migration pattern, with the breeding grounds being in either south-central Europe or northern Africa, assuming that the main breeding season was similar to current times (June and July). As the ice retracted gradually after the LGM in Europe, larger areas of suitable habitat became available, and migration developed southward during the austral summer (January to March) and northwards during summer in the northern hemisphere (June to August). The breeding-only, hindcasted suitable habitat is fragmented and reduced during the GLM (21), in most cases more than 50% (22). Even if enough habitat was available to sustain a population, a substantial decrease in population size would be expected. We find that seasonality in Africa has shifted during periods with glacial expansions, creating the opportunity for intra-African migration. The availability of additional breeding habitat within Africa likely allows for sustaining the growing population size during the cold period despite a decrease in extent of habitat. A simulation study (23) indicated that this could be more general across species, but our study reveals how seasonal modeling in migratory songbirds with complex migration patterns is necessary to disentangle the full range of migration dynamics in the past.

Understanding the flexibility of migratory species to adapt to changes in the environment is essential for predicting how they will respond to future changes in climate and land use (14). We find a high predicted presence in the present. The increase in the suitable area seems to level off during the Holocene. In the period 10,000 B.P. to 4,000 B.P., northern East Africa went from a cold and dry climate to warm and humid conditions (35), while temperatures started rising in Europe (36). Given that past climatic conditions during this time are representative of future climatic scenarios, we might expect that smaller areas of suitable habitat will be available for red-backed shrikes in the future. Furthermore, a recent study found that red-backed shrikes may face difficulties in tracking resource peaks throughout the annual cycle during future climate scenarios (7).

Nevertheless, following the changes in optimal spatiotemporal migration over time would likely be possible given that migratory behavior is an adaptable trait that can change in just a few generations. This includes whether migratory behavior is expressed or suppressed as well as alteration of migration direction and distance (16, 17). For example, a novel migration pattern in blackcaps, which involved changes in both distance and direction, has apparently evolved within a decade (16). Nevertheless, migration patterns are surprisingly conserved in the longest-distance migratory species tracked, and we expect a certain degree of hard-wiring to be adaptive given the complex schedules involved. Shifts from warm interglacials with wide distributions to restricted distributions during glacial maxima occurred gradually, allowing gradual adaption. In contrast, transitions from cold to warm periods occurred fast. Nevertheless, the resulting expanded range would allow a gradual change of migration behavior to exploit the newly available habitat.

Seasonal migration has a considerable impact on the distribution of biodiversity across Earth both at short and long timescales. The large-scale changes in distributions suggested here in migrants involve recurrent, long-term dispersal over long distances both during breeding and nonbreeding seasons, resulting in a larger long-term distribution than current patterns reveal. This could potentially provide mechanisms for speciation in combination with suppression of migration behavior for populations that became isolated as has been suggested in Madanga ruficollis, a forest-adapted endemic bird of the island Buru, Indonesia, which is closely related to the tree pipits of the genus Anthus (37), and it could be more widespread than previously thought (38).

We conclude that the migratory lifestyle of red-backed shrikes is a successful evolutionary adaptation, which allows exploiting seasonal vegetation dynamics, sustaining sizeable populations through time. As such, we suggest that at least for the longest-distance migrating species tracking ephemeral resources, labile migration patterns over time are probably widespread and that more generally migratory lifestyles can be beneficially adapted to various climatic conditions experienced during glacial cycles in the past and expect this to also be the case in the future. The framework provided here outlines an approach for a detailed understanding of past migration and movement as a response to seasonality. Combined with increasing knowledge from population genomics, this would improve our understanding of the processes that may have shaped movements and seasonal distributions of migratory birds in the past.

Materials and Methods

Modeling Spatiotemporal Distribution.

We used the r-package sdm version 1.0–32 (39) to model the spatiotemporal distribution of red-backed shrikes during the annual migration cycle from their European breeding grounds to their nonbreeding grounds in Sub-Saharan Africa and back.

Presence Data.

We derived spatiotemporal data published by Pedersen et al. (40) on annual migrations from archival light-level geolocators of red-backed shrikes (n = 24, 31 tracks) tracked from three breeding sites in southern Scandinavia (data available from ref. 41). Loggers were deployed in 2009 through 2015 and retrieved in 2010 through 2017. Geolocators measure and store light levels, which, upon retrieval of the tag, can be converted into latitude and longitude estimates from derived day length and local noon and midnight (42). Due to the relatively low spatial resolution of longitudes and latitudes derived from geolocator data, estimated as 50 ± 34 km and 143 ± 62 km (mean ± 95 CI), respectively (43), we restricted our sampling of presence points to individual mean estimates of stationary periods. Average individual positions from stationary sites were sampled on a monthly scale. During breeding periods (June and July) when birds were known to be at the breeding site, we included known breeding area positions for each individual. This resulted in a total of 432 positions across the annual cycle.

Pseudoabsence Data.

As the extent of the background area affects estimates of SDMs (44), the extent of the study area should encompass the occurrence of the species but not exceed it (45). Due to the wide distribution of breeding and nonbreeding ranges of the species (46), we defined the extent of the study area to cover two continents (−40° to 75° latitude, −20° to 60° longitude), using a geographical latitude–longitude projection on the global reference system WGS84 datum. Offshore islands (including Yuzhny, Greenland, Iceland Island, the Canary Islands, Cape Verde, and Madagascar) were removed from the background as they are unlikely to be part of the species' distribution. As the tracking data are, by default, presence only and most model algorithms need presence–absence data, we sampled pseudoabsences from the background (44, 47). For each month, we sampled five times the number of unique presences chosen randomly across the study area, avoiding spatiotemporal overlaps with the presence data. This resulted in 131 to 186 pseudoabsences per month sampled during the months August through April, while in May, June, and July, when individuals are at their known breeding sites, 55, 15, and 20 pseudoabsences were sampled, respectively.

Environmental Data.

We used a global dataset of downscaled and bias-corrected monthly precipitation, temperature, and net primary productivity (NPP) for the past 120,000 y (48, 49). Climate data in this dataset were derived by combining simulations of the HadCM3 and HadAM3H climate models with global observational records, while NPP estimates were obtained by applying the curated climate data to the Biome4 vegetation model (50). The data match empirical reconstructions of temperature, precipitation, and vegetation well (48). From the NPP data, we estimated surplus NPP for each cell as the relative measure of NPP for each month (scale: −1 to 1). The three variables (temperature, precipitation, and surplus NPP) were chosen based on a biological perspective of the drivers of songbird migration in terms of resource availability (7) and with respect to the availability of simulated data of the past.

We hindcasted back to 120,000 B.P., and for every 1,000-y interval until 21,000 B.P. and 2,000-y interval until 120,000 B.P., monthly averages of 30 y were calculated. For each unique presence and pseudoabsence position, we extracted spatiotemporal information on the environmental predictors.

Model.

The annual spatiotemporal model was run with temperature, precipitation, surplus NPP, their cubic polynomial, and interaction terms of all predictor variables. We used the following algorithms: generalized linear model (glm), boosted regression tree (brt), random forest (rf), and Maxent. We also used a bootstrapping procedure to resample data into 30 replications. We then fitted the models for each algorithm and replication, resulting in 120 models (4 algorithms × 30 replications = 120 models). As the model output may rely on the choice of algorithm (51), we created an ensemble of models (52) that used a weighted mean to combine the predicted values generated by the 120 models. The performance values of the model, measured by AUC (area under the receiver operating characteristic curve), were used as the weights in the ensemble procedure. We projected the model output on each month of the current climate data (0 B.P.) and then hindcasted the model for each month in the 72 time steps back to 120,000 B.P.

Analyses.

The model performance of each algorithm was assessed using the threshold-independent AUC statistic (53). Our model framework resulted in a good fit, with mean AUC values above 0.9 except for June (0.53) and August (0.88) (glm: 0.88, brt: 0.90, rf: 0.97, and maxent: 0.90). We plotted maps of predicted suitability of the time periods during specific periods of particular interest in terms of climatic conditions (35): 0 B.P., 10,000 B.P., 13,000 B.P., 21,000 B.P., 70,000 B.P., and 120,000 B.P. () using the r-package rasterVis version 0.43 (54). Red-backed shrikes optimally track seasonal vegetation surplus throughout their annual cycle (7). To investigate how migration has likely changed over time, we extracted cells with the highest-estimated monthly suitability to match the current spatiotemporal migration patterns. Monthly maps of the 5% highest suitability recreated the current spatiotemporal patterns well (Fig. 1). For the 5% of cells with the highest suitability, we estimated the mean latitude and longitude for each month and time step (patterns were similar based on the 10% of cells with highest suitability; ). If migration had ceased during the past glaciation cycles, we would expect to find an overlap in the estimated monthly latitudes throughout the annual cycle at the LGM. Furthermore, we assessed whether overlap in predicted presence occurred between breeding (June) and wintering (January) depending on suitability threshold. To investigate whether the amount of the available suitable area had changed over time, we created a binary projection using a threshold based on the True Skill Statistic, which maximizes sensitivity and specificity, for each month to define occurrence and absence. We then summed up the area of predicted occurrence for each month and time step and on average for each time step. All analyses described above were run in the statistical software R version 3.4.2 (55).

Modeling.

The full annual model is based on the Scandinavian population for which we also have population size data. Additionally, we modeled breeding and wintering range at present and during the LGM for individuals across a wide latitudinal and longitudinal span of the breeding range of red-backed shrikes from Pedersen et al. (40). Suitable areas were modeled and hindcasted using the annual data but separately for breeding based on positions in June and July only and wintering based on positions in December through February only.

The full annual model predicted to a large extent occurrence throughout the species’ current range in Europe. This appears reasonable given that our annual model is based on full annual spatiotemporal habitat use and that all populations at present are migratory with similar complex migration patterns (40). The extent of present potential breeding and wintering ranges in the Palearctic based on breeding-only and wintering-only data were rather similar to the full annual model. Because the breeding area in the annual model was similar to the species current range, we find it overall justified to generalize our results at the population level to the species level.

The full annual model additionally indicated suitable breeding habitat in Africa. The overprediction of areas of suitable habitat at present in Africa could be viewed as a challenge for the accuracy of our hindcasting of suitable habitat. The distribution-wide, breeding-only model predicted very limited and fragmented suitable breeding habitat during the LGM. The apparently large population size during the LGM indicates that the shrikes likely used the additional potential breeding areas available in Africa at that time. For these reasons, we believe that the scenario based on the full annual modeling is the most likely. Regardless, both models predicted that migration continued through cold periods.

Hindcasting Migration in Another Long-Distance Migrant.

To investigate whether the seasonal migration patterns inferred back in time were unique for red-backed shrikes among the longest-distance migrants, we used the same framework for hindcasting seasonal distributions of thrush nightingales Luscinia luscinia, for which we also had tracking data (7). Compared to red-backed shrikes, the model of thrush nightingales tended to overpredict the breeding area to a larger extent, but it still reproduced the current spatiotemporal migration pattern.

Historical Effective Population Size.

We estimated the historical effective population size (Ne) using both autosomes and mitochondrial DNA of the red-backed shrikes. We sequenced and assembled the genome for a male individual (scaffold N50 of 24 Mb; ) and sequenced 14 other samples () to an individual mean coverage ≥ 23 ×. The L. collurio individuals were all from eastern Denmark, 13 from Grib Skov, 1 from Møn, and 1 from Gedser. They were caught from May to July during 2009 through 2018.

After assembling the mitochondrial genome of the sample used to build the reference genome, we recovered the consensus sequences for all individuals and extracted fragments of mitochondrial fragments cytochrome oxidase I (COI) and cytochrome B (CYTB). The mitochondrial gene phylogenies equally divided the individuals sampled in Denmark between the two clades previously described by Pârâu et al. (34) ().

We reconstructed mitogenomic demographic histories independently in two phylogeographic groups (Lanius clade 1 and 2) using the Bayesian Skyline Plot (BSP) (56). BSP analyses were performed in BEAST version 1.8.4 using a strict molecular clock with a substitution rate of 2 × 10−8 substitutions/site/y (57) and HKY85 substitution model with gamma distribution. The chain length was established at 1 × 108 iterations, with samples drawn every 10,000 Markov chain Monte Carlo steps. After a discarded burn-in of 25%, BSPs were reconstructed and visualized in a plot with Tracer version 1.6.0 (58). The timing of the estimated increase in Ne is uncertain as it depends on the substitution rate. When using substation rates between 1 × 10−8 and 3 × 10−8, the beginning of the increase ranged from 100,000 B.P. to 30,000 B.P. Effective population sizes estimated using pairwise sequentially Markovian coalescent (PSMC) curves covering mainly before the period of hindcasting are presented in .

DNA Extraction, Library Preparation, and Sequencing.

Tissue samples were extracted using the MagAttract HMW DNA kit (Qiagen) according to the manufacturer’s indications. Blood samples were extracted using the same kit with some modifications. Prior to lysis incubation, samples were mixed by pulse-vortexing three times at the highest speed setting. The final elution volume in the AE buffer was 150 μL. Extractions were quantified using a Qubit dsDNA High Sensitivity assay (Life Technologies).

One 10× genomics linked-read library was constructed (59) for a male individual (sample ID 143784) and sequenced on BGI-seq 500 platform in 150 paired-end modes for the de novo assembly. In total, 153.0 Gb (125.9× coverage) of raw reads were produced. DNA samples from the other 14 individuals were sequenced on BGI-seq 500 platform in 100 paired-end modes for the resequencing. For each sample, one pair-end library with insert size 200 to 400 bp was constructed. On average, 48.5 Gb (39.9× coverage) of raw reads were produced for each individual. A series of filtering steps were applied for these resequencing reads prior to the downstream analyses using SOAPfilter 2 package (version 2.2):

Removing reads with more than 10% of N bases;

Removing reads with more than 40% low-quality bases (Phred score <=10);

Removing reads with undersize insert size; and

Filtering out the PCR duplicates (if read1 and read2 of the same paired-end reads were identical, we considered them duplicates).

Genome Assembly, Repeat Annotation, and Mapping.

All raw reads were introduced into Supernova (59) software (version 2.0.1) under the pseudohap style to assemble the reference genome of L. collurio. After removing the scaffolds with “N” up to 80% from this preliminary assembly, GapCloser (60) (Version 1.12) was used to close the intrascaffold gaps. The statistics of the final L. collurio assembly are listed in . Genome completeness was measured with BUSCO (61).

For the repeat annotation, Tandem Repeats Finder version 4.07b (62) was used to identify the tandem repeats, and both the homology-based and the de novo approaches were applied to identify the transposable elements. The homology-based repeat annotation was done by RepeatMasker version open-4.0.7 (63) (with parameters “-nolow -no_is -norna -engine ncbi -parallel 1”) at the DNA level based on the Repbase library (version: 20170127). The de novo repeat annotation was done by RepeatModeler open-1–0-8 (64) with default parameters to first build a de novo repeat library for each assembly. Further, the de novo repeat library with RepeatMasker-open-4.0.6 (63) were further used to predict repeats for the final assembly. All of the above results were then merged into a union set.

The mitochondrial genome was assembled with NOVOPlasty (65) using a COI nucleotide sequence as a seed (National Center for Biotechnology Information [NCBI] ID: MH938014.1) and the mitochondrial genome from Lanius schach (NCBI ID: NC_030604.1) (66) as a reference. This resulted in an assembly that is 16,820 bp long (average insert size: 239 bp; average organelle coverage: 135×). Scaffolds corresponding to the mitochondria were identified via BLAST and replaced by the full mitochondrial sequence in the reference genome.

The reads from the 14 resequenced individuals were mapped to the reference genome with Burrows-Wheeler Aligner BWA-MEM (version 0.7.12-r1039). Reads showing a mapping hit were further filtered for mapping quality >25. PCR duplicates were removed with Picard MarkDuplicates (version 1.95; http://picard.sourceforge.net) and local realignment around indels was done with GATK (67). Sequencing error was assessed in ANGDS (68) using the option “-doError 0,” which estimates the joint type-specific error rates for multiple samples using the minor allele frequency as prior (described in ref. 69). Only scaffolds larger than 10 Mb were used in downstream analyses (), and regions annotated as repeats by RepeatMasker (63) were masked with BEDtools (70). Scaffolds corresponding to the avian sex chromosome Z were identified as those having half the average coverage in females ().

Mitochondrial Phylogeny.

Consensus sequences for the resequenced individuals were obtained by choosing the most common base per position (-doFasta 2 in ANGSD, 68). Fragments from CYTB and COI were downloaded from the NCBI (GenBank IDs in ) (34, 66, 71–76). Trees were inferred using the Neighbor-Joining method (77) with the Jukes–Cantor distance (78) and 500 bootstrap runs (79). The analysis included 37 COI nucleotide sequences and 39 CYTB nucleotide sequences (GenBank IDs are shown for sequences obtained from NCBI; the asterisk denotes samples from the present study). All ambiguous positions were removed for each sequence pair. There were a total of 606 positions in the final COI dataset and 567 positions in the CYTB dataset. Evolutionary analyses were conducted in Molecular Evolutionary Genetics Analysis (MEGA) X (80).