Forecasting climate change response in an alpine specialist songbird reveals the importance of considering novel climate

Species persistence in the face of climate change depends on both ecological and evolutionary factors. Here, we integrate ecological and whole‐genome sequencing data to describe how populations of an alpine specialist, the Brown‐capped Rosy‐Finch (Leucosticte australis) may be impacted by climate change.

genomic variants associated with environmental data. We modelled future climate change impacts on habitat suitability using ecological niche models (ENMs) and impacts on putative local adaptation using gradient forest models (a genetic-environment association analysis; GEA). We used the metric of niche margin index (NMI) to determine regions of forecasting uncertainty due to climate shifts to novel conditions. Results: Population genetic structure was characterized by weak genetic differentiation, indicating potential ongoing gene flow among populations. Precipitation as snow had high importance for both habitat suitability and changes in genetic variation across the landscape. Comparing ENM and gradient forest models with future climate predicted suitable habitat contracting at high elevations and population allele frequencies across the breeding range needing to shift to keep pace with climate change. NMI revealed large portions of the breeding range shifting to novel climate conditions.

Main conclusions:
Our study demonstrates that forecasting climate vulnerability from ecological and evolutionary factors reveals insights into population-level vulnerability to climate change that are obfuscated when either approach is considered independently. For the Brown-capped Rosy-Finch, our results suggest that persistence may depend on rapid adaptation to novel climate conditions in a contracted breeding range. Importantly, we demonstrate the need to characterize novel climate conditions that influence uncertainty in forecasting methods.

| INTRODUC TI ON
Global climate change is dramatically affecting biodiversity, and extinction rates are accelerating across taxonomic groups (Urban, 2015). Alpine organisms that already inhabit the upper elevational reaches can be at particular risk of climate change driving upslope range shifts due to reduced potential to shift their range (Freeman et al., 2018;Sekercioglu et al., 2008); however, this risk may be tempered in regions that provide an abundance of microclimates (Seastedt & Oldfather, 2021). If range shift is not feasible, a species' long-term persistence in the face of climate change will likely depend on evolutionary or behavioural adaptation (Aitken et al., 2008;Forester et al., 2018;Hoban et al., 2016;Hoffmann & Sgró, 2011). Advances in ecological genomics are elucidating the genomic architecture of local adaptation (Hämälä & Savolainen, 2019;Savolainen et al., 2013;Tigano & Friesen, 2016) and providing insight into population-level responses to climate change Dauphin et al., 2021;Fitzpatrick et al., 2021;Fitzpatrick & Keller, 2015;Rellstab et al., 2016;Ruegg et al., 2018). While common garden experiments are widely recognized as the best method for identifying signals of local adaptation (de Villemereuil et al., 2016;Kawecki & Ebert, 2004), ecological genomic approaches provide an alternative in species where common garden approaches are infeasible due to constraints related to life history and/or conservation status (i.e. threatened or endangered status).
Ecological niche models (ENMs) are used to assess vulnerability to climate change by forecasting the distribution of climatic conditions that characterize an organism's current range (Guisan & Thuiller, 2005;Pacifici et al., 2015). Given the variety of terminology used in the literature surrounding ENMs, we will follow the guidelines set out by Sillero (2011) that ENMs model an organism's ecological niche and the resulting output maps forecast habitat suitability.
Genomic offset is a complementary approach to predicting climate vulnerability and provides a relative measure of the magnitude of evolutionary adaptation required for a population to track changing climate conditions (Capblancq et al., 2020;Fitzpatrick & Keller, 2015;Rellstab et al., 2021). Genomic offset is based on identifying geneticenvironment associations putatively underlying local adaptation and predicting future adaptive genetic composition based on the current genetic-environment associations (e.g., with gradient forest models; Fitzpatrick & Keller, 2015). However, genomic offset predictions may ignore key ecological factors (e.g., habitat suitability) that would affect persistence, especially for organisms with ranges that are experiencing drastic environmental changes. While genomic offset has predominantly been assessed independently of ecological factors (Capblancq et al., 2020but see Chu et al., 2021;Gougherty et al., 2021;Nielsen et al., 2021), vulnerability to climate change is a multifaceted problem that should be assessed with multiple methodologies and data sources (Dawson et al., 2011). Integrating genomic offset and ecological niche models would provide an understanding of the ecological factors shaping where populations could persist, and the evolutionary factors underlying the amount adaptation required to persist there.
The objective of our study was to combine methods for predicting population-level response to climate-driven disruptions to habitat suitability and genomic adaptation to improve forecasting of climate vulnerability. We addressed this objective using the Browncapped Rosy-Finch (Leucosticte australis), an alpine-obligate species endemic to the Southern Rocky Mountains (Wyoming, Colorado and New Mexico) and part of a broader species complex notable for specializations to alpine sky islands and arctic tundra (Johnson et al., 2020). While climate change broadly results in species shifting distributions poleward and upward in elevation (Chen et al., 2011;Parmesan & Yohe, 2003), the Brown-capped Rosy-Finch has limited potential for poleward range shift given the isolation of the Southern Rocky Mountains from other high-elevation mountain ranges and the presence of congeneric species that already inhabit those mountain ranges. Furthermore, Brown-capped Rosy-Finches already occupy nesting cliffs at the highest elevations (above 3350 m) of the Southern Rocky Mountains, which limits the possibility for major upslope range shifts, although they occupy lower elevations during the winter months (Johnson et al., 2020). Recent genetic studies have suggested that mountain ranges do not function as geographic barriers to dispersal for the North American Rosy-Finch complex (Black Rosy-Finch [Leucosticte atrata], Grey-crowned Rosy-Finch [Leucosticte tephrocotis], Brown-capped Rosy-Finch) given the level of ongoing gene flow among these species (Drovetski et al., 2009;Funk et al., 2021). Ongoing gene flow among Brown-capped Rosy-Finch populations may be an important component that mitigates genomic offset and prevents genetic isolation.
Here, we outline a process to assess climate vulnerability that considers evolutionary (e.g. genomic offset) and ecological factors (habitat suitability; Figure 1). We aim to answer the question: How can estimates of genomic offset and habitat suitability be combined to improve forecasts of climate vulnerability? Using genome-wide sequence data, we assessed population genetic structure and estimated levels of inbreeding and genetic diversity to describe spatial genetic variation and appropriately inform subsequent genetic-environment association (GEA) analyses Funk et al., 2019). We performed environmental variable selection to identify a subset of uncorrelated predictors for use in the GEAs and ENMs. We developed ecological niche models (ENMs) using the environmental predictor data and presence-absence data from the citizen-science database eBird (Sullivan et al., 2009). Additionally, we identified a subset of genomic variants associated with the environmental data and used

K E Y W O R D S
Brown-capped Rosy-Finch, climate change, conservation genomics, ecological niche models, genetic-environment association, genomic offset, Leucosticte australis these data to model allelic turnover across the landscape with gradient forest (Ellis et al., 2012;Fitzpatrick & Keller, 2015). Using ensembles of global climate models for two time periods, 2041-2070and 2071-2100(AdaptWest Project, 2021, we then forecast climate vulnerability in relation to genomic offset and habitat suitability. We demonstrate a novel application of the niche margin index (Broennimann et al., 2021) to highlight uncertainty in genomic offset predictions due to novel climate conditions. Specifically, this study aimed to (1) characterize the magnitude of genetic change required to track climate change and where populations could persist to minimize genomic offset; (2) predict climatedriven habitat suitability shifts into the future and (3) compare the underlying climatic drivers of, and spatial vulnerability to, genomic offset and habitat suitability. The integration of these approaches will provide a better understanding of evolutionary and ecological factors underlying species response to climate change and improve our ability to forecast climate change impacts on biodiversity.

| Field sampling and sequencing
We sequenced feather and blood samples from 116 individuals spanning 11 sites across the Brown-capped Rosy-Finch breeding distribution ( Table 1). Samples were collected during the breeding season of 2017 and 2018. Individuals from the Lost Man Lake and Independence Lake sites were combined as a single sampling unit for subsequent analyses based on their proximity (<1 km) and the low sample sizes (5 and 1 individuals, respectively). Engineer Mountain and Horseshoe Basin sites were also in close proximity (<5 km), but we retained them as separate sampling units due to the larger number of individuals per site (8 and 18 individuals, respectively).
We extracted DNA from blood samples using the standard protocol for Qiagen DNEasy Blood and Tissue Kits and we modified the protocol to maximize DNA yield from feathers. Whole-genome sequencing libraries were prepared following modifications of Science samples DMNS52416 and DMNS52417). The reference genome was annotated with the most recent zebra finch annotations available (NCBI GCA_008822105.2) using the program Liftoff (Shumate & Salzberg, 2021). For the input into all subsequent analyses, we extracted high-quality single-nucleotide polymorphisms (SNPs; Supporting information). F I G U R E 1 Workflow for combining ecological niche modelling and genomic offset for determining populationlevel climate vulnerability. Genomic data and occurrence data inform environmental variable selection by providing the geographic points for which environmental correlation is calculated. Species' life history information informs which variables are selected from correlated pairs, resulting in an uncorrelated set of environmental variables. Candidate adaptive SNPs are obtained through geneticenvironment association outlier analyses using the environmental variables and genomic data. The resulting candidate SNPs are the input into gradient forest models which predict adaptive genetic composition across the landscape. The subset of environmental variables are also used with occurrence data in ecological niche models to habitat suitability. Gradient forest models are used to predict adaptive genetic composition to future climate and the distance with the baseline environment provides the measure of genomic offset. Additionally, the niche margin index is calculated to quantify the extrapolation to novel climate. Ecological niche models are also projected to future climate to provide a measure of future habitat suitability. The integration of these models provides a description of where populations are most likely to persist in the future and the magnitude of genetic change required to persist there. Furthermore, regions of novel climate are depicted to highlight uncertainty in the forecasting method TA B L E 1 Sample location information with environmental data (1991-2020 period)

140.2
Note: ID = four letter abbreviation for sample locations used in the manuscript. Latitude and longitude specify the coordinates for sampling site locations and sample size specifies the number of individuals captured at these locations. MWMT is the mean temperature of the warmest month (°C), PAS is the annual amount of precipitation as snow (mm), SHM is the summer heat moisture index (calculated by dividing MWMT by the mean summer precipitation), and the last column is the elevation of the sampling site (m).
Pi = mean nucleotide diversity across 25000 base pair windows. F = individual inbreeding statistics. Ne = effective population size for locations with sufficient sample size for the linkage disequilibrium method of calculation. *Two sampling sites (Engineer Mountain and Horseshoe Basin) from which individuals were combined as a unit for analyses due to close proximity of the sites.

| Population genetic structure
We performed several analyses to describe geographic patterns of genetic variation. The presence of closely related individuals can skew signatures of population structure, so we used KING (Manichaikul et al., 2010) to identify and remove individuals with up to second-degree relationships (kinship > 0.0884). PCA provides an efficient nonmodel-based method for assessing population structure in high-dimensionality data sets (Patterson et al., 2006). We implemented principal components analysis (PCA) using the R package SNPrelate (Zheng et al., 2012) in R version 3.6.2 (R Core Team, 2019). Additionally, we estimated individual ancestry coefficients with the snmf function in the R package LEA (Frichot et al., 2014;Frichot & François, 2015), and tested a range of clusters from K = 1 to 6 with 100 iterations each.
Finally, we tested for effects of isolation by distance (linearized FST versus log10 geographic distance) with a Mantel test in the R package adegenet (Jombart, 2008).

| Bioclimatic variables
Snow is a major component of weather that shapes alpine communities. Snow cover can insulate soils from extreme cold air temperatures (Neuner, 2014) and also dictate the length of the growing season (Jonas et al., 2008;Keller et al., 2005). In some alpine plant species, reductions in snow cover can result in increased frost damage and decreased plant production (Abeli et al., 2012;Baptist et al., 2010;Inouye, 2000). The Brown-capped Rosy-Finch feeds on a variety of seeds throughout the year and on insects during the breeding season (Johnson et al., 2020;Martin et al., 1961;Packard, 1968;Warren, 1916

| Identifying putative adaptive variants
We used two genetic-environment association (GEA) approaches to identify a set candidate SNPs that are associated with environment.
First, we implemented the multivariate approach of redundancy analysis (RDA) as it performs well for detecting weak, multilocus signatures of selection . We performed RDA using Second, we used latent factor mixed models (LFMM) as a univariate regression model to identify candidate SNPs associated with each of the predictor variables (Frichot & François, 2015). We set the number of K latent factors based on the results from the individual ancestry coefficient results. For each model, we set the false discovery rate to 0.05 and calibrated the p-values by setting the genomic inflation factor to achieve a flat p-value distribution with a peak at 0 (François et al., 2016). LFMM analysis was conducted in R using the LEA package (Frichot & François, 2015). SNPs identified in both RDA and LFMM were used as the candidate SNP set putatively underlying local adaptation. We identified chromosomal position and gene information of the candidate SNPs using the Bedtools "closest" function (Quinlan & Hall, 2010) with the annotated Leucosticte australis genome. We identified candidate genes by selecting SNPs within 10,000 bases from genes of known function and tested for gene ontology enrichment with the chicken (Gallus gallus) genome using the Gene Ontology resource (Ashburner et al., 2000;Carbon et al., 2021;Mi et al., 2019).
Importantly, GEA analyses rely on the assumption that current allele frequencies are at equilibrium with the environment (Capblancq et al., 2020;Lasky et al., 2018). However, populations may experience an adaptational lag associated with historical environmental conditions (Browne et al., 2019). To test the influence of this assumption, we created two candidate SNP sets based on two temporal periods: one that temporally encompassed our sample period  and one based on potential adaptational lag ).

| Geographic distribution of putative adaptive variation
We used the gradient forest algorithm to describe the associations of spatial, environmental and genetic variables (Ellis et al., 2012;Fitzpatrick & Keller, 2015). Gradient forest is a machine learning method developed to model ecological community turnover in relation to environmental gradients by creating separate random forest models for each species (Breiman, 2001;Ellis et al., 2012).
Community turnover is then identified by aggregating environmental predictor importance for each species. This concept has been extended to landscape genetics by substituting allele frequencies at genetic loci for species and modelling adaptive genetic composition across the landscape (Fitzpatrick & Keller, 2015). The turnover functions in gradient forest allow for inference of the environmental predictors driving observed changes in allele frequency (Fitzpatrick & Keller, 2015). We fit gradient forest models to environmental and spatial data as predictors for the nine sampling sites with at least six individuals using the package gradientForest (Smith & Ellis, 2013).
We modelled adaptive genetic variation turnover on the landscape using the candidate SNP set as the response variable. Model tuning was performed on the parameters mtry (random subset of predictors used in random forest) and ntree (number of trees grown in each forest; Hastie et al., 2009). We evaluated model performance with prediction accuracy calculated from the out-of-bag samples (Ellis et al., 2012). We tested model performance of the candidate SNPs against a randomized model of candidate SNP allele frequencies and a SNP set that included putatively neutral loci (Supporting information). Using the top gradient forest model, we interpolated genetic composition across the remaining 1 km 2 cells from the breeding range for which we did not sample genetic data.

| Habitat suitability under climate change
We created ENMs using the ensemble modelling approach in the R package biomod2 (Thuiller et al., 2016; Supporting information).
Presence-absence data were obtained from the eBird Basic Dataset  Lek & Guégan, 1999). Given the focus of our subsequent analyses on temporal forecasting, we aimed to use a set of algorithms with balanced biases and avoided models that tend to project extreme outcomes (Beaumont et al., 2016).

| Genomic offset to climate change
Genomic offset estimates the magnitude of evolutionary adaptation needed for a population to keep pace with climate change (Capblancq et al., 2020;Rellstab et al., 2021). When using gradient forest models, genomic offset is calculated by the Euclidean distance between current genetic composition with the predicted genetic composition based on future environment (Fitzpatrick & Keller, 2015). We cal-

| Quantifying uncertainty in geneticenvironment associations
The niche margin index (NMI) is a metric that characterizes the distance from niche margins with 0 representing the margin, 1 being the maximum value within the niche, and decreasing negative values representing distance outside the niche (Broennimann et al., 2021). tive NMI indicate higher model uncertainty due to extrapolation in the gradient forest models.

| Population genetic structure
Whole-genome sequencing produced genomic data with an average 6x depth of coverage and variant filtering resulted in 429,442 SNPs for subsequent genetic analyses. We removed 12 individuals from the dataset due to relatedness. Visualizing PCA results revealed weak clustering of Pike's Peak individuals from other sampling sites ( Figure S1). The weak PCA clustering of individuals suggests low genetic differentiation among the sites, which was also supported by low pairwise FST values ranging from 0 to 0.042 (mean FST = 0.012; Table S1). The Mantel test did not identify as-  Table 1). The per-individual F inbreeding statistic was also similar across sampling locations (F mean = 0.11, range = 0.02-0.25; Table 1). Effective population size for the five sampling locations that had sufficient sampling size ranged from 108 to 403 (Table 1).

| Identifying putatively adaptive variants
The final uncorrelated environmental variable set consisted of mean temperature of the warmest month (MWMT), precipitation as snow (PAS) and summer heat moisture index (SHM; mean summer temperature divided by summer precipitation), as well as elevation. For RDA, we retained the first MEM (MEM1) spatial predictor for accounting for population structure as it was uncorrelated with the other predictor variables and explained 42.4% of the spatial variation. Model selection in the RDA retained all predictor variables. RDA outlier SNPs putatively associated with climate were identified by loadings on the first constrained axis ( Figure S4). We identified 2040 and 2045 candidate SNPs from the 1961-1990 and 1991-2020 environmental predictor datasets, respectively. In LFMM, we used a lambda of 0.7 to achieve the optimal distribution of p-values for each of the four predictor tests ( Figure S5). With K = 2 latent factors, we identified 4844 and 4502 candidate SNPs from the 1961-1990 and 1991-2020 environmental predictor sets, respectively. Intersecting the RDA and LFMM datasets identified 501 and 436 candidate SNPs for the  and 1991-2020 environmental predictor sets, respectively. Gene ontology enrichment analysis identified 12 genes associated with the biological process glutamatergic regulation of synaptic transmission (Gene Ontology ID: 0051966, p-value = 2.69e-6, false discovery rate = 3.67e-2) and six genes associated with regulation of small GTPase mediated signal transduction (Gene Ontology ID: 0051056, p-value = 2.73e-6, false discovery rate = 1.86e-2; Table S2).

| Geographic distribution of putative adaptive variation and habitat suitability
Our evaluation of tuning parameters in gradient forest models identified the out-of-bag testing accuracy to reach convergence with 100 trees (ntree = 100). Using all predictors in each tree (mtry = 5) achieved the highest proportion of variance explained across the predictors ( Figure S6). The comparison of the two time period predictor sets, with the corresponding candidate SNPs, revealed similar relative predictor importance ( Figure S7). Therefore, we continued all subsequent analyses with the 1991-2020 predictor set and candidate SNPs. With the candidate SNP set, raw predictor importance was ranked in descending order of precipitation as snow (

| Genomic offset and habitat suitability under climate change
The magnitude of genomic offset was highly variable across the

| Quantifying uncertainty in geneticenvironment associations
For the baseline time period, the majority of the geographic region for which we interpolated genetic composition was within or close to the niche margins derived from our sampling sites (Figure 4a). For F I G U R E 2 Performance of gradient forest models. (a) Raw R2 importance values for variables used as predictors in gradient forest model for three different datasets, which are: "All" is the total genomic variant set of 429,442 SNPs, "candidate" is the 436 candidate SNPs associated with the 1991-2020 baseline environment, and "random" is randomized genotypes of the candidate SNPs among the sampling locations.
Using the candidate SNPs, larger raw importance values were obtained with the environmental predictors (precipitation as snow, mean temperature of the warmest month, and summer heat moisture index) than in the other two SNP sets. (b-f) the turnover functions from gradient forest model show the weighted cumulative importance values, which represent the relative importance of a variable in explaining changes in allele frequency.
Here, only (b) precipitation as snow reveals consistently higher importance in the candidate SNP set than the other two datasets

| DISCUSS ION
In this study, we evaluate climate change consequences related to

| Comparing climate drivers of habitat suitability and local adaptation
For the Brown-capped Rosy-Finch, precipitation as snow, mean temperature of the warmest month, and elevation were the strong-  et al., 2021). Importantly, the upward elevational shift of predicted high habitat suitability may not necessarily correspond to a similar scale of actualized range contraction. A key factor that may mitigate climate change risks for alpine species, especially in the Rocky Mountains, is the highly heterogeneous topography of the alpine landscape (Seastedt & Oldfather, 2021). Alpine microtopography can result in thermal refugia along short horizontal distances that mimic air temperature changes of hundreds of meters upslope (Scherrer & Körner, 2010). The American Pika (Ochotona princeps) F I G U R E 4 Identifying the magnitude of climate shift to novel conditions. (a) Calculating the niche margin based on our sampling sites and the niche margin index to future climate revealed large portions of the breeding range shifting to novel climate conditions (purple). The southern portions of the breeding range had the largest geographic areas retaining similar climate conditions (green) to the sampling sites. (b-d) of the three environmental variables in the gradient forest model that change temporally (i.e. excluding elevation), the largest shifts to novel conditions are present in the mean temperature of the warmest month is an example of a small alpine species that can behaviourally adapt to suboptimal thermal regimes using different microhabitats (Millar et al., 2018;Rodhouse et al., 2017). While the thermal tolerance of the Brown-capped Rosy-Finch is unknown, behavioural adaptation to microhabitat use may be an important component of their climate change response. Given that Brown-capped Rosy-Finches nest in cliffs (Hendricks, 1977;Packard, 1968;Sclater, 1912), small changes in nesting site selection (e.g. cliff aspect) could provide dramatic differences in the microhabitat climate. Research into Rosy-Finch microhabitat usage and physiology would provide much needed additional information regarding predicted response to climate change.
The amount of precipitation as snow appears to have biological importance for both local adaptation and the realized niche for the Brown-capped Rosy-Finch (Figure 2b; Table S3). In alpine plant communities, snow cover can have a large effect on flower abundance and the evolution of adaptive traits to reduce frost damage (Inouye, 2000). In turn, this could affect Brown-capped Rosy-Finches foraging in the breeding season as they feed on available insects and seeds from a wide-range of plant species and families (Johnson et al., 2020;Packard, 1968 Branch et al., 2022). White-tailed Ptarmigan (Lagopus leucura) is another alpine specialist with low genetic differentiation and range-wide adaptive divergence potentially associated with diet (Fedy et al., 2008;Zimmerman et al., 2021). For Browncapped Rosy-Finch, further elucidating the connections between gene functions and local adaptation (e.g. linking genotypes and phenotypes) is an important next step in understanding the effects of climate change.

| Geographic patterns of climate vulnerability
Our Furthermore, these results underscore the importance of using multiple measures of vulnerability for informing conservation and management (Dawson et al., 2011;Rellstab et al., 2021). For organisms that inhabit regions experiencing large climate shifts, even the lowest genomic offset values may indicate relatively large allelic shifts required for a population to retain optimal genetic-environment associations.
Alpine climate conditions are changing dramatically in the Southern Rocky Mountains, especially in relation to snowpack (Pederson, Gray, Ault, et al., 2011; and summer temperature increases (Pepin et al., 2022).
In our study, NMI results suggest that the central and northwest portion of the Brown-capped Rosy-Finch breeding range are shifting to novel climate conditions (Figure 4a). Of the bioclimate variables most tied to habitat suitability and putative adaptive variation, the amount of precipitation as snow is decreasing across the breeding range in the future (Figure 4b), and the mean temperature of the warmest month is dramatically increasing (Figure 4c). However, our characterization of change in these specific bioclimate variables is based on the ecological niche model predictions of range from eBird citizen science data. Importantly, citizen science data for this organism may be more likely to be collected at lower elevations that are more accessible to observers than the higher elevation portions of the breeding range. This sampling bias could over-(or under-) estimate the current distribution of the breeding range, as well as the distributions of climate values across future time periods (Figure 4bd). Nonetheless, our NMI measures, which were solely based on climate distance from the climate niche defined by our sampling sites, reveal large climate shifts across high-elevation portions of the range (Figure 4a).

| Considerations in forecasting genomic offset
Recent reviews have highlighted a number of key assumptions and limitations that need to be addressed in the ongoing development of genomic offset methods for effective use in conservation (Capblancq et al., 2020;Rellstab et al., 2021). Genomic offset approaches assume that similar future conditions will result in similar genetic composition (space-for-time assumption). While this assumption may be problematic (e.g. multiple genetic architectures underlying an adaptive optimum), novel future conditions further increase the uncertainty of population response due to predicted genetic composition from unobserved environmental conditions. To address part of the uncertainty temporal extrapolation, we used the niche margin index highlight these regions of extrapolation to future climate conditions ( Figure 4a). Given the reliance of gradient forest methods on temporal extrapolation from non-linear turnover functions, potentially to novel conditions, we strongly recommend future studies to provide some measure of this uncertainty. Another similar assumption is that populations are at adaptive equilibrium with the temporal period during which genetic-environment associations are being tested.
Long-lived species (e.g. trees) are particularly prone to violate this assumption given that populations may have been established centuries ago with different selection pressures .
For the shorter-lived Brown-capped Rosy-Finch, we tested for the potential influence of adaptational lag by comparing environmental predictors between two baseline environmental periods. While our results suggested limited differences in genetic-environment associations with these two periods, additional study into the role of effective population size and genetic drift in adaptive (non)equilibrium in this this system may be insightful (Láruson et al., 2022).

| CON CLUS IONS
Here, we show that the Brown-capped Rosy-Finch faces climate threats across their breeding range from changing habitat suitability and disruptions of genetic-environment associations. Future persistence may depend on rapid adaptation to novel climate conditions in a contracted breeding range. Expanding future research to forecast climate threats across the fall and wintering range would facilitate an assessment of climate vulnerability across the full annual cycle. We also note the importance of identifying the potential for behavioural adaptation to alpine microrefugia that may mitigate climate change threats. The results of this study highlight the importance of combining multiple methods to characterize climate vulnerability in a more nuanced manner than provided by any of the methods alone.

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

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available at https://doi.org/10.5061/dryad.stqjq 2c4r. The analysis scripts are available at https://github.com/mgdes aix/bcrf-climate.

PEER R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/ddi.13628.