The role of dispersal limitation and reforestation in shaping the distributional shift of a forest herb under climate change

Forest herbs might be unable to track shifts in habitat suitability due to rapid climate change and habitat fragmentation. In this study, we quantified the role of dispersal limitation and the potential mitigating effect of large‐scale reforestation on the redistribution of the herbaceous forest plant species Primula elatior under climate change.


| INTRODUC TI ON
Global temperature and precipitation patterns are rapidly changing and are expected to severely affect the distribution of many plant species . Land temperatures in different European regions are projected to increase further by 1.4 and 4.2℃ under Representative Concentration Pathway (RCP) 4.5 scenario and by 2.7 to 6.2℃ under the RCP 8.5 scenario (by 2071-2100, compared to 1971-2000) (European Environment Agency, 2020). Furthermore, climate models have shown a high probability of increasing intensity and duration of droughts by the end of the century (Cook et al., 2020). As a consequence, a general latitudinal range shift of plants is expected from south to north .
Europe is among the most forest-rich regions in the world, with forest covering 34% (EU27) of its total land area. However, forests and their understorey vegetation are heavily impacted by anthropogenic influences, increasing their sensitivity to global changes (European Commission, 2011). The understorey in temperate forests is often species-rich, provides food and shelter for many other organisms and mediates ecosystem processes such as tree regeneration, nutrient cycling and evapotranspiration (Landuyt et al., 2019). Forest herbs require very specific habitat conditions and are often highly dispersal limited. Many forest herbs will therefore be unable to track climate-induced shifts in their habitat distribution (Honnay et al., 2002). In-depth modelling of how global change drivers affect specific forest herbs on a local, regional and continental scale have the potential to guide mitigation strategies to safeguard the diversity and functioning of temperate forests. However, integrated modelling approaches that predict changes in the habitat distribution of forest herbs under global change, and at the same time assess dispersal routes are challenging and require improved predictive capacity.
Species distribution models (SDMs) have provided valuable insights into range shifts of certain plant species under climate change scenarios (Guisan & Thuiller, 2005). Yet their accuracy can be limited because input data often do not cover the entire current extent of a species´ range and because predictions have too coarse resolutions. Furthermore, inferences in terms of local extinction and colonization are generally based on range shifts alone (Chave, 2013;Fordham et al., 2012;Guisan et al., 2017). To reliably use SDMs for biological conservation, ecological restoration and impact mitigation, it is necessary to use a predictor resolution and response accuracy which accurately reflect the species' ecological niche breadth (Connor et al., 2018), ideally covering its whole range (Lembrechts et al., 2018;Thuiller et al., 2015). This is especially true for habitat specialists such as herbaceous forest plant species, where the potential for establishment is determined by a complex interplay between biotic and abiotic interactions. Integrating landscape-scale variables that partly capture topoclimatic processes, such as fine-scale topography (Lassueur et al., 2006), land use classes (Sirami et al., 2017) and the presence of riparian zones (Moore et al., 2005), in SDMs may therefore considerably increase their accuracy (Clerici et al., 2013).
Apart from their specific habitat requirements, many forest herb species are also known to be strongly limited in their dispersal, hampering the colonization of new suitable habitat in areas where habitat fragmentation has surpassed a critical level (Ehrlén & Eriksson, 2000). To accurately assess whether habitat fragmentation impedes future range expansion, it is essential to assess the dispersal potential between currently occupied habitat patches and future suitable habitats (Bateman et al., 2013;Franklin, 2010).
Recent advances in dispersal modelling allow for assessments of species-specific dispersal modes through genetic optimization of landscape connectivity (Peterman et al., 2019). Optimization methods rely on the expected relationship between landscape connectivity and genetic distance for selecting the most likely dispersal scenario. Whereas these optimization methods have successfully been applied to animal populations (Littlefield et al., 2017;Maiorano et al., 2019), they have rarely been used for habitat connectivity analysis in plants (Dickson et al., 2019). This is remarkable, because the importance of seed dispersal limitation in mediating range shifts is well known for specialist plant species (Ehrlén & Eriksson, 2000;Honnay et al., 2002;Skov & Svenning, 2004). Furthermore, accurate dispersal modelling of plant species is important in order to guide mitigation strategies such as assisted migration (Lunt et al., 2013) and connectivity restoration (Krosby et al., 2010) in the coming decades.
Changes in forest cover and composition will affect understorey habitat availability and may have a large impact on the dispersal and colonization potential of forest herbs undergoing climate change (Canadell & Raupach, 2008;Millar et al., 2007). Furthermore, future changes in forest cover and composition will directly affect habitat availability for understorey species (Nieto-Lugilde et al., 2015).
Reforestation policy targets, as set by the European Commission, and natural regeneration on abandoned farmland are therefore expected to strongly affect the extent of species distribution shifts induced by climate change (European Commission, 2019;Guo et al., 2018).
Although some studies have examined climate change effects on afforestation success (Duque-Lazo et al., 2018) and tree species distributions (Buras & Menzel, 2019), we are unaware of any studies that have integrated afforestation efforts in modelling the distribution of forest understory species under climate change scenarios.
Primula elatior was used as a study case because it is widespread and highly representative for other forest herbs from alluvial deciduous forests. Due to its specific habitat requirements, intolerance to desiccation, absence of seed dispersal mechanisms and selfincompatibility, it may be specifically sensitive to climate change and habitat fragmentation (Honnay et al., 2002;Taylor & Woodell, 2008).
Here, we aimed to assess how climate change and afforestation/reforestation efforts would affect the amount of accessible habitat for  (Taylor & Woodell, 2008). Successful introductions of P. elatior have occurred in Norway, which could suggest that the current range is restricted by dispersal limitation and not by the cold edge of its thermal niche (collected samples; Alm & Often, 2009).
Primula elatior typically occurs in sub-Atlantic and medio-European oak or oak-hornbeam forests of the alliance Carpinion betuli and to a lesser extent in alluvial forests with European alder and ash of the alliance Alno-Padion (Hennekens et al., 2010;Leuschner & Ellenberg, 2017). The species is known to have limited colonization capacity of newly established forests .

| Species occurrence data and pseudoabsence generation
From a total of 26 species occurrence databases (Appendix S1), we selected occurrences that were (a) only recently documented (post-2000), (b) of high precision (<100 m accuracy), and (c) from natural environments (no garden escapes or introduced individuals). The resulting 24,267 records were partitioned into spatially independent training and validation data, with a distance of at least 5 km between occurrences of the training dataset and occurrences of the validation dataset. Furthermore, each dataset was spatially thinned with a 10 km minimum distance between observations to avoid spatial autocorrelation. This distance was based on semi-variograms for each explanatory variable (Aiello-Lammens et al., 2015;Veloz, 2009). This resulted in seven training and validation subsets of 1,579 and 790 occurrences, respectively (see Figure 1 for an overview of the workflow). Filtering of available true absences (<100 m coordinate accuracy and post-2000 observation) resulted in a spatially biased dataset (Chytrý et al., 2016), and therefore, a presence-background approach was used. Background data were generated with spatial profiling (or geographic restriction) to (a) avoid false negatives, (b) exclude background data from areas outside of the species range, and (c) increase sensitivity (or the true-positive rate) of the modelling procedure (Barbet-Massin et al., 2012). More specifically, a spatial profiling range was determined ranging from 2.5 km distance from observed occurrences to 24.3 km, which was the maximum observed distance between occurrences in Europe. Spatially independent (10 km) background points were generated with 50%, 100% and 200% (maximum

| Explanatory variables
A total of 13 independent variables (Spearman r < .6) were included at 100 m resolution in the modelling procedure (Table 1; Appendix S2).
Land use, elevation, aspect, degree of soil wetness and a river network map (Copernicus Land Monitoring Service et al., 2020) were included to model the effect of landscape-scale properties on the current species distribution (Lembrechts et al., 2018). The degree of soil wetness was based on the Water Wetness Probability Index (WWPI) and was determined independently from the vegetation cover. The river network database was transformed to an Euclidian distance raster, representing the distance of each cell (in metres) to the river features, and distances were log 10 -transformed. The landscape-level tree cover (originally at 1 km 2 resolution) of the three most important co-occurring tree genera (Quercus, Carpinus, Fagus) were also included (Brus et al., 2011). Physical properties of the soil (n = 5; Hiederer, 2013) were reduced in dimensionality with a Principal Component Analysis (PCA) to obtain two independent soil Principal Components (   was based on the predictive capacity (AUC/AICc) of non-collinear combinations (r < .6; 1,077 models; Appendix S3), and three ecologically relevant and independent climatic variables were included in the final model: temperature seasonality (Bio 4), precipitation of wettest month (Bio 13), and the precipitation of warmest quarter (Bio 18). Climate data, landscape-level tree species cover and soil data were disaggregated to 100 m resolution with bilinear interpolation to match the resolution of the other datasets (Fordham et al., 2012).

| Model building and projections
Nine distinct algorithms and ensemble models were evaluated on the extent of the Netherlands and a maximum entropy algorithm (Maxent) yielded the best results based on AUC. True absence Maxent models were compared with a presence/background strategy on the same extent and the latter performed better based on AUC and AICc (Appendix S4). The selection of nonlinear functions of environmental variables (features) and the settings of the regularization multiplier, which controls overfitting/complexity by regulating the chosen feature class intensity, are known to have a strong effect on Maxent predictions and should be carefully chosen (Merow et al., 2013). Parameterization of feature classes (linear, product, quadratic, hinge and threshold) and regularization multipliers were evaluated by ranking models ( features were selected with a default (1) regularization multiplier as optimal parameters. The use of a default regularization multiplier is justified because we had a large amount of records, a diversity in the environmental data over the whole range, and less complex features than the default (Phillips & Dudík, 2008).
The resulting Maxent models were used to predict the relative suitability (presences relative to background data) for P. elatior in Europe for the current and projected climate (RO I), and for a forest restoration scenario (see section "Assessing mitigation effects of envisioned large-scale afforestation"). Bioclimatic projections according to three greenhouse gas scenarios in 2050, namely RCP 2.6, RCP 4.5 and RCP 8.5, were used to analyse future biogeographical shifts of the species distribution. RCP variables were based on averages from 11 general circulation models (Appendix S2).
Maxent's relative suitability output was transformed with a complementary log-log (cloglog) function to an estimate of occurrence probability ranging between 0 and 1 (Fithian et al., 2015;Phillips et al., 2017), and the resulting projections based on the 7 subsample datasets were averaged to obtain robust results. To determine a habitat suitability threshold for delineating suitable patches in an otherwise unsuitable landscape matrix (Nenzén & Araújo, 2011), we used the average projected probability (meanPred; see Appendix S5 for explored alternatives) of all 24,267 occurrences in our dataset (0.714; 95% CI [0.712, 0.715]). The meanPred threshold has been advised to reduce the false positive rate when species prevalence is moderate or high, which is essential to model dispersal between patches (Liu et al., 2013). All habitat above this threshold was considered suitable and patches were delineated based on connected cells according to the 8-neighbours rule (queen's case) with "landscapemetrics v. 1.5.1" (Hesselbarth et al., 2019). Thus, habitat patches are here defined as areas of suitable habitat consisting of contiguous cells with a minimum habitat suitability.
To validate the binary classification of patches, the true skill statistic (TSS) was calculated, which is based on the true-positive rate (TPR; within the habitat patch or an edge raster cell) of all occurrences and the true negative rate of a matching random sample of 24,267 absences (mean of 10 random absence generations). The TSS (TPR+TNR-1) ranges from −1 to 1, with values <0 indicating a performance which is no better than random (Allouche et al., 2006).

| Modelling dispersal and migration events to determine the accessible area
We modelled potential dispersal and migration events based on circuit theory, which describes the landscape as a resistance matrix (cells on a geographical raster grid) and which determines multiple pathways of low resistance between suitable habitat patches (nodes). This measure of isolation between pairs of nodes was then expressed as effective resistance (McRae et al., 2008). Dispersal was here defined as seed movement between habitat patches within a specific scenario, thereby allowing functional connectivity between habitat patches. Migration was defined as repeated seed dispersal events across successive generations, thereby allowing functional connectivity between a current habitat patch and a habitat patch that will be present in the future (2050). The future patch accessibility was determined by both the potential for migration and the potential for dispersal between projected habitat patches (future isolation). The current patch accessibility was determined only by the dispersal potential (current isolation) under the assumption that the current projected distribution is stable and not yet affected by climate change. Because land use, distance to rivers and elevation were expected to affect seed dispersal probabilities in P. elatior, these variables were used to define the resistance matrix. For example, rivers were expected to facilitate dispersal and thus decrease the effective resistance, while anthropogenic land use and high elevation (upslope) was expected to increase the effective resistance.
For elevation, the assumption was made that low resistance matrix values at low elevation were a proxy for downhill dispersal, and high resistance matrix values at high elevation were a proxy for uphill dispersal. to exclude populations in non-equilibrium states. We used single nucleotide polymorphism (SNP) skimming (Wessinger et al., 2018) to construct SNP matrices containing 406 pairwise Nei's genetic distances (Nei, 1978). The SNP matrix was constructed with the R packages "adegenet v.2.1.2", "vcfR v.1.8.0", "tidyverse v.1.2.1" and the genetic distances were calculated with "stAMPP v.1.6.1". These 29 populations captured variation in land use and river distance but not in elevation. To nevertheless account for elevation as a potential dispersal-inhibiting factor, we retrieved genetic data from Konečná et al. (2019), who evaluated colonization of subalpine habitats in central European mountains by P. elatior. This dataset contains a pairwise F st matrix for 16 populations or 120 pairwise distances and was used to separately optimize the resistance matrix of elevation. The genetic optimization (232 models; Appendix S7) was based on the log-likelihood rank of standardized maximum-likelihood populationeffects (MLPE) mixed models, using genetic distance as the response and log 10 -transformed effective resistance as the predictor (Clarke et al., 2002). Additionally linearized partial Mantel tests were used, removing the effect of log 10 -transformed Euclidean distance from the log 10 -transformed effective resistance (Cushman et al., 2006).
The distinct genetically optimized resistance matrices (land use, distance to rivers and elevation) were parameterized by iteratively assigning weights (from 0.1 to 1) to each resistance matrix (110 models), resulting in a new resistance matrix which was then used for connectivity analysis (Koen et al., 2012). The resistance matrices for the projected forest restoration scenarios were adapted to reflect the genetically optimized land use resistance matrix (Appendix S7).
We modelled whether dispersal between habitat patches could take place as a function of the effective resistance between them.
Because suitable habitat patches can either function as sources for further dispersal (stepping stones or source patches) or dispersal dead-ends, we also quantified the within-patch dispersal potential for each suitable patch based on the effective resistance between patch edge and patch centroid. Habitat patches where within-patch dispersal from edge to centroid is highly unlikely (dispersal dead-ends) were then excluded as source patches or stepping-stones. Migration potential was then quantified based on the effective resistance between currently suitable source patches and projected suitable habitat patches in 2050 (see Appendix S8 for the effect of source patch exclusion on migration probabilities). The potential for dispersal (between and within habitat patches) and migration was based on the isolation by resistance relationship. The sigmoid isolation by resistance relationship (Figure 2) was determined by a self-starting nonlinear least square model with the genetic similarity as response and the logarithmic effective resistance as predictor. Genetic similarity was defined as a function of Nei's genetic distance (y) and ranges between 1 and 100: The isolation by resistance model was used to predict genetic similarities based on the calculated effective resistances for dispersal Genetic similarity = 1 − (y − min (y)) max (y) − min (y) × 100 F I G U R E 2 Sigmoid isolation by resistance relationship between the genetic similarity and the log 10 -transformed effective resistance (see Appendix S7 for the untransformed relationship with genetic distance) of Primula elatior in Europe. Genetic similarity is a function of Nei's genetic distance and ranges between 1 and 100 with high values indicating high similarity. This isolation by resistance relationship determines the dispersal potential, migration potential and the accessibility of habitat patches. Dispersal limitation takes place at effective resistance values where the slope starts to decrease (log 10 -transformed effective resistance of approx. 1.25). The threshold (dashed red line) is defined at the centre of the sigmoid curve (inflection point) where the decrease in genetic similarity is at its highest point and dispersal limitation is certain. The linear decrease of genetic similarity at the higher end of the log 10 -transformed effective resistance (right side of the x-axis) suggests a migration-drift equilibrium (Van Strien et al., 2015) Threshold:

| Assessing shifts in patch distribution and configuration
To assess shifts in the patch distribution, we evaluated the me- Patch accessibility is a binary metric determined by the dispersal and migration potential, as explained in the last section. The PRI is an index accounting for patch size, similar to the patch proximity index proposed by Gustafson and Parker (1994), but with nearest-neighbour patch distance replaced by the effective resistance. The PRI is an indicator of metapopulation dynamics, with high values indicating stable metapopulations (Gustafson & Parker, 1994). The PRI value of a habitat patch was calculated by identifying each habitat patch (i) whose edge lay at least partially within the proximity buffer (5,000 m) of the focal patch centroid being indexed. The PRI was calculated based on the area in hectares (A) and the edge-to-edge effective resistance of dispersal (R) from patch i to its nearest neighbour: RCP Scenario represents a categorical variable indicating the respective climate change pathway (3) & Zimmermann, 2000). The explained variance of each fixed term (scenario, geographic components and their interaction) was partitioned to disentangle climate change effects from geographic effects. A more in depth analysis of explained variance drivers, regional effects, the relative contribution of area and dispersal probabilities to PRI, and the relative contribution of dispersal and migration probabilities to patch accessibility can be found in Appendix S8.

| Assessing mitigation effects of envisioned large-scale afforestation
To analyse the potential mitigation effect of anticipated forest restoration targets in the EU, we constructed a land use change sce- and an intersect of each separate layer was taken with the LUISA forest reference scenario, which depicts the projected land use change by 2050 (Lavalle, 2014 Mitigation was a binary categorical variable that determines whether the species distribution model for 2050 was created using the baseline land use or using the reforestation land use scenario.

| Drivers of the habitat suitability and distribution
The predicted suitability for P. elatior is high from the foot of the Pyrenees and throughout the Atlantic zone, but decreases towards Denmark (Figure 3). Furthermore, suitability hotspots can be

| Accessibility of suitable habitat through dispersal and migration
Land use, river distance and elevation contributed 45.5%, 36.4% and 18.2%, respectively, to the resistance matrix, and effective resistances had a high explanatory power for genetic distance (MLPE: Response ∼ RCPScenario × Mitigation

F I G U R E 3 Predicted suitability of
Primula elatior in Europe. The density curve displays the cumulative suitability as a function of latitude (y) and longitude (x). The predicted suitability ranges from 0 to 1 and is displayed in the colour key bar R 2 marginal = .76, R 2 conditional = .92, and the partial Mantel statistic was 0.77). The predicted accessibility of habitat patches decreased from 94.1% (current projected distribution) to 72.6% for RCP 2.6, to 72.8% for RCP 4.5 and to 70.3% for RCP 8.5 (SE < 0.1 in all scenarios). The low migration potential resulted in a reduction of the accessible distribution area (RO II) in all scenarios, where 33,510 km 2 (−40.2% compared to the current scenario) of the projected distribution area was accessible for RCP 2.6 compared to 28,274 km 2 (−49.5%) for RCP 4.5 and 23,735 km 2 (−57.6%) for RCP 8.5 (Figure 4; Figure 5). The limited migration potential thus resulted in 15.6 ± 1.7% (mean ± SD of climate change scenarios) of the projected total distribution in 2050 that was not accessible.

| Expected shifts in the habitat distribution and configuration
Projected habitat loss was most severe in continental France, and habitat gain was most prevalent in the temperate maritime climate zones of France, England and Norway (Figure 4) [198.3, 203.4] m (chi-square = 2,638.6, p = <.001, df = 6). However, the reduction in median elevation was most likely related to the shift to a more maritime climate and was therefore not considered further (Appendix S8).
The distribution shift of habitat patches resulted in a reduced patch accessibility in the north and in the more temperate maritime regions (RO III; Table 2; Figure 6). The predicted average PRI of accessible habitat patches decreased from 210.4 in the current projected distribution to 158.3 for RCP 2.6, 116.1 for RCP 4.5 and 118.2 for RCP 8.5, with a standard error smaller than 0.1 for each scenario.
Furthermore, the PRI strongly decreased in the south under climate change (see Table 2 for test statistics and Figure 6 for predicted means and their confidence interval).
Overall, the median latitude of habitat patches is thus expected to shift 183.2 ± 34.8 km (mean ± SD of climate change scenarios) F I G U R E 4 Habitat loss (red), gain (green) and habitat relicts (unchanged; blue) of Primula elatior in Europe for RCP 2.6 and RCP 8.5, with (right) and without (left) forest restoration mitigation. Dispersal and migration limitations are taken into account for the displayed habitat change. All geographical calculations are based on the ETRS89extended/LAEA Europe and the scale is visualized in 10 6 metres on the left (latitude) and on the bottom (longitude). WGS84 graticules are indicated in degrees (curved projection) on the right (latitude) and on the top (longitude). The embedded figures display the proportional turnover (%) of habitat in loss, gain and relict (unchanged habitat), respectively (y-axis) in relation to the latitude (x-axis) in 10 6 metres (LAEA Europe) northwards and 58.1 ± 9.3 km (mean ± SD) to more maritime climate regions. Furthermore, the PRI is expected to reduce by 37.8 ± 11.3% (mean ± SD) under climate change.

| Mitigation effects of large-scale afforestation
Reforestation mitigated 115.6% (net gain) of the total distribution area loss for RCP 2.6, 73.7% for RCP 4.5 and 55.9% for RCP 8.5. The mitigation of the accessible distribution area loss amounted to 73.9% for RCP 2.6, 50.3% for RCP 4.5 and 38.2% for RCP 8.5 (RO IV). The median northward latitudinal shift was only slightly mitigated for RCP 2.6 (1.3% reduction; p = .037) and not significantly mitigated for RCP 4.5 (5.9%, p = .108). However, 6.9% of the median northward shift was mitigated for RCP 8.5 (p = <.001). The median predicted shift from continental to more maritime regions was considerably mitigated, with a 35.6% reduction for RCP 2.6 (p = <.001), 34.8% reduction for RCP 4.5 (p = <.001) and 30.0% reduction for RCP 8.5 (p = <.001; see Appendix S8 for post hoc result tables and an overview of scenario medians and their 95% confidence intervals).
Overall, reforestation mitigated 81.7 ± 30.7% (mean ± SD of climate change scenarios) of the total distribution area loss, 54.1 ± 18.2% (mean ± SD) of the accessible habitat area loss and 33.5 ± 3% (mean ± SD) of the shift in the median habitat patch continentality. More habitat became available in the north and more maritime regions under reforestation management ( Figure 5 and Appendix S8). Therefore, the predicted mean (±SD) patch accessibility decreased in 2050 to 65 ± 0.9% compared to 72 ± 1.5% without reforestation management (Figure 7). However, reforestation mitigated 49.5 ± 4.2% (mean ± SD) of the loss in PRI in the accessible distribution area due to climate change (see Table 2 for test statistics and Figure 7 for predicted confidence intervals).

| Suitability drivers
The high contribution of landscape-scale variables in our species distribution model, with land use, distance to rivers and elevation as main predictors, points at the importance of ecosystem characteristics and landscape-scale processes for the occurrence of P. elatior. The importance of distance to rivers confirms its known ecological niche of rich mesic soils with loam deposits and high organic matter content, typical of floodplain forests (Taylor & Woodell, 2008). Furthermore, the importance of local and landscape-scale conditions for the occurrence of P. elatior is in line with other studies examining such drivers on forest herbs (Bernhardt-Römermann et al., 2015;Greiser et al., 2020;Valdés et al., 2015). Even though local factors are important drivers, climate still has a considerable contribution to the species' distribution (21.7%).
A clear optimum in temperature seasonality implies a dependence on annual variation in temperatures, but also a sensitivity to temperature extremes. This coincides with earlier findings that exposed a strong phenological sensitivity in P. elatior . The positive relation between the precipitation of the warmest quarter and P. elatior occurrences confirms the intolerance to desiccation of the species (Whale, 1983). Climate change effects are projected to be more severe in southern continental regions (Fick & Hijmans, 2017), and therefore, it is likely that conditions would shift outside of the species' climatic niche optimum in these regions. In the more maritime and northern regions, on the other hand, new available habitat is likely to become available because the climatic niche would shift into the species' optimum ( Figure 4).
Expected habitat loss in the continental south exceeds habitat gain in northern and maritime regions and the total distribution area is likely to reduce considerably by 2050 (46.4 ± 13.9%; mean ± SD). This coincides with the coarse median suitability loss of forest herb flora distributions in Europe predicted by Skov and Svenning (2004). However, one limitation of our study is that we did not take the climatic buffering effect of forests directly into account (Lembrechts et al., 2018;Zellweger et al., 2020). Furthermore, species responses to macro-climatic changes might be attenuated by plant phenotypic responses and/or local adaptation of fitness-related traits (Benito Garzón et al., 2019;Razgour et al., 2019). Integrating both temperature measurements of forest and riparian microclimates, and the adaptive potential of populations in predictions could further improve accuracy in the future.

F I G U R E 5
Scenario-specific summaries of the projected total distribution area (TA), accessible area (AA), within-patch dispersallimited area (WLA), migration-limited area (MLA) and dispersallimited area (DLA) of Primula elatior in Europe

| Predicting spatio-temporal dispersal patterns
Due to the low predicted migration potential of P. elatior to new available habitat in 2050 ( Figure 5), it is estimated that 15.6 ± 1.7% (mean ± SD) of the projected total distribution area will be inaccessible for colonization. The low migration potential of P. elatior is reflected by its inherent life-history traits, such as the long period until first reproduction (~3 years), specific germination requirements and absence of specific dispersal mechanisms (Taylor & Woodell, 2008;Verheyen et al., 2003). Honnay et al., (2002) and Jacquemyn et al., (2002) both empirically determined the dispersal potential of P. elatior in fragmented landscapes and found that colonization after 40-50 years was virtually zero when habitat patches were more than 1,000 m apart. This roughly coincides with modelled dispersal rates of 24 to 95 m per year in forest species without adaptations to animal dispersal vectors . Correspondingly, the mean migra-  (Engler et al., 2012), but this process would further complicate the genetic optimization process. In addition, efforts that link field-verified colonisation rates over time (e.g. individual-based spatially explicit modelling ;Landguth et al., 2010) with genetic divergence could improve the temporal estimation of dispersal models. Taking these limitations in account, our study integrates the genetic and landscape structure on a continental scale which results in highly accurate isolation-by-resistance predictions on a habitat patch scale. The functional dispersal patterns of P. elatior emphasize the importance of streams as efficient corridors TA B L E 2 Climate change effects on the accessibility and proximity resistance index (PRI) of projected Primula elatior habitat patches in Europe. "Between scenario" climate change effects are determined by the projected scenario (factor), geographic effects are determined by latitude and continentality (continuous) and the distribution shifts (within scenario) are determined by the interaction of geographic effects with scenario's. N depicts the total amount of projected habitat patches in each scenario and values correspond to z-statistics (Accessibility) and t-statistics (PRI) of the fixed terms. Effects on patch accessibility (binary) are based on the total distribution area and the proximity resistance index (PRI) is based on the accessible area. The intercept is the current scenario and the minimum latitude is corrected to the most southern patch (latitude-min{latitude}). Residual variance is largely related to regional differentiation in observed patterns (Appendix S8). The reforestation section depicts the overall mitigation effect and the scenario-specific mitigation effects are determined by the interaction between climate change scenario and a reforestation factor (binary). The intercept depicts the corresponding climate change scenario without reforestation. Each row and each column in the reforestation section corresponds to a distinct model and the intercept corresponds to the overall mitigation model. Scenario-specific mitigation effects are calculated on data subsets (future scenario with and without reforestation) and Chrysosplenium oppositifolium (Verheyen et al., 2003). Validation of our approach for these and other forest herbs could further optimize management decisions to protect and restore forest understorey communities in a context of global change.

| The impact of climate change on habitat distribution and configuration
Climate change will likely result in a northward shift of the median patch distribution (183.2 ± 34.8 km; mean ± SD) and to more maritime regions (58.1 ± 9.3 km; mean ± SD). The low migration potential to new available habitat is estimated to result in a general reduction of patch accessibility in these regions ( Figure 6 and Appendix S8). This indicates that P. elatior will not be able to track the changing climate, and decrease in dispersal probabilities between habitat patches (Appendix S8) will likely cause a general reduction of the metapopulation stability (as determined by PRI) in these regions (Table 2; Figure 6).
P. elatior is known to quickly lose genetic diversity when the metapopulation stability decreases (Jacquemyn et al., 2009), and therefore, southern remnant populations are estimated to still be threatened in the long term.
Interestingly, populations that were introduced in Norway during the late 19th century are currently in the naturalisation process (Gederaas et al., 2012). The ten regions currently occupied there, are in the close vicinity of habitat gain in the future, and could drive local migration patterns. Therefore, current naturalizing populations could become part of the shifting range. Overall, this illustrates the value of our modelling strategy for guiding assisted migration scenarios (Hunter-Ayad et al., 2020). Combining SDMs with dispersal models thus have a strong potential for guiding conservation and restoration efforts, but field validation remains critical (Laliberté & St-Laurent, 2020).

F I G U R E 6
Interaction effects between climate change (scenario) and the geographical position of Primula elatior habitat patches (as determined by latitude and continentality) on patch accessibility and the proximity resistance index. The predicted means and their 95% confidence interval (colour fill) relate to the within-scenario effects in Table 2. The latitude is represented in meters based on the LAEA89 coordinate system and continentality is represented in km distance from the coast. The patch accessibility is based on the total distribution area and the proximity resistance index is based on the total accessible area. The northernmost known native populations (green) can be found around 3.7 × 10 6 m latitude (LAEA89) or 57° North (WGS84

| Mitigating climate change effect of forest herbs
Reforestation will likely mitigate most of the projected habitat loss due to climate change (81.7 ± 30.7%; mean ± SD) and a substantial amount of the accessible habitat area loss (54.1 ± 18.2%; mean ± SD). The mitigation of the shift in the median patch latitude is limited because mitigating effects in the south of the range are accompanied by habitat gain in the north of the range ( Figure 5).
Therefore, patch accessibility is not positively affected by reforestation. However, habitat gain in continental regions in the centre of the distribution ( Figure 5) results in a considerable mitigation of the median patch continentality shift. These areas are generally in the vicinity of existing populations and are therefore accessible for migration. Habitat gain in the vicinity of these remaining populations results in a considerable mitigation of the loss in metapopulation stability (49.5 ± 4.2%; mean ± SD), as determined by the PRI. However, it is important to note that colonization success of reforested and afforested habitat patches is strongly dependent on past land use patterns. Successful colonization is often delayed by decades to centuries and short-term mitigation effects could therefore turn out to be lower than anticipated (Baeten et al., 2010;Naaf & Kolk, 2015).

| CON CLUS ION
The sensitivity of Primula elatior to the change in seasonal climate patterns and rainfall during the warmest quarter will likely result in a considerable loss of the total distribution area by 2050. Furthermore, due to its limited migration capacity, it is estimated that the species will be unable to track the shifting climate. The metapopulation stability will likely decrease across the range and southern populations will be most affected. Reforestation is estimated to considerably mitigate the loss of the distribution area, and metapopulations are predicted to remain more stable. However, to completely alleviate habitat loss it seems required to integrate strategies for climate change mitigation (RCP 2.6), reforestation, restoration of ecological connectivity and assisted migration.

This research was funded by the Flemish Research Foundation
(FWO project G091419N). Genetic data were collected in accordance with the Nagoya protocol on access to genetic resources (NOR: TREL1902817S/128). We would like to thank local biodiversity agencies for providing occurrence data. Furthermore, we would like to thank Koenraad Van Meerbeek for useful comments.

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

PE E R 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.13367.

DATA AVA I L A B I L I T Y S TAT E M E N T
Training and validation data for the species distribution modelling, projections across space and time, the isolation by resistance model and corresponding data, and shapefiles containing patch-specific characteristics are deposited in a Dryad data repository: https://doi. org/10.5061/dryad.msbcc 2fxx.

R E FE R E N C E S
Aiello-Lammens, M. E., Boria, R. A., Radosavljevic, A., Vilela, B., & Anderson, R. P. (2015). spThin: An R package for spatial thinning F I G U R E 7 Predicted 95% confidence intervals of habitat patch accessibility and metapopulation stability, as determined by the proximity resistance index (PRI). The Y-axis depicts the projected current distribution ("Current") of Primula elatior, the projected distribution in 2050 under climate change ("2050") and the projected distribution in 2050 under climate change with large-scale reforestation/afforestation ("2050 mitigation"). The predicted mean and confidence intervals relate to the GLM models used to assess mitigation effects of reforestation with RCP Scenario × Mitigation as explanatory variable