Ethics statement
Animals’ care was in accordance with institutional guidelines. Ethical permit was issued by Malmö-Lund djurförsöksetiska nämnd 5.8.18-00848/2018.
Field work
We carried out the field work in Sweden during four breeding seasons (2018–2021). Adult male willow warblers were captured in their breeding territories using mist nets and playback of a song. From each bird, we collected the innermost primary feather from the right wing. From the birds that returned with a logger we also collected ~20 μl of blood from the brachial wing vein. The blood was stored in SET buffer (0.015 M NaCl, 0.05 M Tris, 0.001 M of EDTA, pH 8.0) at room temperature until deposited for permanent storage at −20 °C. We deployed Migrate Technology Ltd geolocators (Intigeo-W30Z11-DIP 12 × 5 × 4 mm, 0.32 g) and used a nylon string to mount them on birds with the “leg-loop” harness method as outlined in our previous work24. The mass of the logger relative to that of the bird was on average 3.3% (range 2.7–3.8%).
The tagged birds were ringed with a numbered aluminum ring, and two, colored plastic rings for later identification in the field. In total, we tagged 466 males (349 in 2018 and 117 in 2020) at breeding territories. During the first tagging season (2018), birds were trapped at 17 locations (average 22 birds per site; range 7–30) distributed across Sweden (Fig. S1). Three of the sites were in southern Sweden to document migration routes of allopatric trochilus and three sites were located above the Arctic circle to record migratory routes of allopatric acredula, whereas the remaining (239) loggers were spread over 11 sites located in the migratory divide. Given the observed densities and distribution of hybrids after analyzing returning birds in 2019, we deployed 117 more loggers at one single site (63.439°N, 14.831°E) in 2020. We successfully retrieved tracks from 57 birds tagged in 2019 and 16 from birds tagged in 2021. In search for birds with loggers, we checked circa 3000 willow warbler males and covered an area of at least 0.5 km radius around each site the year after tagging.
Geolocator data treatment
The R package GeoLight (version 2.0)25 was used to extract and analyze locations from raw geolocator data. All twilight events were obtained with light threshold of 3 lux. The most extreme outliers were trimmed with “loessFilter” function and a K value of 3. We used GeoLight’s function “getElevation” for estimating the sun elevation angle for the breeding period: these sets of locations were used to infer the positions for autumn departure direction. In addition, we carried out a “Hill-Ekström” calibration for the longest stationary winter site during the period before the spring equinox. Winter calibration produced location sets that better reflected the winter coordinates of the main winter site in sub-Saharan Africa26. We reduced some of the inherent geolocation “noise” by applying cantered 5-day rolling means to the coordinates. The equinox periods were visually identified by inspecting standard deviations in latitude. Latitudes from equinox periods were omitted (on average autumn equinox obscured data for 45 days (range 25–68). For the main winter site, we used the longest period at which bird stayed stationary and from which in all cases begun the spring migration (mean = 118, SD = 23 days). Timing of autumn departure was estimated by manual inspection of longitudes and latitudes plotted in time series. To estimate at which longitude the birds crossed the Mediterranean, we extracted the longitude when birds crossed latitude 35 N° (Mediterranean crossing longitude). For 29 birds, it was possible to directly extract the longitude at crossing latitude 35 N°. For the rest of the cases, the birds had not reached latitude 35 N° before the latitude was obscured by the equinox, we calculated the mean longitude of 10 days from the onset of fall equinox as a measure of the Mediterranean crossing. This measurement correlated highly with the winter longitude (r = 0.78, p = 2.8 × 10−16). To control for the birds relative breeding site longitude, we extracted the departure direction (1°–360°) relative from the tagging site to the location where the birds crossed the Mediterranean (departure direction). The departure data was of circular type (measured in 360°), however the variance did not span more than 180° degrees (range 151°–224°). Therefore, we proceeded with analyses using linear statistics. Geographic distances and departure direction were calculated using R package “geosphere” (version 1.5-10). Complete set of positions of each individual bird with equinoxes excluded is presented in Supplementary Data 1.
Laboratory work and molecular data extraction
We extracted DNA from blood samples following the ammonium acetate protocol16. Genotyping for divergent regions on chromosome 1 (InvP-Ch1) and chromosome 5 (InvP-Ch5) was done using a qPCR SNP assay16, which is based on one informative SNP per region (SNP 65 for chromosome 1 and SNP 285 for chromosome 5). Probes and primers were produced by Thermo Fisher Scientific and were designed using the online Custom TaqMan® Assay Design tool (Table S4). We used Bio-Rad CFX96™ Real-time PCR system (Bio-Rad Laboratories, CA, USA) and the universal Fast-two-steps protocol: 95 °C, 15 min—40*(95 °C, 10 s–60 °C, 30 s, plate read. Both regions contain inversion polymorphisms that restrict recombination between subspecies-specific haplotypes and contain nearly all the SNPs separating the two subspecies13. For each region, we scored genotypes as either “Tro” (homozygous for trochilus haplotypes), “Acr” (homozygous for acredula haplotypes) or “Het” (heterozygous). The method that we used to assess the presence of MARB-a is based on a qPCR assay that quantifies the copy number of a novel TE (previously known as AFLP-WW212) that has expanded in acredula. The quantification of repeats by this method has been shown to be highly repeatable (R2 = 0.88) when comparing estimates obtained from DNA in blood and feathers15. We used the forward (5′-CCTTGCATACTTCTATTTCTCCC-3′) and reverse (5′-CATAGGACAGACATTGTTGAGG-3′) primers developed by Caballero-López et al.15 to amplify the TE motif. For reference of a single copy region we used the primers SFRS3F and SFRS3R27. We diluted DNA to 1 ng/μl−1 and used a Bio-Rad CFX96™ Real-time PCR system (Bio-Rad Laboratories, CA, USA) with SYBR-green-based detection. Total reaction volume was 25 μl of which 4 μl of DNA, 12.5 μl of SuperMix, 0.1 μl ROX, 1 μl of primer (forward and reverse), and 6.4 μl of double distilled H2O. We ran quantifications of the single copy gene and the TE variant found on MARB-a on separate plates with the following settings: 50 °C for 2 min as initial incubation, 95 °C for 2 min X 43 (94 °C for 30 s [55.3 °C SFRS3 and 55.5 °C for TE, 30 s] and 72 °C for 45 s). Each sample was run in duplicate and together with a two-fold serial standard dilution (2.5–7.8 × 10−2 ng). Allopatric trochilus have 0–6 copies whereas allopatric acredula have 8–45 copies15; a bimodal distribution was also confirmed in this new data set (Fig. S2). Accordingly, for the present analyses, we split the data in two groups: birds with ≤6 TE copies and birds with >7, translating into absence or presence of MARB-a, with the former assumed to be homozygous for the absence of MARB-a and the latter heterozygous or homozygous for the presence of MARB-a. Data from two investigated willow warbler families suggest a Mendelian inheritance pattern and provide support for our interpretation of how TE copy numbers reflect the three genotypes (Table S5). Moreover, the TE copy numbers within the hybrid swarm have a distribution similar to a combination of allopatric trochilus and acredula, further supporting that the copies are inherited as intact blocks (haplotypes). However, a precise distinction between heterozygotes and homozygotes on MARB-a is still not possible15.
Statistical analysis
We used linear models with departure direction, winter longitude, migration distance and departure timing as response variables and the three genetic markers: MARB-a (a factor with two levels), InvP-Ch1 (a factor with three levels) and InvP-Ch5 (a factor with three levels) as explanatory variables. Models were constructed with R base package “stats”. We reported Type II ANOVA for models with more than one explanatory variable and no interactions and type III ANOVA results for models with interaction term by using R package “Car” (version 3.0-12)28. We initially constructed mixed effect models with timing of departure and tagging year as random factors however, this delivered singular fits due to insufficient sample sizes across categories. Normality of residuals was checked with a Shapiro–Wilk test. For carrying out circular statistics on autumn migration direction we used the R package “circular” (version 0.4-93). Watson’s U2 pairwise comparisons of different groups delivered the same results as linear models (Table S2 and Fig. S5). Circular means were identical to conventional linear means in our data set, which we take as another evidence that linear models are appropriate for the analysis of our data (Table S3 and Fig. S5). Maps in Figs. 1 and 2b and S1, S3 and S4 were created with R package “ggplot2” (version 3.3.6) using continent contours from Natural Earth, naturalearthdata.com/. Heat gradient over the maps in Fig. 1a–d were created with R package “gstat” (version 2.0-8) and the inverse distance weighting power of 3.0. Circular plots were created with ORIANA (version 4.02). All analyses were carried out with R version 4.1.1 (R Core Team 2021).
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
Source: Ecology - nature.com
