Evolutionarily distinct lineages of a migratory bird of prey show divergent responses to climate change

Abstract

Accurately predicting species' responses to anthropogenic climate change is hampered by limited knowledge of their spatiotemporal ecological and evolutionary dynamics. We combine landscape genomics, demographic reconstructions, and species distribution models to assess the eco-evolutionary responses to past climate fluctuations and to future climate of an Afro-Palaearctic migratory raptor, the lesser kestrel (Falco naumanni). We uncover two evolutionarily and ecologically distinct lineages (European and Asian), whose demographic history, evolutionary divergence, and historical distribution range were profoundly shaped by past climatic fluctuations. Using future climate projections, we find that the Asian lineage is at higher risk of range contraction, increased migration distance, climate maladaptation, and consequently greater extinction risk than the European lineage. Our results emphasise the importance of providing historical context as a baseline for understanding species' responses to contemporary climate change, and illustrate how incorporating intraspecific genetic variation improves the ecological realism of climate change vulnerability assessments.

Fig. 1. Map showing sampling localities for lesser kestrel genetic data were obtained and the geographic distribution of occurrence records.

a Localities (coloured circles; n = 16) from which we obtained double-digest Restriction-Site Associated DNA (ddRAD) or mitogenome data; acronyms and sample size for each locality are reported in Supplementary Table 1. b Map of breeding (dark-shaded circles) and non-breeding (light-shaded circles) occurrence records (original occurrence data pooled for 2.5 arc-minute grid cells) used for species distribution modelling of Western (orange circles) and Eastern (blue circles) evolutionarily significant units (ESUs). Lesser kestrel breeding and non-breeding distribution ranges are shown in brown and blue, respectively. Background maps were obtained from the rnaturalearth v.0.3.2 R package. Lines delimiting countries are shown to facilitate map interpretation and do not necessarily represent accepted national boundaries. The lesser kestrel illustration (male) is used with permission from Martí Franch ©. Data underlying all components of Fig. 1 are provided at 10.5281/zenodo.14988067.

Fig. 2. Genetic differentiation and gene flow among lesser kestrel populations.

a Admixture component profiles for each sampling locality at K = 2 (best K value based on cross-validation; Supplementary Fig. 11) based on double-digest Restriction-Site Associated DNA (ddRAD) data, showing the Western (orange) and Eastern evolutionarily significant units (ESUs) (blue). Supplementary Table 1 shows acronyms and details of sample size for each sampling locality. b Principal Component Analysis (PCA) based on 27,853 unlinked single-nucleotide polymorphisms (SNPs) highlighting the two ESUs (orange: Western ESU; blue: Eastern ESU) and the four fine-scale genetic clusters identified by fineRADStructure (Supplementary Fig. 6; coloured ellipses: Iberian peninsula in red, Italian and Balkan peninsulas - including Turkey - in yellow, Israel in green, Asia in blue). c EEMS-predicted barriers to gene flow (orange) showing a main barrier between the Western and Eastern ESUs. Pie charts show the frequencies of the three major mitochondrial haplogroups for each sampling locality (n = 89 individuals). For each of the fine-scale genetic clusters (coloured boxes; colours defined in b), the arrows show the predicted fraction of immigrating individuals per generation (proportional to the size of the arrow; numbers shown close to  each arrow) from other clusters (prediction based on BayesAss3-SNPs analysis). d Maximum-likelihood population tree inferred in Treemix using Falco tinnunculus (TIN) as an outgroup. e Residuals of the observed versus predicted squared allele frequency difference inferred in Treemix, expressed as the standard error of the deviation. Residuals above zero represent populations that share more genetic variation than predicted by the best-fit tree, potentially due to gene flow or shared ancestral genetic variation. Negative residuals represent populations that share less genetic variation than predicted by the best-fit tree. Data underlying all components of Fig. 2 are provided at 10.5281/zenodo.14988067. Background maps were obtained from the rnaturalearth v.0.3.2 R package. Lines delimiting countries are shown to facilitate map interpretation and do not necessarily represent accepted national boundaries.

Fig. 3. Ecological differentiation between Western and Eastern lineages of the lesser kestrel.

a Reconstructed breeding and non-breeding ranges for Western and Eastern evolutionarily significant units (ESUs) obtained through species distribution models based on selected bioclimatic variables, largely matching the known distribution ranges of the species (Fig. 1, Supplementary Fig. 14). b (Top) Habitat use for each combination of lineage and season calculated as the proportion of seven main land cover categories (pooled from 22 original land use categories from the Copernicus Land Monitoring Service; Supplementary Table 3) in all 2.5 arc-minute grid cells with breeding or non-breeding occurrence records. (Bottom) Heatmap of habitat preference/avoidance for each of the seven main land cover categories based on a sign test of habitat selection. White cells (f ~ 0.5) indicate habitats used proportionally to their availability, red cells (f > 0.5) indicate preferred habitats and blue cells (f < 0.5) indicate avoided habitats. c–d Climatic niche comparison between Western and Eastern ESUs in a two-dimensional space defined by the first two axes of a principal component analysis (PCA) of available climatic conditions across breeding (c) and non-breeding (d) ranges of the Western (orange contour lines) and Eastern (blue contour lines) lineages. The x-axes represent a gradient of increasing coldness (warmer climates in darker red, colder climates in darker blue), whereas the y-axes represent a gradient of increasing precipitation (beige to darker blue). The solid and dashed contour lines represent 100% and 75% of the available (background) climate, respectively. Coloured areas represent climatic niches (kernel densities of the climatic conditions at occurrence records) of Western (orange) and Eastern (blue) ESUs, with darker colours denoting higher densities and transparency adjusted to facilitate the visualisation of overlaps. Schoener’s D index of niche overlap (0 = no overlap, 1 = full overlap) and the p-values of niche equivalency tests performed with the ecospat.niche.equivalency.test function from the ecospat R package (option overlap.alternative = “lower” and 1000 random permutations) are reported. Principal components 1 (PC1) and 2 (PC2) were flipped in panel d to facilitate the comparison of climatic axes shown in panel c. Data underlying all components of Fig. 3 are provided at 10.5281/zenodo.14988067. Background maps were obtained from the rnaturalearth v.0.3.2 R package.

Fig. 4. Geography and climate explain spatial genomic variation of lesser kestrels.

a Pairwise genetic distance (Φ_ST/1 – Φ_ST) correlates positively with geographic and climatic distance. Mantel tests and their associated p-values (one-sided) are reported. Background colours reflect the density of points (blue indicating low density and red indicating high density) and show a discontinuity consistent with a scenario of two distant and differentiated genetic clusters. Distances between localities in the same or different evolutionarily significant units (ESUs) are shown as triangles and circles, respectively. b Principal component analysis (PCA) of 61 climate-associated single-nucleotide polymorphisms (SNPs) showing that the two putative adaptive units (AUs) coincide with the identified evolutionarily significant units (ESUs). c Allelic turnover functions relative to the two highest ranking bioclimatic variables from a gradient forest (GF) analysis of 61 climate-associated SNPs, i.e., precipitation of the coldest quarter (BIO19) and temperature annual range (BIO7). Y-axis values report the cumulative importance of SNPs in the GF models, which reflects the total amount of allele frequency turnover across the environmental gradient. Thin lines show allelic turnover functions for each of the 61 candidate SNPs. Thick blue and red lines show allelic turnover functions across all candidate SNPs and thick black lines across all putatively non-adaptive reference SNPs. Higher turnover values for candidate SNPs compared to neutral SNPs evidence the stronger association of candidate SNPs with climate. Circles at the top represent sampling localities coloured based on the ESU they belong to (orange: Western; blue: Eastern) ordered along the BIO19 and BIO7 gradients. d (Left) Hierarchical clustering of associations between bioclimatic variables (columns) and allele frequencies for the 61 climate-associated SNPs (rows) (Spearman correlation, absolute values). Bioclimatic variables associated with temperature and precipitation are coloured in red and blue, respectively. (Right) Tiles are coloured in grey when the candidate SNP is found within a gene identified, through a literature search (Methods), as being related to local adaptation (first column), phenotypic traits important for local adaptation (second column), adaptation to urban environments or domestication (third column), and/or stress response (fourth column) in vertebrates. Data underlying all components of Fig. 4 are provided at 10.5281/zenodo.14988067.

Fig. 5. Climate fluctuations, demographic history, and hindcasted breeding and non-breeding distribution ranges of lesser kestrels.

a Temperature anomaly as inferred from the EPICA (European Project for Ice Coring in Antarctica) Dome C ice core; b DIYABC best-supported scenario showing divergence and admixture times (mean and 95% confidence interval [CI] from 1000 out-of-bag testing samples using a set of broad priors drawn from uniform distributions [Supplementary Table 9]) of four fine-scale genetic clusters (Iberian peninsula, Italian and Balkan peninsulas, Israel and Asia; Fig. 2b); c effective population size (N_e) changes through time for the Western evolutionarily significant unit (ESU) estimated by MSMC2 (steps with bootstraps shown as faded steps); d Bayesian Skyline Plots (BSP) showing N_e changes through time (ribbons indicating 95% highest posterior density [HPD] intervals with lines indicating the median) obtained from mitogenome sequences; e extent of predicted breeding and f non-breeding ranges. The grey shaded area across panels represents the Last Glacial Maximum (LGM). In panels d–f estimates for the Western and Eastern ESUs are shown in orange and blue, respectively. In panels e–f estimates are shown for the last 20,000 years (2,000-years steps) and dotted lines connect estimates for 130 kya (when the two ESUs had not diverged yet) to estimates for 20 kya; g Predicted breeding and non-breeding ranges for Western and Eastern ESUs at selected timepoints; the timepoints shown in panel g are highlighted with grey vertical lines in panels (a–f). Data underlying all components of Fig. 5 are provided at 10.5281/zenodo.14988067. Background maps were obtained from the rnaturalearth v.0.3.2 R package.

Fig. 6. Forecasting lesser kestrel distributional responses to climate change.

a Predicted breeding (dark-shaded colours) and non-breeding (light-shaded colours) range in the present and in the future (2041-2070 and 2071-2100), for Western (orange) and Eastern (blue) evolutionarily significant units (ESUs), using an ‘extreme warming’ future climate (UKESM1-0-LL; SSP5-8.5). Green circles show the centroids of breeding distributions for each of the ESUs and the distance between centroids is shown above the line connecting them. Geographic isolation between the two ESUs is expected to increase in the future. b Trends from 6 kya to the future for: (top) the extent of breeding range, (middle) the extent of non-breeding range, and (bottom) migratory distance, showing divergent patterns between Western and Eastern ESUs. Migration distance is calculated as the minimum distance between the breeding distribution and the non-breeding distribution range centroids. For future time periods, estimates obtained from the ‘extreme warming’ future climate are shown as dashed lines, while those from the ‘moderate warming’ future climate (GFDL-ESM4; SSP3-7.0) are shown as solid lines. Data underlying all components of Fig. 6 are provided at 10.5281/zenodo.14988067. Background maps were obtained from the rnaturalearth v.0.3.2 R package.

Fig. 7. Genetic offsets for lesser kestrel evolutionary lineages.

a, b Density plots of genetic offsets in 2041–2070 (a) and 2071–2100 (b) for all 2.5 arc-minute grid cells with breeding occurrence records within the Western (orange) and Eastern (blue) evolutionarily significant units (ESUs), showing higher offsets for the Eastern ESU. Median values for each ESU are shown as coloured vertical lines. c–d Genetic offsets across the current breeding range in 2041–2070 (c) and 2071–2100 (d) based on projections using an ‘extreme warming’ future climate (UKESM1-0-LL; SSP5-8.5; estimates for a ‘moderate warming’ future climate are shown in Supplementary Fig. 20). The highest offsets are in the western and northern sectors of the Eastern ESU’s current range. Data underlying all components of Fig. 7 are provided at 10.5281/zenodo.14988067. Background maps were obtained from the rnaturalearth v.0.3.2 R package.