Restricted dispersal and inbreeding in a high‐elevation bird across the ‘sky islands’ of the European Alps

High‐elevation specialist species are threatened by climate change and habitat loss, and their distributions are becoming increasingly reduced and fragmented. In such a context, dispersal ability is crucial to maintain gene flow among patches of suitable habitat. However, information about dispersal is often lacking for these species, especially for those taxa that are usually considered as good dispersers such as birds. We adopted a landscape genomics approach to investigate dispersal in a climate‐sensitive high‐elevation specialist bird. Our aims were to assess the levels of gene flow within a wide mountain area, and to assess the effects of geographic distance and landscape characteristics on dispersal, by testing the isolation by distance (IBD) hypothesis against the isolation by resistance (IBR) hypothesis.

cleotide polymorphism (SNP) data by ddRAD sequencing.We then calculated site-and individual level genetic distances and individual inbreeding coefficients.To test IBD versus IBR, we related genetic distances to both geographic distances and different measures of landscape resistance by using maximum likelihood population effects models.
Results: Gene flow among breeding areas was partly restricted, and we found support for IBD, indicating that geographic distance limits snowfinch dispersal.Spatial patterns of genetic distances suggested that philopatry strongly contributed to determine the observed IBD.High inbreeding coefficients in several individuals indicated frequent mating among relatives.
Main Conclusions: Restricted dispersal and frequent inbreeding within 'sky island' systems can also occur in highly mobile species, because their potential ability to cover very large distances can be counteracted by high philopatry levels that are likely related to high dispersal costs.IBD and philopatry will increasingly hinder snowfinch dispersal among suitable areas within the future more restricted and fragmented breeding range, increasing the risks of local extinctions.

| INTRODUC TI ON
In mountain regions, the rapid variation of abiotic factors across short distances leads to a high variety of habitats within relatively small areas and this often results in highly patchy species distributions (Cadena et al., 2012;Körner & Ohsawa, 2006).For example, species that are strictly connected to high-elevation open areas often inhabit island-like patches of suitable habitat surrounded by unsuitable areas such as forested mountain slopes and valley bottoms.Moreover, they often show specific adaptations to the intrinsic harsh conditions of these habitats, making them vulnerable to environmental changes (e.g., Brambilla, Cortesi, et al., 2017;Imperio et al., 2013;Martin & Wiebe, 2004); highly specialized species with narrow niche breadth are indeed especially vulnerable to climate and habitat changes (Davey et al., 2012;Pearce-Higgins et al., 2015).
High-mountain ecosystems are particularly threatened by the interaction of (i) anthropogenic land use changes, such as increasing leisure activities and abandonment or intensification of livestock grazing (e.g., Brambilla et al., 2016;Gehrig-Fasel et al., 2007), and (ii) climate change, due to the effects of stronger temperature increase with increasing elevation (Liu & Chen, 2000;Pepin et al., 2015).
Given their rapid response to habitat and climate changes (Desrochers, 2010;Hallman et al., 2022;Lantz & Karubian, 2017), birds are an ideal biological model to investigate the consequences of such changes in high-mountain ecosystems.Accordingly, mountain birds are expected to change their distributions in order to track suitable climatic conditions and habitat (Reif & Flousek, 2012).On the one hand, this could lead to gains in suitable areas for mountain generalist, warm-adapted species (e.g., Ceresa et al., 2021), while on the other hand it could reduce the distributions of high-elevation specialist and cold-associated species because of uphill shifts toward a smaller ground area with increasing elevation in pyramidal mountain systems such as the European Alps (Elsen & Tingley, 2015).For some of these species, uphill shifts with reduction of the overall elevational range (e.g., Hallman et al., 2022), abandonment of breeding areas at lower elevations (Knaus et al., 2018) and population declines (e.g., De Gabriel Hernando et al., 2022;Lehikoinen et al., 2019;Rete Rurale Nazionale & LIPU, 2021) have been already detected in several regions.Furthermore, predictions of species distributions based on future climatic scenarios forecast strong breeding range reductions within a few decades for several high-elevation specialist birds (Brambilla, Caprio, et al., 2017;Brambilla, Rubolini, et al., 2022;De Gabriel Hernando et al., 2021).
In this context of increasingly patchy and reduced suitable habitat, the dispersal ability, (i.e., the capability to exchange individuals among different breeding areas) and the resulting gene flow are of crucial importance for population viability.Sufficient dispersal levels indeed reduce the risks of inbreeding and genetic drift (Frankham et al., 2010;Kvist et al., 2011) and allow recolonizations of vacant habitat patches after local extinctions (or colonization of new available areas), promoting maintenance of wider breeding ranges and larger populations.However, the available information about dispersal ability of high-elevation birds is extremely scarce, possibly because birds are perceived as extremely vagile and therefore they are supposed to be scarcely subject to possible limits to dispersal due to landscape structure (Kozakiewicz et al., 2018).Most studies about dispersal and inbreeding in high-elevation species are indeed focused on plants, or on animals that are usually considered as poor disperser such as ground-dwelling arthropods and reptiles (e.g., Atkins et al., 2020;Beckers et al., 2020;Morgan & Venn, 2017;Peyre et al., 2020;Slatyer et al., 2014;Tovar et al., 2020).Thus, at present it is not possible to properly evaluate the consequences of the ongoing and future habitat loss and fragmentation on the population viability of these threatened bird species.Such lack of knowledge is possibly due also to the logistical challenges of field work in mountain areas, which exacerbate the intrinsic difficulties of investigating dispersal in highly mobile organisms like birds (Cayuela et al., 2018;Paradis et al., 1998).Very large study areas are necessary to avoid describing only short-distance dispersal (Paradis et al., 1998), and this implies a huge sampling effort in studies based on mark-resight of single individuals.Studies based on satellite telemetry are limited as well, due to their high costs and the weight of the equipment that is not suitable for small birds.
Thus, information about dispersal can be obtained indirectly by means of population genetics, i.e. by comparing the degree of genetic differentiation among individuals or groups of individuals sampled in different areas (Frankham et al., 2010).This approach allows not only to assess levels of gene flow, but also to investigate which factors influence dispersal.For example, it is possible to detect if dispersal is affected by the geographic distance among reproductive areas (isolation by distance hypothesis, IBD), rather than by the characteristics of the landscape matrix separating them (isolation by resistance hypothesis, IBR).Landscape characteristics can indeed hinder dispersal or, conversely, facilitate it (e.g., Ceresa et al., 2015Ceresa et al., , 2023;;Klinga et al., 2019), even in extremely mobile species such as long-distance migratory birds (García et al., 2021).
Occurrence of IBD, IBR or unrestricted dispersal have different implications for conservation.IBR patterns emphasize the importance of dispersal corridors with specific characteristics to facilitate animal high-mountain conservation, isolation by distance, isolation by resistance, maximum likelihood population effects (MLPE) models, Montifringilla nivalis, white-winged snowfinch movements (Klinga et al., 2019), while IBD implies the need to consider geographic distances among habitat patches, e.g., to evaluate risks of isolation or the probability of (re)colonization of individuals after habitat restoration (e.g., Taylor et al., 2021;Wright et al., 2015).
In this study, we investigated dispersal of the white-winged snowfinch Montifringilla nivalis (hereafter snowfinch) within a wide area of the central-eastern European Alps (Figure 1).This spe-  Brambilla et al., 2018;Resano-Mayor et al., 2019).The species is less specialized during the non-breeding period, though still associated with high mountain areas (Bettega et al., 2020;Resano-Mayor et al., 2017); during this period, the snowfinch usually performs erratic movements, often in large flocks, looking for food resources and some birds may also migrate over short distances (Resano-Mayor et al., 2017;Resano-Mayor et al., 2020).While the breeding ecology and the habitat preferences of this species have been the object of several studies (Bettega et al., 2020;Brambilla et al., 2018;Brambilla, Cortesi, et al., 2017;Resano-Mayor et al., 2019), information about dispersal is still lacking.Given its sensitivity to climate change (Brambilla et al., 2018;Schano et al., 2021;Strinella et al., 2020), the evidence of an ongoing range contraction and local declines (Brambilla & Delgado, 2020;Knaus et al., 2018;Lardelli et al., 2022;Patrinat, 2019), as well as strong forecast declines for the next decades (Brambilla, Rubolini, et al., 2022;De Gabriel Hernando et al., 2021), the snowfinch is among the most threatened mountain birds in Europe and a highpriority species for conservation (Brambilla, Caprio, et al., 2017;Brambilla, Rubolini, et al., 2022;Schano et al., 2021).Yet, it is still Europe is likely due to the currently unknown population trends (BirdLife International, 2021;Brambilla, Bettega, et al., 2022).The European Alps host the largest European snowfinch population and their importance for the species survival will further increase in the future, because no other mountain range in the continent has a similar extension of high-elevation areas, where the suitable conditions for the species will persist in spite of climate change (i.e, climate refugia; Brambilla, Rubolini, et al., 2022).
By investigating snowfinch dispersal across a wide Alpine area, with characteristics representative of the entire mountain chain, our main goal was to assess the connectivity in a population that is crucial for the conservation of the species at a continental scale.Specifically, we aimed (i) to assess levels of gene flow within our study area, using high-resolution genomic data obtained through ddRAD (double digest restriction site-associated DNA) sequencing (Peterson et al., 2012) and (ii) to determine if snowfinch dispersal is influenced by the geographic distance among breeding areas or by the characteristics of the landscape matrix separating them, by testing the isolation by distance hypothesis (IBD) against the isolation by resistance hypothesis (IBR), using both habitat suitability and landscape connectivity models to obtain different estimates of landscape resistance.

| Habitat suitability and landscape connectivity models
To obtain landscape resistance surfaces, we first modelled habitat suitability for the snowfinch by adopting the maximum entropy approach (MaxEnt; Phillips et al., 2006), using the SDMtool package (Vignali et al., 2020) in R 4.2.1.1 (R Core Team, 2022).This widely used approach allows the use of presence-only data collected through different protocols (Elith et al., 2011).We used 951 snowfinch occurrence records (years 2010-2020, Figure S1) collected during the snowfinch breeding season (15 May-31 July) in multiple previous studies and surveys carried out in the Trentino-South Tyrol region by the research institutes involved in this study (Anderle et al., 2022;Brambilla et al., 2018;Brambilla, Cortesi, et al., 2017;Ceresa et al., 2021;Ceresa & Kranebitter, 2020;Chamberlain et al., 2016).Although the bird distribution data and bird DNA sampling were limited to the Trentino-South Tyrol region and a sector of the Stelvio massif in Lombardy, we modelled habitat suitability and landscape connectivity for a wider area (see Figure 1), as dispersal can also occur by crossing neighbouring regions.This area covers a large portion of the snowfinch breeding range in the Alps (25%-30%; Figure 1) and is representative of the whole Alpine range of the species, because it includes both inner-Alpine and more peripheral massifs.The procedure followed to model habitat suitability is described in the Supplementary materials (section S1).
We estimated landscape connectivity for the snowfinch within our study area by means of a potential connectivity approach, which integrates distribution models and connectivity models based on the output of the former (Rödder et al., 2016).We used a method based on the circuit theory, i.e., the Omniscape algorithm (Landau et al., 2021), to consider all potential connections among patches.This algorithm is especially recommended for habitat suitability maps (McRae et al., 2016) and represents a 'spatial generalization' of the more largely used Circuitscape method (Hall et al., 2021;McRae et al., 2013).We ran Omniscape in Julia 1.6 (Bezanson et al., 2017), using the inverse of habitat suitability from MaxEnt output as resistance raster and following the procedure and settings described in the Supplementary materials (section S2).

| Landscape resistance surfaces
To investigate the effects of landscape characteristic on snowfinch dispersal, we used two different landscape resistance surfaces, the first one representing the inverse of habitat suitability and the second one the inverse of landscape connectivity (calculated from MaxEnt and Omniscape output, respectively).In both cases, we calculated values of accumulated landscape resistance among all individuals by means of the Run_gdistance function of R package ResistanceGA 4.2.4 (Peterman, 2018).Relating multiple resistance surfaces to the observed genetic distances is a widely used approach, but it is usually adopted to compare the effects of different landscape elements (e.g., Amos et al., 2014;Ruiz-Gonzalez et al., 2014), or to consider landscape changes through time (e.g., Miller et al., 2018).Differently, we compared two conceptually different layers to assess if landscape resistance was better represented simply by habitat 'unsuitability', or by an approach which explicitly tries to describe potential movements through the considered landscape.

| Sampling and DNA extraction
We obtained blood samples from 85 snowfinches captured with mist-nets during the 2021 and 2022 breeding seasons (mid May-mid July) at 7 sampling areas (Figure 1).Our sample included 69 adult birds (2nd year individuals or older) and 16 recently fledged juveniles.The sampling areas were located between 24 and 106 km from each other (mean distance ± SD: 61 ± 35 km).This is a sufficiently large distance range to study dispersal in passerine birds, because their dispersal distances very rarely exceed 100 km according to the information available for many species (e.g., Paradis et al., 1998;Rushing et al., 2021;Winkler et al., 2005).Within each sampling area, we captured birds at several different sites up to a few kilometres from each other.We sampled between 9 and 30 birds per area, excepted a logistically challenging sampling area (Vizze) where we only captured 2 birds.Almost all captured adults showed well developed incubation patches (females) or cloacal protuberances (males), and juveniles were very recently fledged individuals according to plumage characteristics and behaviour.This confirmed that we sampled reproductive individuals at their breeding sites or their offspring.We obtained blood samples (20 μL) by puncturing the brachial vein (Owen, 2011) and we stored them in ethanol.We extracted DNA using E.Z.N.A.® Tissue DNA Kit (Omega Bio-Tek) following the manufacturer's protocol.

| ddRAD sequencing and initial data filtering
Library preparation and sequencing, as well as initial raw data analysis and SNP calling, were carried out at IGA Technology Services (Udine, Italy) as described in the Supplementary materials (section S3).
We then used the program PLINK 1.90 beta (Chang et al., 2015) to filter SNPs for linkage disequilibrium (−indep function) and missing data, setting minimum allele frequency at 0.01 and excluding individuals with >20% missing genotypes (with functions -maf and -mind, respectively).Through this procedure, we obtained a dataset including 27,072 SNPs and all 85 snowfinches, given that no individual showed excessive missing data; genotyping rate was 0.968.
We then converted the PLINK dataset into the different formats required in downstream analyses using PGDSpider 2.1.1.5(Lischer & Excoffier, 2012).Waples & Anderson, 2017).Therefore, we performed our analyses by using both the full dataset and a reduced dataset (N = 74, details in Table S3) created by including only one bird within each group of first-degree related individuals and removing only juveniles, because the co-occurrence of closely related adults in the same area is very likely a consequence of the philopatry of these individuals rather than of non-random sampling.Further details about this purging procedure are provided in the Supplementary materials (section S4).For all the downstream analyses, the results were similar among the two datasets, therefore we report the results obtained for the full sample, while those of the reduced dataset are reported in the Supplementary materials (section S12).

| Genetic distances and population structure
As a first data exploration, we calculated basic statistics such as expected and observed heterozygosity (H e and H o ) in the program Arlequin 3.5 (Excoffier & Lischer, 2010) and individual methodof-moments inbreeding coefficients (F) using the PLINK function -het.Given that F estimates can be influenced by demographic history (Polašek et al., 2010), and may partly reflect also ancient, shared ancestry besides recent inbreeding, we further deepened our approach by examining runs of homozygosity (ROH) within our sample.Short ROH segments indicate loss of genetic diversity due to historical founder effect or bottlenecks, while the occurrence of long ROH indicate recent inbreeding (e.g., Martin et al., 2023;Pilot et al., 2014).Following Martin et al. (2023) We then assessed genetic differentiation at both site-and individual levels, by calculating pairwise F ST values among sampling areas using Arlequin and pairwise identity-by-state (IBS) distances among individuals in PLINK (function -distance 1-ibs).We also investigated the overall genetic population structure because this can help to identify especially important barriers to dispersal (e.g., Ceresa et al., 2018;Millions & Swanson, 2007).Using and comparing different statistical methods is especially recommended when population structure is weak (Frosch et al., 2014;Kraus et al., 2016), which is a possible scenario within our study system.Therefore, we adopted two different approaches: a discriminant analysis of principal components (DAPC) implemented in R package adegenet 2.1.8(Jombart, 2008;Jombart & Ahmed, 2011), and the Bayesian approach implemented in the program Structure 2.3.4 (Falush et al., 2003;Pritchard et al., 2000).In Structure, we fitted a model with population admixture and correlated allele frequencies to estimate the most likely number of distinct genetic clusters (K), including spatial information about sampling sites (LOCPRIOR procedure).We carried out 5 independent runs for each value of K between 1 and 10, with a burn-in period of 10,000 iterations and 10,000 Markov chain Monte Carlo replications.We adopted this range of K to explore a wide range of possible scenarios, given the lack of previous information about snowfinch population structure within the study area.To better identify the actual number of genetic clusters, we used the Structure results to calculate the ad hoc ΔK statistic (Evanno et al., 2005).DAPC was carried out by means of the find.clusters function of R package adegenet 2.1.8,which identifies genetic clusters through the K-means algorithm (Jombart et al., 2010).We ran the analysis by retaining all the principal components to avoid the loss of information, while the only advantage of using only a part of them would have been the reduction of computational time (Jombart & Collins, 2015).We set K = 10 as the maximum number of clusters and kept the default number of iterations to be used in each run of the algorithm (N = 100,000).

| Testing isolation by distance and by resistance
As a first step to test our IBD and IBR hypotheses, we explored the relationship between pairwise geographic distances/landscape resistances and individual genetic distances by using Mantel and partial Mantel tests (Amos et al., 2014;Cushman et al., 2013).

Given the possible interpretation problems of Mantel correlation r
in landscape genetics (Meirmans, 2015), we further deepened our comparison of IBD and IBR hypotheses by means of maximum likelihood population effects (MLPE) models (Clarke et al., 2002), which outperform other regression-based methods for model selection in landscape genetics (Shirk et al., 2018).We only carried out individual-level analyses, given the low number of sampling areas.
We performed Mantel and partial Mantel tests using the mantel function in R package ecodist 2.0.9 (Goslee & Urban, 2007) with 10,000 permutations.We used both 1/suitability and 1/connectivity matrices to test IBR (hereafter, IBR SUIT and IBR CONN ); when IBR SUIT and IBR CONN were significantly supported by Mantel tests (p < 0.05), we considered the model obtaining the highest Mantel correlation r as the best description of landscape resistance (Amos et al., 2014).When both IBD and IBR were significantly supported, we used partial Mantel tests in a causal modelling framework to assess if variation of genetic distances was better described by IBR or by IBD (Amos et al., 2014;Cushman et al., 2013), i.e., we calculated the relationship between landscape resistance and genetic distances after accounting for the effect of geographic distances, and vice versa.
We fitted MLPE models by using the MLPE.lmmfunction in R package ResistanceGA 4.2.4,keeping the default specification 'scale = TRUE' in order to allow comparisons among the effects of different explanatory variables.We fitted a set of competing models using alternatively Euclidean geographic distances (IBD model), 1/suitability cost distances (IBR SUIT model) and 1/connectivity cost distances (IBR CONN model) as explanatory variables, and genetic distances as response variable.This allowed us to compare our hypotheses while avoiding including highly correlated explanatory variables (landscape resistances and geographic distance) within the same model.We also fitted a null model with no distance or landscape resistance effect.We ranked the models based on their Akaike Information Criterion, and we considered as substantially supported those models with ΔAIC < 2 (Burnham & Anderson, 2002).We then calculated bootstrapped 95% confidence intervals of explanatory variables' effects with the tidy function of the R package broom.mixed0.2.9.4 (Bolker & Robinson, 2022), and we considered a variable effect to be significant when 95% confidence intervals did not include 0. For each model, we also calculated R 2 by using function r.squaredGLMM in R package MuMIn 1.47.1 (Bartoń, 2019).

| Habitat suitability and landscape connectivity models
Distribution modelling led to an accurate and robust description of breeding habitat suitability for the snowfinch in the study area: accuracy statistics showed no decline in model accuracy and discriminatory ability, when tested on the independent dataset (TSS train = 0.61; TSS test = 0.69; AUC train = 0.87; AUC test = 0.91).
Landscape connectivity modelling and the resulting cumulative current map indicated high connectivity only along the main mountain ridges, with wide areas of low connectivity corresponding with valley floors and forested mountain slopes (Figure 2).The resistance surfaces obtained from both habitat suitability and landscape connectivity depict a higher landscape resistance of valley floors and woodlands, compared to the very low resistance of high-mountain areas (Figure 3).While the two surfaces may seem quite different according to their visual representation, the respective matrices of cost distances among sampled birds were highly correlated (r = 0.98, p < 0.001).Cost distances among sampling sites are summarized in Table S2.

| Genetic distances and population structure
IBS individual genetic distances ranged between 0.079 and 0.165 (mean ± SD: 0.156 ± 0.007; Figure 4).Pairwise population differentiation among sampling areas was low, but always statistically significant, except for a few cases of non-significant F ST values relative to the Vizze sampling area (Table 1), where we only sampled 2 birds.Cima d'Asta, a peripheral breeding area (see Figures 1 and 2), showed a higher differentiation with the other areas than all other pairs of sites (Table 1).
Population structure analyses did not provide evidence of population clustering in our study area.Using the program Structure, we obtained the highest likelihood for K = 5, but ΔK (which cannot be calculated for K = 1) did not show a clear single peak (Table S6).Using the reduced dataset, the results changed partly (Table S9), consistently with the reported influence of family groups on this kind of analysis (Anderson & Dunham, 2008), but also in this case we did not obtain clear ΔK peaks (Table S9).According to DAPC analysis, K = 1 obtained the lowest value of Bayesian Information Criterion (BIC) although the difference with K = 2 was very reduced (ΔBIC = 1.43).
However, using the reduced dataset the difference with K = 2 was clearer (ΔBIC = 2.40).In addition, the true number of K is often indicated by an 'elbow' in the curve depicted by BIC values as a function of K (Jombart & Collins, 2015), and we did not find this pattern for K = 2 or for higher numbers of clusters in both the full and the reduced sample (Figures S3 and S5).
For several individuals, high values of the inbreeding coefficient F indicated mating between closely related birds (Figure 5).For 17 out of 85 individuals, F exceeded 0.0625, i.e., indicated a relationship of first-cousins or closer (e.g., Polašek et al., 2010); in some more individuals, F approached this threshold (Figure 5).These individuals occurred in all sampling areas (excepted Vizze), but were more   common in the peripheral breeding area Cima (F > 0.0625 in 8 out of 14 snowfinches), which was reflected in a clearly higher mean F (mean ± SD: 0.076 ± 0.047) than across the entire sample (0.039 ± 0.036).For 3 individuals, 2 of which sampled at Cima d'Asta, we obtained F > 0.125 (half siblings or closer).Mean F for each sampling area, as well as observed and expected heterozygosities are provided in Table S4 in the Supplementary materials.ROH segments confirmed the occurrence of recent inbreeding in our study area: we found long ROH (> 1 Mb) for 62 out of 85 individuals, and for 40 birds the longest ROH exceeded 4 Mb (Table S5).All but one individual with F > 0.0625 showed ROH >4 Mb, and in some birds, we found very long ROH segments, ranging between 10 and 17 Mb (Table S5).
Both maximum and total lengths of long ROH were highly correlated with F coefficients (r = 0.74, p < 0.001 and r = 0.84, p < 0.001, respectively).Consistently with the higher F values at Cima d'Asta, birds from this site showed a higher number of long ROH and proportionally more individuals with ROH >4 Mb (11 out of 14) than in the other sampling areas (Table S5).

| Testing isolation by distance and by resistance
Mantel tests significantly supported IBD as well as IBR SUIT and IBR CONN, with a higher Mantel correlation r for IBD (Table 2).
According to the partial Mantel test, IBD was still significantly supported after accounting for IBR (both 1/suitability and 1/connectivity), while the inverse did not occur (Table 2), which clearly indicates a stronger support for IBD than for the two IBR models.
According to MLPE modelling, IBD was the most supported hypothesis, with a large difference in AIC values with both IBR SUIT (ΔAIC = 55.4) and IBR CONN (ΔAIC = 79.3).The difference in AIC between IBD and the null model was even larger (ΔAIC = 295.9).
Marginal R 2 values were higher for IBD than for IBR SUIT and IBR CONN (Table 3).Geographic distances and both 1/suitability and 1/connectivity cost distances all showed a significant positive effect on genetic distances (95% CI did not include 0), but the effect was stronger for IBD (95% CI: 0.00157-0.00210)than for IBR SUIT (95% CI: 0.00143-0.00196)and IBR CONN (95% CI: 0.00134-0.00184).The snowfinch is highly mobile during the non-breeding period (Bettega et al., 2020), performing erratic movements and occasionally short-distance migrations (Resano-Mayor et al., 2017, 2020).However, high flight efficiency and the capability to cover large distances during migration do not necessarily imply high dispersal distances or rates: dispersal can indeed be limited in even highly mobile birds, such as long-distance migrants, in case of high philopatry (e.g., Ceresa et al., 2016;García et al., 2021;Hansson et al., 2002;Rönkä et al., 2021).The relationship between migration habits and dispersal ability in birds is still not fully understood, and it is possible that they are decoupled (Chu & Claramunt, 2023).
Limited dispersal due to philopatry likely occurs also in our study system, given the detection of closely related individuals within the same breeding areas (Figure 4).This pattern is evident also when considering the reduced dataset, where first-degree relationships involving juveniles were removed (Figure S6), supporting high philopatry in our study areas.This suggests that snowfinches might have to face high dispersal costs and risks, such as, e.g., lower survival probability or breeding success.During breeding, the snowfinch has very specific requirements for both nesting sites and foraging areas (Brambilla et al., 2018), therefore a good knowledge of a habitat patch could strongly favour locally born individuals over immigrants.Philopatry apparently strongly contributes to shape the observed IBD pattern, determining a strong increase in genetic distances within relatively short geographic distance, followed by a plateau for distances > 40 km (Figure 4).
Although this pattern would be better investigated through a less clustered sampling (and, therefore, more continuous distance classes), it is consistent with a scenario where many individuals are philopatric, but dispersing ones can cross large distances; this is a common pattern in passerine birds, described by the fat-tailed distribution of their dispersal distances (e.g., Ceresa et al., 2016;Paradis et al., 2002;Van Houtan et al., 2007)  the MLPE model for IBD (anyway higher than for IBR): as MLPE are linear models, they may have partly underestimated the variance explained by geographic distances.
The high inbreeding coefficients found in many snowfinches reveal frequent mating between related individuals (third-degree, sometimes second-degree relationships) in the considered breeding areas, as confirmed also by ROH lengths.In combination with our results about population differentiation and structure, this suggests that the exchange of individuals among breeding areas is sufficient to generate gene flow (very few ones per generation should be sufficient, see Slatkin, 1987), but not large enough to provide sufficient opportunities to mate with unrelated individuals.Inbreeding may lead to inbreeding depression, with deleterious effects such as reduced survival and breeding success (e.g., Hemmings et al., 2012;Niskanen et al., 2020).We found that mean F was the highest at the most peripheric breeding area (Cima d'Asta), but offspring of related individuals were found across the entire study area, suggesting possible fitness reductions due to inbreeding depression across large parts of the Alpine range.
The combination of IBD, philopatry and inbreeding, jointly with the forecasted decrease of breeding range (Brambilla, Caprio, et al., 2017;Brambilla, Rubolini, et al., 2022), depicts an unfavourable scenario for the snowfinch.Within a more reduced and more patchy distribution in the future, a weak tendency to disperse and the effect of distance to dispersal will increasingly hinder the exchange of individuals among breeding areas, increasing the risks of local extinctions especially in the more isolated patches (lack of immigrants, reduced fitness due to inbreeding, interaction with stochastic events), and reducing the probability of recolonizations after the local extinctions.Therefore, the forecasted breeding distribution based on the occurrence of suitable environments under future climatic scenarios for the next decades is probably too optimistic.
The species could rapidly disappear from many habitat patches through the aforementioned processes, even when they remain suitable for breeding.It is possible that the patterns we observed, in synergy with climate change and habitat loss, were already involved in the recent reduction of the Alpine breeding range of the species, occurring especially in the lower and more peripheric massifs (e.g., Knaus et al., 2018;Lardelli et al., 2022).
Our results highlight the need to maintain the continuity of snowfinch breeding range, in order to facilitate dispersal.One way to buffer against climate change could be a targeted sward management through grazing, which may help to compensate for earlier snowmelt by improving food availability during nestlings rearing (Brambilla et al., 2018).More generally, maintaining extensive grazing in alpine pastures is crucial to avoid shrub encroachment, which erodes the extension of snowfinch foraging areas near their lower elevational limit; at the same time, overgrazing should be avoided, because it causes grassland and soil degradation (Garcia-Pausas et al., 2017) and can be detrimental for grassland birds (Brambilla, Gustin, et al., 2020;Pavel, 2004).In addition, avoiding the creation of new infrastructures in high-mountain areas is very likely beneficial to this species.Although snowfinches also breed on buildings and winter tourism provides additional food resources during the non-breeding period, ski-pistes and related infrastructures severely impact high-alpine grasslands, affecting arthropod abundances during nestlings rearing (Rolando et al., 2007) and decreasing foraging habitat suitability for the species (Brambilla et al., 2018).Therefore, grassland restoration at disused touristic areas/infrastructures and ski-pistes revegetation (Caprio et al., 2011) are other important measures to reduce habitat loss and degradation.Based on our results, these actions are recommended not only in the currently suitable large breeding areas, but also in relatively small patches hosting few breeding pairs.These patches may be used as stepping-stones to other occupied patches/massifs, contributing to the overall breeding range continuity and promoting long-distance dispersal (e.g., Saura et al., 2014).This is especially important in those Alpine areas where the breeding range is more scattered, both currently and in the predicted future scenarios, e.g. in the Dolomites and more southern massifs and in the outermost north-east of the species distribution in Austria (see Figure 1 and S7).Promoting connectivity and avoiding local extinctions in these peripheral breeding areas is crucial to counteract the overall Alpine breeding range reduction, and would help to maintain a larger population size and higher genetic diversity.To this aim, the detailed habitat suitability model provided in this study, and the coarser scale Alpine-wide models provided by Brambilla et al. (2022a), represent tools that can be used to evaluate the risk of isolation of breeding areas and to identify potential stepping-stones patches.
Our study is focused on within-mountain range dispersal, and IBD/IBR patterns could be different at larger scales, e.g., continentalscale dispersal among different mountain chains.Future studies at such a scale would be especially useful to assess if the more scarce and scattered population of lower mountain ranges receive immigrants from the main breeding areas such as the Alps.The gene flow limitations and IBD found already at the regional scale in our study suggest that the very large distances separating different mountain ranges may strongly limit dispersal.
Our results show that, within the scattered and discontinuous ranges of high-elevation species, restricted dispersal and frequent inbreeding can occur also in highly mobile organisms.
While the snowfinch is potentially able to cover very large distances (like during the post-breeding erratic movements or shortdistance migrations), such ability is apparently counteracted by a highly philopatric behaviour.Also other high-elevation birds with patchy distributions and very specific habitat requirements could show similar dispersal patterns as observed in the snowfinch in this study, such as gene flow limitations through IBD or IBR, philopatry and mating among relatives.Signs of high philopatry have indeed been reported also for the rock ptarmigan (Bech et al., 2009) and for the water pipit Anthus spinoletta (Ceresa et al., 2023), suggesting high dispersal costs also in these species.The water pipit also showed a clear IBR pattern, with a stronger effect of landscape resistance where alpine grasslands are less extended (Ceresa et al., 2023).Jointly with the results of the present study, cies is tightly connected to high mountain open areas above the treeline and its range includes some European mountain chains (Cantabrian mountains and Pyrenees, Alps, Apennines, mountains of Corsica and of the Balkan peninsula; Brambilla, Resano-Mayor, et al., 2020) and the main mountain systems of western and central Asia, up to Mongolia and western China (BirdLife International, 2022).Habitat needs during breeding are very specific, requiring both adequate nest sites like rock crevices or artificial structures and arthropod-rich foraging areas, represented by short grasslands and snow patches (Alessandrini et al., 2022; classified as 'least concern' in most European countries and in the global and European IUCN red lists (BirdLife International, 2021).Such underestimation of the threatened status of the species in F I G U R E 1 Study area considered to investigate snowfinch dispersal within the European Alps: (a) DNA sampling sites (red dots) and region where we collected bird occurrence data to model habitat suitability and landscape connectivity (Trentino-South Tyrol, north-eastern Italy, blue line); (b) area covered by the European Alps (light grey), fine-scale species distribution calculated in this study (green areas), and large-scale species distribution across the rest of its Alpine breeding range (blue areas), as modelled by Brambilla et al. (2022a).
Based on field observations, we suspected that closely related individuals occurred within our sample (full siblings, parentsoffspring), especially because recently fledged juveniles are still fed by adults and follow them in the foraging areas.Whether purging or not closely related individuals from population genetics data set (and the procedure to follow) is still a relatively controversial issue, and both maintaining large family groups or indiscriminately removing close relationships imply the risk to obtain biased estimates(Anderson & Dunham, 2008;Rodríguez-Ramilo & Wang, 2012; , who studied a passerine bird with similar generation length than the snowfinch, we considered ROH >1 Mb as long ROH, indicating last common ancestor <50 generations, and ROH >4 Mb were considered as an indication of recent inbreeding.We looked for long ROH in our sample by means of the PLINK function -homozyg with the default program settings, considering the original SNPs dataset (not filtered for MAF and LD; Meyermans et al., 2020) which included 64,612 SNPs.

F
I G U R E 2 Maps of (a) habitat suitability (according to MaxEnt modelling) and (b) landscape connectivity (Omniscape modelling) for the snowfinch in the study area.Blue dots represent DNA sampling sites.
Landscape resistance surfaces obtained from (a) habitat suitability and (b) landscape connectivity models.Blue dots represent DNA sampling sites.

F
I G U R E 4 Geographic pattern of individual genetic distances among snowfinches: (a) scatterplot showing the occurrence of more closely related individuals captured within the same study areas (see the IBS values between c. 0.08 and 0.14 for distances approaching 0 km) and the smoothing spline (black line) interpolating the data; (b) detail of the smoothing spline, with Bayesian confidence intervals in grey.

TA B L E 1
Sample sizes (N) and genetic differentiation (F ST values) among sampling areas.Our results about genetic differentiation and structuring indicate the occurrence of gene flow among all the breeding areas considered; gene flow is anyway partly restricted, as highlighted by the significant pairwise F ST values.Both Mantel test and MLPE models indicate that, at least at the scale we considered, geographic distance is involved in restricting snowfinch dispersal, irrespective of the characteristics of the landscape matrix among breeding areas.
. Many of the low IBS values observed within breeding areas (N = 23 for IBS < 0.14) are indeed distances among adult birds, i.e., individuals that already had the opportunity to disperse.In addition, IBD acting already at relatively short distances is fully consistent with significant F ST values also among those sampling sites that are close to each other (i.e., 20-40 km).The clearly nonlinear IBD pattern depicted in Figure 4 may also partly explain the relatively low marginal R 2 of F I G U R E 5 Method-of-moments individual inbreeding coefficients F across sampling areas.The dashed line represents the F value indicating a relationship of first-cousins (F = 0.0625).TA B L E 2 Results of Mantel and partial Mantel tests, used to test the isolation by distance (IBD) and the isolation by resistance (IBR) hypotheses based on individual genetic distances.
IBR SUIT and IBR CONN indicate, respectively, the use of the inverse of habitat suitability and of the inverse of landscape connectivity to represent landscape resistance.Maximum likelihood population effects models, ranked according to the Akaike Information Criterion (AIC) and R 2 (marginal and conditional).IBR SUIT and IBR CONN indicate, respectively, the use of the inverse of habitat suitability and of the inverse of landscape connectivity to represent landscape resistance.
TA B L E 3