Genomic analyses point to a low evolutionary potential of prospective source populations for assisted migration in a forest herb

Abstract Climate change is increasingly impacting temperate forest ecosystems and many forest herbs might be unable to track the changing climate due to dispersal limitation. Forest herbs with a low adaptive capacity may therefore benefit from conservation strategies that mitigate dispersal limitation and evolutionary constraints, such as assisted migration. However, assisted migration strategies rarely consider evolutionary constraints of potential source populations that may jeopardize their success. In cases where climate adaptation is overshadowed by competing evolutionary processes, assisted migration is unlikely to support adaptation to future climates. Using a combination of population and landscape genomic analyses, we disentangled local adaptation drivers and quantified the adaptability and vulnerability to climate change of the self‐incompatible deciduous forest herb Primula elatior. Southern populations displayed a sharp genetic turnover and a considerable amount of local adaptation under diversifying selection was discovered. However, most of the outlier loci could not be linked to climate variables (71%) and were likely related to other local adaptation drivers, such as photoperiodism. Furthermore, specific adaptations to climate extremes, such as drought stress, could not be detected. This is in line with the typical occurrence of forest herbs in buffered climatic conditions, which can be expected to reduce selection pressures imposed by climate. Finally, populations in the south of the distribution area had increased sensitivity to climate change due to a reduced adaptive capacity and a moderate genetic offset, while central European populations were sensitive due to a high genetic offset. We conclude that assisted migration from southern source populations could bear significant risk due to nonclimatic maladaptation and a low adaptive capacity. Regional admixture and restoration of ecological connectivity to increase the adaptive capacity, and assisted range expansion to suitable habitat in the north might be more appropriate mitigation strategies.

constraints of potential source populations that may jeopardize their success. In cases where climate adaptation is overshadowed by competing evolutionary processes, assisted migration is unlikely to support adaptation to future climates. Using a combination of population and landscape genomic analyses, we disentangled local adaptation drivers and quantified the adaptability and vulnerability to climate change of the selfincompatible deciduous forest herb Primula elatior. Southern populations displayed a sharp genetic turnover and a considerable amount of local adaptation under diversifying selection was discovered. However, most of the outlier loci could not be linked to climate variables (71%) and were likely related to other local adaptation drivers, such as photoperiodism. Furthermore, specific adaptations to climate extremes, such as drought stress, could not be detected. This is in line with the typical occurrence of forest herbs in buffered climatic conditions, which can be expected to reduce selection pressures imposed by climate. Finally, populations in the south of the distribution area had increased sensitivity to climate change due to a reduced adaptive capacity and a moderate genetic offset, while central European populations were sensitive due to a high genetic offset. We conclude that assisted migration from southern source populations could bear significant risk due to nonclimatic maladaptation and a low adaptive capacity. Regional admixture and restoration of ecological connectivity to increase the adaptive capacity, and assisted range expansion to suitable habitat in the north might be more appropriate mitigation strategies.

K E Y W O R D S
assisted migration, climate change, forest herb, genetic offset, landscape genomics, local adaptation

| INTRODUC TI ON
Ongoing climate change will likely have increasingly severe effects on biological diversity (Butchart et al., 2010). Global surface temperatures increased by 1°C compared to preindustrial levels in the last decade (2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019), and land temperatures in Europe have increased even faster by on average 1.8°C (European Environment Agency, 2021). Due to these rapidly increasing temperatures together with increasing drought frequencies (Grillakis, 2019), more than half of the plant species in Europe are expected to become vulnerable or threatened by 2080 (Thuiller et al., 2005). Specifically, rapid climate-change induced shifts of their potential distribution will exceed their adaptive capacity or their ability to migrate to newly available habitat (Kubisch et al., 2013). The strong fragmentation of forest habitats on the European continent makes forest herbs particularly susceptible to climate change due to local genetic erosion and loss of adaptive potential, and due to their limited dispersal capacity (Dullinger et al., 2015;Naaf et al., 2021;Svenning et al., 2008;Van Daele et al., 2021). Assessment of their adaptive capacity is therefore required to predict species-specific vulnerabilities to climate change, and to devise mitigation strategies against climate change-induced local extinctions across their range (Bussotti et al., 2015;Razgour et al., 2019).
When plant populations experience increasing environmental stresses due to climate change and are unable to adapt or migrate, assisted migration could prevent their (local) extinction (Hoegh-Guldberg et al., 2008). Plant populations can be relocated through assisted population migration within their range (assisted gene flow) or outside the range (assisted range expansion) where the climate has become suitable. This requires the identification of suitable source populations and target sites. Source populations should have beneficial adaptive traits and sufficient genetic variation allowing future evolutionary responses (Aitken & Whitlock, 2013). Target sites for assisted gene flow are populations with a high sensitivity to climate change, that is, low adaptive capacity and/or a degree of maladaptation to the future climate (Capblancq et al., 2020;Holderegger et al., 2006;Rellstab et al., 2021). Furthermore, in case of assisted gene flow, careful consideration of the genetic relatedness between source and target populations is required to prevent potential outbreeding effects (Vandepitte et al., 2010), detect fitness trade-offs (Ågren et al., 2013), and avoid disruption of local adaptation to other environmental factors. Specifically, the introduction of genotypes adapted to environmental conditions varying substantially from those at the introduction sites increases the probability of maladaptation and can imperil evolutionary and migration potential (Wadgymar & Weis, 2017). As a result, such evolutionary constraints may impact the success of assisted migration aiming to facilitate adaptation to warming climates.
Local adaptation is shaped by divergent selection pressure, gene flow and demographic processes, which together drive phenotypic differentiation (Orsini et al., 2013;Savolainen et al., 2013). Local adaptation along geographical or environmental gradients often results in the gradual turnover of allele frequencies or phenotypes. In randomly mating species with large population sizes such clines are generally manifested gradually, whereas species with limited dispersal may be featured by a strong genetic turnover in specific regions of the gradient (Savolainen et al., 2013). The duality of the role of gene flow in local adaptation is that it can either increase genetic diversity and therefore increase the adaptive potential, or disrupt local adaptation by introducing maladaptive alleles (Akerman & Bürger, 2014;López-Goldar & Agrawal, 2021). When gene flow has been historically low, local adaptation may be more pronounced but genetic diversity could be eroded through genetic drift. Furthermore, there is the danger of outbreeding depression when individuals from locally adapted populations are introduced in more northern populations with the aim to complement them with genetic variants preadapted to warmer climates (Frankham et al., 2011). Forest herbs have been shown to frequently manifest distinct fine-grained signatures of local adaptation (De Kort et al., 2020;Garrido et al., 2012;Herrera et al., 2017), which may cause strong maladaptation to environmental factors other than climate, following assisted migration. Assisted migration of forest herbs, typically characterized by limited dispersal capacity and gene flow, thus requires careful population genomic analyses to pinpoint suitable source populations, but more sophisticated analyses are required to disentangle climate adaptation from other local adaptations and to evaluate sensitivity to climate change (Aitken & Whitlock, 2013;Vanhove et al., 2021).
Local adaptation to climatic drivers is species dependent (De Frenne et al., 2011), can occur at multiple scales (Csilléry et al., 2014;Pluess et al., 2016;Rellstab et al., 2017), and multiple adaptations in plant species can be related to distinct climatic drivers (Franks et al., 2014;Leroy et al., 2020;Mahony et al., 2020). However, the cost of acclimation to drought and temperature stress can result in fitness and metabolism trade-offs (Reich et al., 2003;Thiel et al., 2014;Vanwallendael et al., 2019). On a genetic level, distinct stressors generally result in different selection pressures on distinct gene groups (Lotterhos et al., 2018). The relation between adaptive genes and climate drivers gives information on existing climate adaptations, but when the climate changes it is necessary to know the mismatch between the current and required genomic composition to thrive under the novel conditions. Maladaptation to the anticipated future climate (i.e. genetic offset) can be determined as the shift in the adaptive genetic component required to match the expected climatic changes (Capblancq et al., 2020). Populations with high genetic offset thus require a high (genetic) adaptive capacity to deal with climate change, here defined as the climate-adaptive genetic variation which determines the potential of populations to change the phenotypic expression of functional traits to environmental change. When the genetic offset is integrated with the adaptive capacity it is possible to predict the sensitivity to climate change of specific populations. Knowledge of the adaptive capacity as well as the degree of maladaptation to climate is a major prerequisite for assisted migration to succeed. Maladaptation to environmental stressors other than climate is another major, but frequently overlooked, determinant of the success of assisted migration (Fitzpatrick & Keller, 2015;Rellstab et al., 2021). Even though some studies integrated the adaptive capacity, genetic offset to climate change, and signatures of maladaptation to environmental stressors in tree species (Jia et al., 2020;Martins et al., 2018;Rellstab et al., 2016), we are unaware of studies that evaluated these factors in forest herbs in order to assess the potential and risks of assisted migration and other potential mitigation strategies. This is essential as forest herbs might experience a reduced amount of evolutionary pressures to climatic factors due to the micro-climatic buffering capacity of forest ecosystems (Zellweger et al., 2020) and an increased amount of other local adaptation drivers (Baeten et al., 2015;De Kort et al., 2020;Van Daele et al., 2021).
Here, we aimed to assess signatures of selection (allele frequencies deviating from neutral expectations as a consequence of selective pressures imposed by the environment) and climate change sensitivity (as determined by the genetic offset and the genetic adaptive capacity) in Primula elatior, a self-incompatible herb species representative for European moist deciduous forests. Due to its limited dispersal capacity, P. elatior is unlikely to track the shifting climate (Honnay et al., 2002;Van Daele et al., 2021;Whale, 1983), rendering its persistence under climate change predominantly dependent upon its adaptive capacity (Bussotti et al., 2015). We sampled 29 Primula elatior populations along a large latitudinal gradient to generate single-nucleotide polymorphisms (SNPs) frequency data, using a genome-skimming approach (Wessinger et al., 2018). We aimed to

| Study species and data collection
Primula elatior subsp. elatior mainly occurs in oak or oak-hornbeam forests in sub-Atlantic and continental Europe and its distribution ranges from southern France to Northern Denmark (Leuschner & Ellenberg, 2017). P. elatior is highly dispersal-limited due to the absence of morphological seed dispersal adaptations and has specific germination requirements (Taylor et al., 2008). Its flowers exhibit heteromorphic reciprocal herkogamy and the species is selfincompatible (Keller et al., 2016). Pollen flow is highly limited and rarely exceeds 150 m (van Rossum et al., 2011).
To capture genome-wide signatures of selection, we sampled leaf material from 29 large populations (number of flowering individuals ranging from 430 to >1000 individuals) along a latitudinal gradient (c. 1645 km), largely covering the extent of the species' distribution (Figures 1 and 2). Leaf material from nine to 10 individuals at each site, spaced >100 m apart, were sampled and stored in silica gel, mounting to a total of 285 sampled individuals. Total genomic DNA was isolated from dry leaf tissue using a Plant DNA Extraction Kit (Norgen Biotek, 2015).
One mismatch or missing data point was allowed in the barcode read (5,343,298 ± 1,598,102 average total reads per sample ± SD). Sequencing adapter remnants were clipped from all raw reads and reads with a final length < 20 bases were discarded (2,671,608 ± 799,040 average quality trimmed read pairs per sample ± SD). Adapter-clipped reads were quality trimmed by removing sequencing errors and trimming was focused on the 3′-end to get a minimum average Phred quality score of 10 over a window of ten bases (2,512,362 ± 747,790 average adapter clipped read pairs per sample ± SD). Merged quality trimmed reads were error corrected (201,600,000 total read pairs) using 'Musket v.1.0.6' with a 21 k-mer size for correction (Liu et al., 2013). Furthermore, digital normalization of error corrected reads (148,503,886 total reads) were performed with the 'normalize_by_median.py' script from 'khmer v.1.1' with a kmer size of 32 and a coverage cut-off of 80 (Crusoe et al., 2015).
Quality trimmed reads were aligned against Arabidopsis as reference genome using 'BWA-MEM v.0.7.12' (Burrows-Wheeler Aligner; Li, 2013;Li & Durbin, 2009). Variant discovery and genotyping of samples was executed with 'Freebayes v.1.0.2-16' with a diploid setting and reads with more than two mismatches were excluded from the dataset (Garrison & Marth, 2012). Annotations of variant effects on annotated genes and transcripts were performed using 'SnpEff v.4.31' (Cingolani et al., 2012). The predicted genes and transcripts from the genome annotation were used to predict downstream functional effects of the variants.

| Data filtering, imputation and uncertainty
For the filtering of bi-allelic SNPs, a minimum read depth of 5 per SNP and a minimum allele frequency across all samples of 5% (Min. MAF = 0.05) were used as threshold. Furthermore, genotypes that were not observed in at least 10% of samples (i.e. at least 29 samples) were removed from the dataset. We selected SNPs with less than 50% missing data within each of at least 11 populations.
We made five subsets to maximize the number of SNPs, depending on the amount of populations in which a SNP with less than 50% missing data occurred: all 29 populations (3773 SNPs with 3.9% To determine library effects on downstream analyses, three increasingly restrictive SNP matrices, based on the SNP sequencing errors, were constructed (Appendix S1.1, Table S1). First, SNPs that were erroneously sequenced in two out of the three duplicated samples were excluded from the primary dataset (805 SNPs excluded, mean ± SE of the false discovery rate (FDR) = 21.8 ± 0.4). For a second dataset we used the SNPs with calls in only two of the duplicated samples and excluded SNPs that were erroneous in one of these sam- matrix against this putatively error-free dataset allowed assessing the impact of library effects on the results. We found that the genetic diversity and signatures of selection were very similar between datasets (Appendix S1.1, Figure S1) and therefore the SNP dataset with a FDR of 21.8 ± 0.4 was used for all downstream analyses.
The datasets were imputed for analyses that cannot handle missing data (Bayescan, RDA and sNMF and gradient forest). Specifically, the VCF data format was transformed to the genlight data format (containing alternative alleles 0;1;2) and the imputation was based on the rounded population mean. Some SNPs still contained missing data because they had a high missing data percentage on a population level but a low missing data percentage over the whole range (N = 279). The specific SNPs in these populations were then imputed based on a mean of three distinct regions (based on the sNMF analysis described in Section 2.3). The mean imputation uncertainty of the final dataset was 3.7% (Appendix S1.2, Figure S2, Tables S2 and S3).
Finally, we determined the FDR on a SNP level for each subsequent analysis, by repeating each methodology with distinct parameters (see Sections 2.4 and 2.5 for details), to estimate the cumulative uncertainty of climate-adapted SNPs.

| Population genomics
We combined two methods to identify loci that deviate significantly from background genetic structure, and are therefore assumed to reflect a locus-specific imprint of diversifying selection. First, using 'PCAdapt v.4.3.3' in R (Luu et al., 2017), the K number of principal components were selected based on scree plots and the amount of clustering when PC axes were compared. PCAdapt can avoid problems related to linkage disequilibrium (LD) with a clumping strategy, and less important SNPs were removed based on a window radius of 200 (distance between two SNPs on the same gene) and a squared correlation threshold of 0.1. The mean length of all predicted genes (245,806) was 309.8 ± 254.2 and 330.8 ± 335.6 in the final dataset.
A radius window of 200 removed most detectable LD and results were highly similar compared to larger window sizes (Appendix S1.3, Figures S3 and S4). When all five SNP subsets were evaluated with PCAdapt, a total of 3143 SNPs (28.2% of the filtered dataset) were discarded due to LD. Second, this dataset without LD was used in 'BayeScan v.2.1' (Beaumont & Balding, 2004), which relies on a Bayesian approach that estimates the posterior probability that a given locus is under selection, thus reducing false positives under a variety of demographic scenarios (Foll & Gaggiotti, 2008). Because we utilized multiple smaller subsets, we designated prior odds of 10 for the neutral model, which correspond to the outlier ratio of PCAdapt in subset 1 (11.5%). A total of 100 pilot runs with a length of 5000 were executed, followed by 100,000 iterations with a burnin length of 50,000 (Lotterhos & Whitlock, 2014). The FDRs were determined based on q values and ranged from 0.1% to 50%, with a precision of 0.1%. Deviation from neutrality can be determined based on the locus-specific effect (alpha parameter in Bayescan), where a positive alpha value indicates diversifying selection and a negative value indicates balancing or purifying selection (Foll & Gaggiotti, 2008). Considering our focus on climate change sensitivity, only outliers under diversifying selection, as determined by the locus-specific effect in Bayescan (alpha), were evaluated further.
SNPs under LD were removed based on the retained SNPs after the PCAdapt analysis. The number of genetic clusters (K) tested ranged from 1 to 25 for the first subset and 1 to 15 for subset 2-4. The optimal K was selected based on the minimum cross-entropy of each run. Each run was executed with 10 replicates per K, 200,000 iterations and a regularization parameter of 10. Statistical estimates of ancestry proportions were determined based on the sNMF models and depicted graphically.

| Environmental association analysis
To identify genetic variation that is linked to climatic variation while accounting for neutral genetic structure shaped by the spatial organization of our populations, we performed a series of redundancy analyses (RDAs) using the R package 'Vegan v.2.5-7' (Orsini et al., 2013). Specifically, a matrix of Hellinger-transformed minor allele frequencies for each population was analysed against climatic and geographic variables to identify potential drivers of adapta-  (Dray et al., 2021). We used a minimum spanning tree (gb-MEM) to select relevant Euclidian distance connections between populations and converted these to a neighbourhood matrix (Appendix S2, Figure S5). Incorporating the functional connectivity of dispersal modes in MEM, increases the biological realism of spatial vectors (Bauman, Drouet, Fortin, & Dray, 2018;Ver Hoef et al., 2018). Therefore, the cumulative landscape resistances  Table S4) relatively, and were therefore excluded for further analysis. The three most significant MEM (#4, #1 and #5 based on eigenvalue ranking) were included as spatial IBR variables in the RDA analyses (Appendix S2, Figure S6).
To select relevant climate variables we used forward and backward selection of partial RDA (pRDA) models (Borcard et al., 1992) with the MEM as conditional variables, calculated with the ordistep function of the 'Vegan v.2.5-7' package in R (Oksanen et al., 2013).
Potential explanatory candidate variables were the annual mean temperature (bio 1), max. Temperature of warmest month (bio 5), mean temperature of warmest quarter (bio 10) and precipitation seasonality (bio 15). The temperature variables (bio 1, 5 and 10) were correlated (r > 0.6) and were therefore evaluated separately with precipitation seasonality as second climate variable. Resulting models were highly similar due to the high correlation of temperature variables and only the model with the highest explanatory power (R 2 adj.), which contained max temperature of warmest month and precipitation seasonality as fixed variables (Appendix S2, Figure S7), were retained for further analysis (R 2 adj. Was 8.9 ± 2.7 for bio 5 compared to 8.5 ± 2.8 for bio 1 and 8.8 ± 2.6 for bio 10). Potential confounding effects between climate and geography were further evaluated with variation partitioning. The FDRs were determined based on q values and ranged from 0.1% to 50%, with a precision of 0.1%. Only outliers under diversifying selection, as detected by Bayescan, were evaluated further.
To obtain insights in the function of outlier loci, we conducted a test for enrichment of GO terms related to outliers (Bayescan, PCAdapt and pRDA) under diversifying selection with <5% FDR.
We considered significant GO terms of biological processes with a p value lower than 0.05, as detected by the Kolmogorov-Smirnov statistic with classic and elim algorithms (Alexa et al., 2006), and the Fisher statistic with classic and parent-child algorithms (Grossmann et al., 2007). A total of 264 genes in our gene universe had GO func-  (Zhang, 2021).
To calculate the genetic vulnerability to future climatic conditions we calculated the genetic offset with 'gradientforest v.0.1-32' in R, which has shown experimental support for genomic predictions (Fitzpatrick et al., 2021). The gradient forest algorithm calculates the genetic turnover, which reflects the magnitude of genetic distance of many putatively adaptive candidate loci along multiple environmental and geographical gradients (Capblancq et al., 2020;Fitzpatrick & Keller, 2015). The genetic offset reflects the overall genetic distance from the future genetic composition that is required to maintain the current gene-environment relationships (Vanhove et al., 2021). The inherent complexities involved in the calculation of the genetic offset require careful consideration of the potential caveats associated with gradient forest models (Láruson et al., 2022). Residuals were normally distributed and no visual heteroscedacity was detected.
To evaluate the sensitivity to climate change of P. elatior it is necessary to take both the adaptive capacity and the genetic vulnerability into account (Vranken et al., 2021). To this end, a climate sensitivity metric was constructed for each scenario that corrected the genetic offset for the expected heterozygosity of partial RDA outliers in each population: This formula gives a similar weight to both the expected heterozygosity and genetic offset as the range of both metrics are highly similar (0-0.3). The relation of climate sensitivity to latitude was modelled with a linear model and a quadratic term (genetic offset ~ latitude + latitude 2 × scenario + year). Residuals were normally distributed and no visual heteroscedacticity was detected.

| Outlier analysis and population structure
With a FDR threshold of 5%, PCAdapt detected 266 outliers (3.4%) under diversifying selection, while 163 outliers (2.1%) were detected with a FDR below 1%, and 107 outliers (1.4%) with a FDR below 0.1%. A large amount of outliers were thus featured by a decisive signature of diversifying selection (Figure 3). Along the first PC axis three clusters could be identified with a clear genetic differentiation between southern, central and north European populations (Appendix S3, Figure S8). Inference of the population structure through sNMF analysis resulted in five ancestral clusters when all populations were considered (subset 1). SNPs that did not occur in all populations (subset 2-4) were aggregated in 14 distinct clusters but did not yield additional insights. Southern populations had a low amount of mixture but with three abrupt cluster transitions between populations ( Figure 2). Northern populations, on the other hand, had a large amount of mixture and smooth transitions between clusters ( Figure 2, Appendix S3, Figure S8).

| Unravelling environmental drivers of local adaptation
The included isolation-by-resistance MEM with the strongest contribution (MEM 4, R 2 adj. = 2.8%, λ = 3.46, p = 0.001) had strong regional eigenvector value turnover in south Europe and central Europe (similar to the dark blue and red genetic turnover in Figure 2; Appendix S2, Figure S6). The second most contributing eigenvector (MEM 1, R 2 adj. = 2.8%, λ = 3.9, p = 0.001) was mostly related to regional eigenvector value turnover in north Germany, and the final included eigenvector (MEM 5, R 2 adj. = 2%, λ = 3, p = 0.001) was related to regional eigenvector value turnover in south-central Europe (similar to the edges of the red cluster in Figure 2).
Climate had more explanatory power on SNP frequency variance than eigenvectors based on IBR in both the RDA and partial RDA models, and the overlap in explanatory power was minimal ( Table 1). The partial RDA detected 63 outliers (0.8%) with a FDR below 5%, which were thus related to climatic clines, and 53 outliers (0.7%;

Most of the variation in
Appendix S4,   Figure S10).

| Sensitivity to climate change across a latitudinal gradient
The genetic diversity (H e ) of SNP frequencies under diversifying selection driven by climatic clines as predicted by our general linear model (estimated marginal means ± SE = 0.23 ± 0.007) was significantly lower (Z ratio = 4.6, p < 0.001) than the neutral genetic F I G U R E 3 The partial redundancy analysis (pRDA) biplot in panel a illustrates the pRDA scores of the first two axes with Primula elatior single nucleotide polymorphisms (SNPs) as response variables, climate as fixed explanatory variables and isolation by resistance (IBR) moran eigenvector maps as conditional explanatory variables. Decisive outliers had <1% false discovery rate (FDR), very strong outliers <5% FDR, strong outliers <10% FDR and substantial outliers <25% FDR. Only decisive and very strong outliers were used for further analysis. The pRDA biplot in panel b displays the pRDA scores of the sampled population sites along the distribution range in Europe. Numbers indicate the plot ID from south (low) to north (high). The latitude colour scheme indicates the latitude in metres × 10 6 (Coordinate system LAEA89). The statistical evaluation of the pRDA model (Climate | IBR) can be found in Table 1. Our gradient forest algorithm predicted strong genetic turnover of climate related SNP outliers around the 22°C maximum temperature of the warmest month threshold (light blue to dark blue crossover in Figure 6), with three abrupt transitions in southern Europe (Figure 6, Appendix S3, Figure S8). The genetic offset TA B L E 1 Results of the (partial) redundancy analyses (RDAs) to partition among-population genetic variation of Primula elatior into climate, isolation by resistance (IBR), and their combined effects.  The climate outliers of the eight most southern populations had a relatively low genetic offset (Figure 7, red to green in Figure 6), but also had significantly reduced genetic diversity for climate outliers ( Figure 5). Consequently, the overall sensitivity to climate change (as determined by genetic diversity and genetic offset) was high in the south and only starts to decrease in northern France (R 2 = 83.3%).
The sensitivity to climate change was predicted to be significantly indicates that the integrated IBR patterns were good indicators for the demography and gene flow (Hoban et al., 2016), which explained 5% of allele frequencies ( Table 1).
The strong genetic differentiation along the range (Figure 2) favours the use of local genetic material for restoration purposes (Breed et al., 2018;Broadhurst et al., 2008). However, rapid climate change requires considerable adaptive shifts that may not be achieved without assisted gene flow. Here, we found considerable signatures of diversifying selection across the latitudinal gradient, pointing to the occurrence of alleles preadapted to a warming climate somewhere along the gradient. Climatic clines explained a considerable amount of genetic variation between populations (13%), but only a small amount of SNPs were significant outliers under diversifying selection that were likely to be driven by climatic selection F I G U R E 7 The genetic offset and climate sensitivity of Primula elatior climate outlier SNPs (Climate | IBR; Table 1; Figures 3 and 4) as predicted by the gradient forest algorithm. The genetic offset reflects the overall Euclidian genetic distance between the current adaptive genetic turnover ( Figure 6) and the predicted adaptive genetic turnover in 2050 (left)  pressures (0.7%). Furthermore, most of the outliers were biological functions related to metabolic processes and no adaptations to drought or heat stress were detected. However, it has to be noted that the stringent selection procedure, which was used to avoid uncertainties in the data analyses, could have eliminated potentially informative alleles on drought and heat resistance. Complementary sequencing techniques with higher read depth, preferably combined with common garden experiments, could potentially uncover overlooked drought resistance signatures, which are often characterized by a complex polygenic architecture. Nevertheless, growth regulation in response to seasonal factors is known to be important in perennials (Wingler, 2015) and perennial herbs generally have a greater capacity for temperature acclimation than annuals (Yamori et al., 2014). Temperature regulation results in altered carbon dynamics and generates sugar signals that further modulate metabolic pathways involved in biosynthetic and catabolic processes. These metabolic processes can in turn result in local adaptations of relative growth, caused by resource allocation trade-offs in acclimation regulation processes (Gray & Brady, 2016;Wingler, 2015 (Yamori et al., 2014).
P. elatior is highly dependent on the short light phase before tree canopy leaf development in spring (Baeten et al., 2015;Taylor et al., 2008) and photosynthetic adaptations may have a strong relation to phenological adaptations and growth regulation (Poorter et al., 2009;Rothstein & Zak, 2001). However, further research is needed to disentangle the intricate relationships between plant metabolism, photosynthesis, and local adaptation to climatic factors.

| Climate sensitivity
Forest herbs lacking adaptations for long-distance dispersal are estimated to lose genetic diversity at a faster rate under an increasingly changing climate (Alsos et al., 2012). Without adaptive evolution, Primula elatior is expected to lose over half of the total distribution area by 2050 due to climate change and dispersal limitation (Van Daele et al., 2021). Furthermore, southern populations are expected to be most affected and therefore populations will be highly dependent on their adaptive capacity. The potential of P. elatior to adapt to the changing climate will depend on the genetic offset of climate outliers and the adaptive genetic diversity, together shaping its climate sensitivity. We found that the cumulative importance of climate outliers sharply increased between a max. Temperature between 22°C and 23°C, which corresponds to a sharp increase of the adaptive genetic turnover in Central France. The sharp increase in adaptive genetic turnover is likely related to a combination of climate selection pressures (Appendix S2, Figure S7) and low gene flow in these regions (Appendix S3, Figure S9). When environmental change increases, gene flow is needed to introduce preadapted alleles and enable local adaptation under the changed conditions (Blanquart et al., 2013). However, the expected climate shift in this region resulted in a high genetic offset in central Europe and the low migration rate due to dispersal limitations are unlikely to alleviate climate change effects (Van Daele et al., 2021). Overall, southern and central European populations are more sensitive to climate change than northern populations. It is however important to take into account that the appropriate application of the gradient forest algorithm is still under active development and that the genetic offset could be confounded by the polygenic nature of drought resistance (Láruson et al., 2022). Furthermore, the potentially incomplete depiction of the genomic architecture and unverified relationship between the adaptive genetic offset and the adaptive capacity could have influenced the climate sensitivity curve. Further research is needed to evaluate the empirical relationship of the calculated climate sensitivity and plant fitness under climate change. Nevertheless, the low adaptive capacity, likely maladaptation of southern genotypes to longer photoperiods and other local conditions, and potential outbreeding effects could jeopardize successful assisted migration from southern to central and northern European regions. Therefore, conservation measures should prioritize preserving or improving (meta-) population stability trough ecological restoration of the habitat quality and ecological connectivity. Additional measures could include admixture provenancing within population clusters, characterized by low genetic turnover and a wide selection from various environments ( Figure 2; Breed et al., 2013), to improve (south/centre) and maintain (north) the adaptive capacity across the range (Vergeer et al., 2004;Whiteley et al., 2015). After careful monitoring of the admixture provenancing, it could then be decided to introduce offspring into more central European populations to alleviate their high genetic offset. Finally, the high adaptive genetic diversity in northern populations and the low climate offset could enable successful range expansions to projected suitable habitat in Scandinavia. As many perennial forest herbs exhibit similar dispersal modes and life history strategies (Verheyen et al., 2003), experience buffered climatic conditions (Zellweger et al., 2020), and are generally sensitive to photoperiodism (Flynn & Wolkovich, 2018), it is likely that evolutionary constraints may impact the efficacy of assisted migration across forest species. Further research can shed light on the general tendency of maladaptation and climate sensitivity in forest species.

ACK N OWLED G M ENTS
This research was funded by the Flemish Research Foundation (FWO project G091419N). Furthermore, we would like to thank Koenraad Van Meerbeek and Stef Haesen for providing microclimate data in the preliminary analysis.

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

DATA AVA I L A B I L I T Y S TAT E M E N T
Genomic sequencing data (VCF) and georeferenced sampling locations are deposited in a Dryad data repository named 'Genome skim-