# A novel statistical workflow using pollen records and regression kriging to reconstruct the spatially and temporally explicit demographic history of tree species

Azzurra Pistone, Denis Allard, Christoph Schwörer, César Morales-Molino, Willy Tinner, Katalin Csilléry

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Azzurra Pistone, Denis Allard, Christoph Schwörer, César Morales-Molino, Willy Tinner, et al. A novel statistical workflow using pollen records and regression kriging to reconstruct the spatially and temporally explicit demographic history of tree species. Ecological Informatics, 2026, 93, pp.103561. 10.1016/j.ecoinf.2025.103561. hal-05482213

HAL Id: hal-05482213

https://hal.inrae.fr/hal-05482213v1

Submitted on 28 Jan 2026

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_3.png" /></figure>

# A novel statistical workflow using pollen records and regression kriging to reconstruct the spatially and temporally explicit demographic history of tree species

Azzurra Pistone<sup>a,c ID*</sup>, Denis Allard<sup>b ID</sup>, Christoph Schwörer<sup>c ID</sup>, César Morales-Molino<sup>c,d ID</sup>, Willy Tinner<sup>c ID</sup>, Katalin Csilléry<sup>a ID</sup>

## ARTICLE INFO

Dataset link: https://github.com/azzpis/pollen kriging

**Keywords:**

Climate change

Geostatistics

Holocene

Niche

Range expansion

Species abundances

## ABSTRACT

Understanding the effects of past climate shifts on the demography of forest tree species is crucial to assessing their response to ongoing climate change. However, little effort has so far been made to quantify past demographic changes in a spatially and temporally explicit manner at a continental scale. We have developed a novel statistical workflow that integrates two regression kriging models to reconstruct the demographic history of tree species across Europe. Our workflow anticipates spatially the probability of species occurrence (PoO), and interpolates their relative abundances (RelAb) spatially and temporally. Climate variables can be included as covariates. Our approach can accommodate non-stationary species responses to climate, and incorporates the presence of source populations, colonization constraints, and population trends as factors influencing species RelAb. We applied this workflow to European fir species (Abies spp.) since the Last Glacial Maximum (LGM), using fossil pollen records from 241 sites, and simulated paleoclimate data on a 0.41-degree grid and 500-year time bins. Model performance, assessed with cross-validation, demonstrates that including climate as a covariate enhances the spatial heterogeneity. Climate has a positive effect on RelAb interpolation under millennial static spatial distribution structure conditions, while the presence of source populations plays a more important role during rapid demographic processes. Additionally, we applied our workflow to assess future regional changes in the RelAb of Abies spp. under the main future climate scenarios. Our workflow is particularly suited for temperate and boreal tree species and can be used in various downstream analyses.

## 1. Introduction

Many forest ecosystems worldwide are increasingly suffering from the impact of global change, evidenced by large-scale die-back events (Allen, 2009; Hammond et al., 2022; Hartmann et al., 2022), reduced growth, and rising incidence of diseases (Lindner et al., 2010; Sturrock et al., 2011; Bennett et al., 2015). Additionally, tree responses to ongoing climate change may worsen in the future due to adaptation lags associated with their long generation times (Davis et al., 2005; Kremer et al., 2012; Aitken and Bemmels, 2015). While ongoing climate and land-use changes are reaching an unprecedented speed (Masson-Delmotte et al., 2021), the paleoecological record has revealed the past occurrence of dramatic changes in forest vegetation caused by

abrupt climate shifts like the rapid warming events (ca. 3–4 °C in a few decades) that occurred at the onset of the Bølling-Allerød Interstadial and the Holocene in the Northern Hemisphere (Birks and Ammann, 2000; Magri, 2007; Ammann et al., 2013; Ruosch et al., 2016; Lang et al., 2023). Understanding how such past rapid climate shifts triggered demographic changes in forest tree species may help assess their future responses to climate change. However, despite its interest, little effort has so far been made to quantify past demographic changes in forest tree populations from a spatial and temporal perspective. Demographic reconstructions that explicitly quantify changes in population size at a continental scale may improve predictive models and their application in forest management (Pecchi et al., 2019), particularly if they

<footnote><sup>a</sup> Evolutionary Genetics Group, Biodiversity and Conservation Biology, Swiss Federal Research Institute WSL, Zürcherstrasse 111, 8903 Birmensdorf, Switzerland</footnote>

<footnote><sup>b</sup> Biostatistics and Spatial Processes (BioSP), INRAE, 84914 Avignon, France</footnote>

<footnote><sup>c</sup> Institute of Plant Sciences and Oeschger Centre for Climate Change Research, University of Bern, Bern, Switzerland</footnote>

<footnote><sup>d</sup> Grupo de Ecología y Restauración Forestal (FORECO), Departamento de Ciencias de la Vida, Facultad de Ciencias, Universidad de Alcalá, Alcalá de Henares, Spain</footnote>

<footnote>* Correspondence to: Evolutionary Genetics Group, Biodiversity and Conservation Biology Unit, Swiss Federal Research Institute WSL, Zürcherstrasse 111, 8903 Birmensdorf, Switzerland.

E-mail address: azzurra.pistone@wsl.ch (A. Pistone).</footnote>
incorporate temporal changes in species niches driven by competition, adaptation, land use, and colonization lags (Benito Garzón et al., 2007; Elith et al., 2010).

For over a century, fossil pollen records have provided an invaluable wealth of information on past vegetation dynamics (von Post, 1916; Erdtman, 1969; Prentice, 1988). More recently, the increasing abundance of fossil pollen records has led to the development of large databases of chronologically and taxonomically curated and harmonized stratigraphical pollen datasets (Fyfe et al., 2009; Giesecke et al., 2013; Williams et al., 2018; Lang et al., 2023). However, translating fossil pollen assemblages into quantitative vegetation reconstructions is challenging and only a few models have addressed it. First, REVEALS (Sugita, 2007a,b) produces pollen-based reconstructions across various spatial scales based on site type and size, pollen counts, estimates of pollen productivity and their related error, estimates of pollen fall speed, and other parameters involved in the dispersal and deposition functions (Sugita, 2007a,b). Although REVEALS has been applied to reconstruct land cover as well as the distribution and abundance of some species (Githumbi et al., 2022), its applicability is often limited by the lack of data on several of the parameters related to pollen production and transport, particularly relative pollen productivity data (Liu et al., 2022; Serge et al., 2023). Further, the spatial extent of such reconstructions is often local, with exception of subcontinental, continental or northern hemispheric overviews, which however are rarely continuous in space (Githumbi et al., 2022; Schild et al., 2025). This means that REVEALS remains site-specific and cannot provide continuous reconstructions across a continent. Second, STEPPS (Paciorek and McLachlan, 2009; Dawson et al., 2016) is a hierarchical Bayesian model that uses pollen counts from a network of pollen records, along with forest inventory data, to reconstruct spatio-temporal variations in forest composition or the relative abundance of a taxon. This model considers trends in space and time and estimates prediction uncertainties. However, its calibration requires vegetation surveys, which are often unavailable (Trachsel et al., 2020). Finally, the ROPES approach (Theuerkauf and Couwenberg, 2018) is based on pollen influx or accumulation rates (i.e., pollen grains $cm^{-2} year^{-1}$) and accounts for changes in pollen productivity estimates over time. This makes ROPES suitable for unbiased reconstruction at a specific site (Theuerkauf and Couwenberg, 2018). However, this model has never been tested at a continental scale. Indeed, harmonizing pollen accumulation rates across a large dataset is challenging, requiring robust and precise chronologies and, ideally, annually laminated sediments, which are available from few sites. Finally, none of the models mentioned above includes environmental parameters and/or biological processes.

In this study, we developed and tested a new workflow based on geostatistical models that account for spatial and temporal autocorrelation of pollen presence and relative abundances to reconstruct past population size fluctuations. We used *Abies* in Europe as a case study, but the workflow is transferable to any other taxon, provided adjustments are made for its pollen dispersal strategy, and to other continents, given sufficient pollen data density. We focused on the widespread *Abies alba* Mill. (silver fir) and two more narrowly-distributed Mediterranean species, *A. cephalonica* Loudon and their hybrid *A. borisii-regis* Mattf. These species cannot be distinguished from pollen morphology; thus, throughout this paper, we referred to them simply as *Abies* spp. We chose this genus because its resilience to drought and bark beetle attacks makes it one of the best-suited European conifer genera to cope with ongoing global change in the future (Wermelinger, 2004; Lebourgeois et al., 2013; Vitali et al., 2017). Further, the postglacial dynamics of *Abies* spp. in Europe, particularly of *Abies alba*, remain unclear, and our study can help fill this knowledge gap (Liepelt et al., 2009; Ruosch et al., 2016).

In our workflow, we integrated fossil pollen records with climate data to “fill the gaps” in areas and time periods without pollen records. To achieve this, we implemented a geostatistical interpolation approach

using regression kriging (Hengl et al., 2004), a well-established method in environmental science that has been shown to outperform non-geostatistical techniques such as inverse distance weighting (IDW) (Li and Heap, 2011). We applied this method separately to presence/absence and abundance (relative percentage) data to address the high proportion of zero values in the pollen data. Similar spatial kriging of presence/absence data has recently been used to produce species distribution maps from modern pollen data (Oriani et al., 2024). However, applying spatio-temporal regression kriging to abundance data is novel. This geostatistical technique is suited for variables that vary both in space and time, such as weather (Hu et al., 2017; Sekulić et al., 2019) or pollutant concentrations (van Zoest et al., 2019). Specifically, we applied spatio-temporal kriging to interpolate the relative pollen percentages of *Abies* spp., accounting for population trends over time and space. Our workflow has a relatively low computational cost, making it feasible to include additional explanatory variables, such as climate data, competition, and human-impact indices, including the land-use probability index (LUP) (Deza-Araujo et al., 2022). Additionally, the statistical nature of our method enables rigorous cross-validation, quantification of total error, and monitoring through hindcasting. As a final aspect of this study, we applied our workflow to predict the relative abundances of *Abies* spp. under different future climate trajectories.

## 2. Material and methods

### 2.1. Workflow

Given the zero-inflated characteristic of the dataset, we divided the workflow into two components: (1) a binary model that interpolates the probability of occurrence (PoO) and (2) a continuous component that models the positive non-zero values, i.e., the relative abundances (RelAb) (Fig. 1). Both components combine regression with kriging, leveraging auxiliary explanatory variables and spatial (binary component) or spatio-temporal (continuous component) autocorrelation to make more accurate estimates (Fig. 1). We excluded the missing data from the interpolation to avoid the introduction of false zeros and the imputation of the missing values.

First, we modeled the relationship between the variable of interest (absence/presence or RelAb) and the explanatory variables using a regression model: a logistic regression model for the binary component and a generalized linear model with a logit link for the positive values. This step captures the trends driven by known covariates. For the binary data, we modeled each year separately, resulting in 40 regression models. For RelAb, we first applied a logit transformation to approximate a normal distribution, as the data were highly right-skewed. Then, we divided the dataset into four main eras: the Last Glacial Maximum (LGM) (20–19k CalYrBP, Calibrated Year Before Present, with present defined as 1950 CE) and the Late Glacial (19–11.5k CalYrBP; 11.7k CalYrBP rounded to 11.5k CalYrBP throughout the manuscript for time binning reasons), the Early Holocene (11.5–8.5 CalYrBP), the Middle Holocene (8.5–4.5k CalYrBP), and the Late Holocene (4.5–0k CalYrBP). We first fitted separate regression models for each era to accommodate the limited data available for the LGM and the Late Glacial and to reflect the species’ known history. We then extracted the residuals from each model and aggregated them into a single dataset. In the regressions, we used the PoO of the previous time step (derived from the binary component), climate variable sets, and their combinations as explanatory variables (Tables 1 and 2).

After the model fitting, we extracted the residuals, which were assumed to contain spatially (and temporally) structured variation not explained by the covariates. These residuals were then interpolated using kriging (spatial kriging for the binary component and spatio-temporal kriging for the continuous component). This step required the construction and fitting of a variogram, which captures the spatial (or temporal) variability of the observed data as a function of their distance (Oriani et al., 2024). For the binary component, we
<figure><img src="image_4.png" /><figcaption>Fig. 1. Graphical abstract of the workflow.</figcaption></figure>

<table><caption>**Table 1**
Sets of climatic variables used in the binary and continuous components.</caption><thead><tr><th>Variables set name</th><th>Abbreviation</th><th>Used in</th><th>Variables included</th></tr></thead><tbody><tr><td>Variables good for the establishment and survival of trees (Mauri et al., 2022)</td><td>VarSel_Surv</td><td>Binary and continuous components</td><td>Summer temperature, winter and summer precipitation sums, precipitation seasonality (bio15), continentality</td></tr><tr><td>Expert choice for Abies alba (e.g. Bledý et al., 2024)</td><td>VarSel_Exp</td><td>Binary and continuous components</td><td>Mean temperature of the coldest month (bio06), mean temperature of the warmest quarter (bio10), precipitation of the warmest quarter (bio18)</td></tr><tr><td>Variables selected for the binary component through random forest and regression methods</td><td>VarSel_IK</td><td>Binary component</td><td>Total summer precipitation, minimum monthly winter temperature, mean summer temperature</td></tr><tr><td>Variables selected for the continuous component through random forest and regression methods</td><td>VarSel_STRK</td><td>Continuous component</td><td>Temperature of the driest quarter (bio09), mean summer temperature, mean winter temperature, isothermality (bio03)</td></tr></tbody></table>

<table><caption>**Table 2**
Table of all models used in this manuscript. Variables defined in Table 1.</caption><thead><tr><th>Probability of occurrence (PoO)</th><th>Relative abundance (RelAb)</th></tr></thead><tbody><tr><td>PoO ~ 1</td><td>RelAb ~ PoO<br>RelAb ~ PoO + VarSel_Exp<br>RelAb ~ PoO + VarSel_Surv<br>RelAb ~ PoO + VarSel_STRK</td></tr><tr><td>PoO(VarSel_Exp) ~ 1 + VarSel_Exp</td><td>RelAb ~ PoO(VarSel_Exp)<br>RelAb ~ PoO(VarSel_Exp) + VarSel_Exp</td></tr><tr><td>PoO(VarSel_Surv) ~ 1 + VarSel_Surv</td><td>RelAb ~ PoO(VarSel_Surv)<br>RelAb ~ PoO(VarSel_Surv) + VarSel_Surv</td></tr><tr><td>PoO(VarSel_IK) ~ 1 + VarSel_IK</td><td>RelAb ~ PoO(VarSel_IK)<br>RelAb ~ PoO(VarSel_IK) + VarSel_STRK</td></tr><tr><td>/</td><td>RelAb ~ VarSel_Exp</td></tr><tr><td>/</td><td>RelAb ~ VarSel_Surv</td></tr><tr><td>/</td><td>RelAb ~ VarSel_STRK</td></tr></tbody></table>

constructed and fitted a variogram for each time bin independently, using the exponential correlation model, as in previous studies (e.g. Chiles and Delfiner, 2012; Dubrule, 2017; Oriani et al., 2024). For RelAb, which is spatially and temporally explicit, we chose the most appropriate correlation model based on the non-separability index. This index indicates how the spatial and temporal components of a covariance function interact (Allard et al., 2022). We found high non-separability for short distances and full separability beyond 45 km, indicating that spatial correlation depended on the temporal lag only under local conditions (Fig. S1). We thus selected the separable Matern spatio-temporal correlation model (Bevilacqua et al., 2024). Additionally, we employed a spatio-temporal dynamic function that accounts for changes in spatial locations over time to address the missing values in the pollen records. The missing values are thus interpolated using the same method as for other grid cells, but their estimated value benefits from the additional knowledge provided by relative abundances at other time points. We computed both the variogram and the kriging of the RelAb residuals with a limited spatial correlation range of 500 km,
while for the PoO we applied a limit of 1000 km. Importantly, since the temporal interpolation is symmetrical, for a given time bin, both “relative past” and “relative future” data are used for the interpolation in our approach, except for the last time bin, which is based on past data only.

Then, we predicted values at new locations by summing the stochastic (kriged residuals) and deterministic (regression model's fitted values) components of the workflow. To obtain the regression fitted values at new locations in the binary component, we used the same 40 regression models from which the residuals were previously extracted. We then added the interpolated residuals to the fitted values and backtransformed the resulting values using inverse logit, resulting in 40 maps that represented the PoO. For RelAb, we used a single regression model for the LGM and Late Glacial (20–11.5k CalYrBP) and a single model using the first five Holocene time bins (11.5–9.5k CalYrBP, 2500 years). Then we predicted the unknown values at each time bin independently by using the corresponding regression model. Finally, we added the residuals of the corresponding time bin. We aggregated these Early Holocene time bins because populations were still restricted mainly to refugia and data were limited, making it impractical to develop separate models for each time bin. From 9k CalYrBP onward, we trained a separate regression model using data from a six time-bins window and fitted it to the youngest time bin; we then added the interpolated residuals of the corresponding time bin to the prediction. The values for each resulting map were then backtransformed using the inverse logit.

Finally, we cropped the resulting RelAb maps using the PoO surfaces of the corresponding time bin derived from the binary step. We assigned grid cells with a PoO <0.2 as “predicted absences”. This threshold was chosen to minimize the prediction of false occurrences and false absences (see Supplementary Results).

For this study, we predicted the RelAb under three future climate scenarios: ssp126, ssp370, and ssp585 (Karger et al., 2017). To achieve this, we included an additional time step in the kriging interpolation, enabling the prediction of the residual surface based on previous time steps. Additionally, we used Late Holocene data to train the second regression model and fitted it at unknown locations using the predicted climate for the next 100 years under different climate scenarios. We did not crop the final maps.

We performed both kriging methods using the GeoModels package (Bevilacqua et al., 2024) in R 4.4.1 (R Core Team, 2024).

## 2.2. Pollen dataset

We used the selection of European paleoecological sites compiled by Lang et al. (2023), which comprises 323 representative, well-dated, and chronologically highly resolved post-glacial sites (20 000 CalYrBP to present) that cover the entire continent. We selected a subset of 241 sites in which our target species, Abies spp., could have been present in the past (Table S1). We excluded sites where Abies pollen most likely corresponded to Abies sibirica Ledeb., based on distributional criteria (latitude >60°N and longitude >40°E). We also excluded marine cores and those from the Balearic Islands. Finally, we excluded 16 sites due to the unavailability of the original datasets (Table S1).

Seven sites had consecutive, parallel cores that we checked for chronology and merged to achieve greater time coverage and reduce pollen detection uncertainty. The quality and quantity of pollen and climate data are key to our proposed purely statistical workflow. In our study, 172 sites had samples covering 20–11.5k CalYrBP, and 69 sites covered only the Holocene, with a median of one or two samples per time bin (defined below) from 20k to 15k CalYrBP and a median of three samples per time bin from 14.5k to 0 CalYrBP. Abies spp. occurred in 108 of the 241 sites, and the empirical limit (i.e., continuous pollen curves; Lang et al., 2023) was reached around 8.5k CalYrBP in most of them (Fig. 2A and B). Before the onset of the Holocene (11.5k CalYrBP), only a few sites recorded the presence of Abies spp. located in the

refugial zones in the Balkans (six sites), Apennines (four sites), and southern Alps (two sites) (Fig. 2A).

Different species may be over- or underrepresented in pollen assemblages due to their ecological characteristics. Thus, we applied correction factors to weigh pollen records according to species pollen productivity (Liu et al., 2022; Lang et al., 2023) (Table S2, values from Lang et al. (2023)). Additionally, we applied an empirical limit, following Lang et al. (2023). The empirical limit is defined as the presence of the relevant pollen taxon, in this case Abies, in at least five consecutive samples that cover at least 1000 years, with a maximum gap (absence) of one sample (Lang et al., 2023). It is interpreted as evidence of the establishment of local (initial) populations. Here, we set all observations preceding the empirical limit to zero. Then, we divided the dataset into 500-year time bins, resulting in 40 bins for 0–20k CalYrBP. The large size of the time bins was justified by the need to have as much pollen as possible in each bin (Chevalier et al., 2020; Githumbi et al., 2022; Lang et al., 2023) and to avoid the introduction of too many missing values, which could impact the regression component results by providing too few data points. Moreover, previous studies showed that a 500-year time bin is a suitable compromise between data availability and prediction quality (Finsinger et al., 2017; Giesecke et al., 2019; Mottl et al., 2021). For each time bin, we calculated the percentage of Abies pollen (y) relative to the sum of terrestrial tree, shrub, and herb pollen as:

$$y(s, t_{\Delta}) = \frac{\sum_{i=t_1}^{t_2} c_{Abies(s,t_i)} \times \alpha_{Abies}}{\sum_{i=t_1}^{t_2} \sum_{k=1}^{K} c_{k(s,t_i)} \alpha_k} \times 100$$

where y = Abies pollen %, s = site, t_Δ = t_1 - t_2 = time bin, α = correction factor for species k, k = species, K = species of trees, shrubs and herbs found at site s, and c = pollen count. If a bin did not contain any sample at a site (e.g., due to a lack of deposits or damage), it was classified as missing data. We refer to the relative percentages of Abies spp. as relative abundances (RelAb).

## 2.3. Explanatory variables

We tested a set of bioclimatic variables derived from the down-scaled climate model datasets CHELSA-TraCE21k and CHELSA v2.1 (Karger et al., 2017, 2020, 2021, 2023) as covariates in the regression kriging. The spatial and temporal resolutions of these simulated data were higher than those of the pollen data, so we averaged the values over 500-year time bins. Further, we reduced the spatial resolution to 0.41 degrees (~32 × 32 km) using the Terra package (Hijmans, 2020) in R. We chose this size to balance the aim of reconstructing regional presence with the diverse site types (e.g., lakes, bogs) and sizes (Githumbi et al., 2022; Lang et al., 2023). Additionally, Liu et al. (2022) found that for the Abies genus, the 75% cumulative influence is reached at 10–20 km from the depositional site, thus a gridcell size of circa ~32 × 32 km would capture the regional signal of pollen, climate, and crucially, their relationships.

Many of the bioclimatic variables were highly correlated with one another; we therefore selected several reduced sets of bioclimatic variables that best described the species' ecological niche (Table 1) (Svenning et al., 2011). First, we selected a set of three variables related to the establishment and survival of the species, based on the available literature on European Abies spp. ecology (e.g., Bledý et al., 2024): mean temperature of the coldest month (bio06), mean temperature of the warmest quarter (bio10), and precipitation of the warmest quarter (bio18). We refer to this set of climate variables as VarSel_Exp. Second, we used variables that are considered the main drivers of establishment and persistence of forest tree species (VarSel_Surv): winter and summer temperature and precipitation, precipitation seasonality, mean annual temperature (bio01), mean temperature of the coldest month, total yearly precipitation (bio12), continentality, a humidity index, growing degree days above 5 °C (GDD), soil pH, and soil organic carbon content
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 2.** Map of sites and distributions of input data. (A) The empirical limit over time and space and the present-day distribution of European *Abies* spp. Ka: thousand years. Distribution of *Abies* spp. from Caudullo et al. (2017) and isochrones from Lang et al. (2023). (B) Number of sites for each time bin. The light gray bars indicate the number of missing data points across all the datasets. (C) Boxplot of *Abies* spp. relative pollen percentages in each time bin. Graphs (B) and (C) were made after the empirical limit application.</figcaption></figure>

(OCC) (Mauri et al., 2022). We ultimately used only summer temperature, winter and summer precipitation, precipitation seasonality, and continentality, due to collinearity and the availability of paleosoil data (Table 1). Finally, we used a machine learning approach with four methods to infer the two sets of variables (VarSel_IK and VarSel_STRK) that best explained pollen presence and RelAb at our sites over the whole 20k-year period. We excluded monthly data from this process to reduce the data's dimensions and avoid overfitting. We calcu-

lated variable importance when possible and retained the variables selected by all four methods. We tested all sets of selected variables for collinearity using the variance inflation factor (VIF), condition index, and variance proportions using the car and mctest packages in R (Supplementary Results). We implemented the machine learning methods random forest, MARS, VSURF, and stepwise regression using R (Liaw and Wiener, 2002; Milborrow et al., 2007; Kuhn and Max, 2008; Kursa and Rudnicki, 2010; Genuer et al., 2013; Probst, 2018).
<figure><img src="image_8.png" /><figcaption>Fig. 3. *Abies* spp. pollen probability of occurrence (PoO) and relative abundance (RelAb) in selected scenarios and years. (A–E) PoO surfaces calculated without using covariates. (F–J) RelAb interpolated values obtained by using only PoO (no climate) as a covariate. (K–O) Interpolated RelAb values calculated using PoO (obtained by integrating climate as a covariate (VarSel_IK)) and climate (VarSel_STRK). All grid cells containing PoO values lower than 20% have been removed. Distribution of *Abies* spp. from Caudullo et al. (2017). High resolution version of this figure is available online.</figcaption></figure>

## 2.4. Cross-validation

We assessed the performance of the workflow’s binary and continuous components using several validation approaches. We evaluated the accuracy of the binary component by examining all fitted variograms (Fig. S2). We also conducted five-fold cross-validation by randomly removing 20% of the locations and calculated Brier’s score (Rufibach, 2010; Allard et al., 2012) to assess prediction error. For the continuous component, we examined both spatial and temporal variograms (Fig. S3, S4) to evaluate model accuracy. We also implemented hindcasting to test the model’s temporal predictive performance. Specifically, we used a rolling time window of six time bins to predict the subsequent time point (Fig. 4B). We also performed five-fold cross-validation on the continuous component. We assessed model performance using Spearman’s correlation and the relative median error (RME) between observed and predicted RelAb values, for both hindcasting and cross-validation (five-fold validation ×100 iterations; Fig. 4C, S5). Folds not including all the time bins were excluded (~7%).

# 3. Results

## 3.1. Using the workflow on a fossil pollen dataset

We used variograms to evaluate our workflow and the differences between climate variable sets. In the binary component, robust variogram exponential fitting, indicating strong spatial autocorrelation, was observed from 9k CalYrBP onward and in some previous time bins (e.g., 14k CalYrBP; Fig. S2). Conversely, adequate variogram fitting in earlier periods was hardly obtained, since the spatial component was weak and allowed the regression component to overly influence predictions, occasionally placing *Abies* presence in unlikely areas, such as at the edge of the Scandinavian ice sheet (Fig. 3K, S2). This prediction far from the known refugial areas is not unexpected, as these are clustered and thus not adequate for kriging interpolation. Additionally, refugia are often characterized by microclimates not captured at our spatial scale. Indeed, robust fitting was associated with a greater number of sites with *Abies* spp. recorded, reflecting its expansion outside the refugia, and with the emergence of a stronger spatial structure (Fig. 2A). VarSel_Surv demonstrated the most consistent exponential variogram among the climate variable sets tested. In contrast, the scenario without covariates tended towards a more linear distribution.

Spatially, the predicted distribution showed minimal differences between scenarios with and without climate variables. Still, those with climate variables tended to produce more heterogeneous PoO maps driven by short-range values in the variogram (Fig. S2) and by climate data. Climate often has a “diminishing” effect on PoO, increasing the heterogeneity of the range boundary and extending areas of low PoO. Additionally, the presence of climate variables led to non-zero PoO values in areas with sparse or no site coverage, including Turkey, North Africa, and the Caucasus. This may indicate potentially suitable areas not yet reached by *Abies* or areas where other *Abies* spp. not included in our dataset, e.g., *Abies clicicica* Carriere, *A. nordmanniana* (Steven) Spach, and *A. subsp. equii-trojani* (Asch. & Sint. ex Boiss.) Coode & Cullen, were likely present based on their current distributions (Fig. 3K–O).

Once we calculated the PoO, we interpolated the RelAb using a spatio-temporal framework. The spatio-temporal variogram showed close agreement between the model and the residuals for the spatial component (Fig. S3). The temporal component also showed a reasonable fit to the residuals, particularly for time bins close together (Fig. S4). However, beyond approximately the 7th–8th time bin (corresponding to 3500–4000 years), the temporal variance increased more rapidly than the fitted variogram model (Fig. S4). Although these are not optimal results, the effect of this deviation is limited since the resulting weights are higher for shorter time intervals. Thus, correlations between close time points have a greater impact on interpolation than those between far-apart time points, a property known as the “screening effect” in geostatistics (Chiles and Delfiner, 2012).

From a spatial perspective, we observed greater variability and uncertainty during the LGM and Late Glacial periods than in the Holocene (see also Section 3.2). Nevertheless, in some scenarios, e.g., VarSel_Surv (PoO climate included), the workflow predicted the presence and RelAb of *Abies* spp. during the LGM and Late Glacial in areas not covered by the sites but where the species’ presence is most likely, e.g., Northern Italy. Starting from the Holocene, results considering the three climate scenarios diverged from those without climate variables. Indeed, in the latter, the species’ RelAb was strictly linked to the site values and distribution. Conversely, when climate was included, species RelAb values were more heterogeneous across space and influenced RelAb predictions in uncovered regions (Fig. 3K–O). As an example, the Adriatic Balkan coastline, which was poorly represented in our dataset (Fig. 2A), was predicted to have high RelAb values when VarSel_Surv
<figure><img src="image_9.png" /><figcaption>**Fig. 4.** Summary statistics over time of the cross validation results. OD: Oldest Dryas; B-A: Bølling-Allerød; YD: Younger Dryas; EH: Early Holocene; MH: Middle Holocene; LH: Late Holocene. (A) Results of the five-fold cross-validation of the binary component. The results for presences and absences are divided. High Brier's scores indicate incorrect detection, while low scores indicate correct detection. (B) Hindcast results of the continuous component. (C) Results of the five-fold cross-validation of the continuous component. (D) Spatial differences in predicted relative abundance (RelAb) excluding versus including climate variables (RelAb ~ PoO (climate excluded) vs. RelAb ~ PoO(VarSel_IK)+VarSel_STRK). Negative values (blue) indicate that the model predicts lower abundances when climate is included. Positive values (red) indicate that the model predicts higher abundances when climate is included. Values are the true differences (%). Distribution of *Abies* spp. from Caudullo et al. (2017).</figcaption></figure>

was used. Another example is the prediction of the species' present-day RelAb and distribution (Fig. 4D). In this last time bin, when climate is included, *Abies* spp. RelAb is higher in mountain regions and at the leading edge of the distribution, and lower in areas not occupied by the species (e.g., the Pannonian Plain) (Fig. 4D).

## 3.2. Cross-validation

Both the number of sites and their spatial distribution affected prediction accuracy. The Brier's scores from the five-fold cross-validation revealed a clear difference in performance between the LGM/Late Glacial and Early, Middle, and Late Holocene. In particular, they indicated a larger prediction error in species presence in the older time bins (Fig. 4A), most likely due to the small number of sites with *Abies* spp. records during this period and by their spatial distribution (Fig. 2A-B). Additionally, the positive error decreased significantly at 9k CalYrBP (Fig. 4A), when the species expansion began and the number of sites registering *Abies* spp. reached ~50 (of the ~200 available). Absences were adequately predicted across all time periods (Fig. 4A). The different sets of climate variables did not improve the prediction (Fig. 4A), indicating a strong influence of the spatial structure on the prediction.

The cross-validation of the positive values indicated an unstable Spearman's correlation and relative median errors (RMEs) for the LGM and Late Glacial (Fig. 4C, S5). During this period, no model consistently

outperformed the others. During the Holocene, cross-validation results showed that RME values across the different scenarios were overall very similar (Fig. S5). While VarSel_Surv performed slightly better and VarSel_Exp slightly worse, the differences were minimal (Fig. S5), and they consistently over-predicted RelAb starting from 5k CalYrBP. The Spearman's correlation remained highly unstable until the beginning of the Holocene (Fig. 4C). In the first time bins of this period, scenarios that did not include climate in the PoO showed the highest correlation between predicted and observed RelAb. The Spearman's correlations of all models decreased markedly in the Middle Holocene, reaching values of around 0.3 (Fig. 4C). Spearman's correlation values then recovered in the Late Holocene, only to decrease in the final time bins (last 1500 years; Fig. 4C). Notably, models that included only PoO (obtained with or without including climate) had the lowest correlation during the second half of the Middle Holocene and in the Late Holocene, with values reaching close to zero in the last time bin (the last 500 years; Fig. 4C).

At local conditions, including climate in the prediction appeared to play a less significant role than PoO alone (Fig. 4B). Hindcasting RME results showed similar outcomes among the models, with VarSel_Surv and VerSel_STRK emerging as the best models in the Late Holocene (Fig. 4B). Regarding correlation, hindcasting reached significantly higher values than cross-validation results (~0.88 vs. ~0.50; Fig. S5). The difference in Spearman's correlation values between cross-validation and hindcasting is expected, due to the distinct validation methodologies employed. In cross-validation, sites are removed entirely from the
dataset during validation, eliminating any demographic effects from other time bins. Conversely, in hindcasting, a site being predicted for a given time bin still benefits from the influence of its past data (i.e., the “legacy effect”). This suggests that, at the local level, the site trend plays a crucial role, whereas it may be less significant when predicting new sites ex novo. Similar to cross-validation, hindcasting correlation decreased, and RME increased in the last two time bins.

### 3.3. Post-glacial colonization history of *Abies* spp. in Europe

Based on regression models and random forest models across the past 20k years, analysis of the 29 simulated climate variables indicated that total summer precipitation, minimum monthly winter temperature, mean summer temperature, and temperature seasonality (bio04) were the most important predictors of *Abies* spp. presence, with minimum monthly winter temperature and total summer precipitation ranking as most important (Supplementary Results). The minimum monthly winter temperature and temperature seasonality showed high collinearity; therefore, temperature seasonality was removed from the analysis. For RelAb, a different set of variables was selected, including the mean temperature of the driest quarter (bio09), mean summer and winter temperatures, and isothermality (bio03) (Supplementary Results). The selected variables showed low collinearity (VIF < 4) (Supplementary Results). Although the variable selection relied on the entire dataset due to limitations at specific time points and was certainly biased towards time points when the species was most abundant (Fig. 2C), these data-driven sets complemented our published expert-based climate variable set.

We explored the demographic history of *Abies* spp. using the VarSel_STRK (PoO obtained by including climate) scenario (Fig. 3K–O). This set of climate variables resulted in intermediate RMEs and high Spearman’s correlations over time (Fig. 4C, S5, S6). During the LGM and the Late Glacial, *Abies* spp. were present in refugia on the Balkan Peninsula, from Greece to the southern Dinaric Alps, and in southern Italy, rarely exceeding 15% in RelAb. Unfortunately, no site in the Northern Apennines and the Pyrenees included in our dataset covered this period; however, the interpolated values suggest the possible presence of the species in the Northern Apennines, Western Alps, and Spain, as well as in other regions such as Turkey and the Italian coast (Fig. S6). The Northern Italian inference is supported by pollen records of *Abies* reaching the absolute limit (i.e., scattered presence) during the LGM in Northern Italy (e.g., Lago della Costa; Kaltenrieder et al., 2009).

While pollen data indicate that both the presence and the RelAb of the species started to increase between the Oldest Dryas (17.5k CalYrBP) and the Bølling-Allerød (15–12.9k CalYrBP), climate variables could not fill the gaps due to missing data (Fig. 3K). In the Younger Dryas (12.9–11.5k CalYrBP), populations’ RelAb were low and stable, with no further expansion until the Holocene (11.5k CalYrBP); conversely, some areas where the species had been present were no longer occupied, such as parts of the Northern Apennines and North Macedonia (Fig. S6). During the Younger Dryas, several climatically suitable areas were present (Fig. 5A), but the lack of nearby and/or sufficiently large source populations may have impeded their colonization (Fig. 2A and C). From the onset of the Holocene, *Abies* spp. started to spread, filling its climate niche by growing populations at locations where they already existed (Fig. 5B). For example, the Northern Apennine and Balkan populations reached RelAb values of ~45% and ~25%, respectively, in the Early Holocene (Fig. 3L–M). Further, the species started to reach cooler areas previously uncolonized (Figs. 3L–M, 5B), such as the Central Alps, where the species’ RelAb reached ~10% at 9k CalYrBP. Moreover, the colonization routes diversified, connecting populations from different glacial refugia (Fig. 3M). During the Middle Holocene, the species’ RelAb continued to increase, especially in wetter and cooler areas and all previously colonized regions (Fig. 5C). In contrast, in warmer and drier regions that contained predominantly

ancient populations, RelAb remained at low levels (~15%; Figs. 3M–N, 5C), occasionally reaching 35% in the Southern Balkans.

Starting from the second half of the Middle Holocene and Late Holocene, *Abies* spp. RelAb declined first in the Italian and Balkan peninsulas, then in the northern part of its distribution range. According to interpolated pollen RelAb values, the decline started around 8.5k CalYrBP in Northern Africa (median max. 3%), then at 7.5k CalYrBP in northern Italy/Po Plain (median max. 4%), in Greece at 6.5k CalYrBP (median max. 4%), in southern Italy and Central Eastern Europe at 6k CalYrBP (median max. 6% and 4.2%, respectively), in Spain and the Pyrenees at 5.5–5k CalYrBP (median max. 3.5%), and in the Alps at 4.5k CalYrBP (median max. 10%). Central Europe showed a decline at 4.5k CalYrBP. Still, the species recovered after only 1000 years, until 1k CalYrBP, when a general decline of the species also occurred in other regions that had not yet shown a decline, such as the Carpathians and the Pannonian Plain. Finally, there were some regions where the species’ RelAb is currently increasing, e.g., the Pontic region and temperate Eastern Europe.

The predicted present-day distribution of the species mainly covered the climate conditions the species occupied in the past, yet many areas suitable for the species emerged as unoccupied (Figs. 3O, 5E). Moreover, the center of the niche, defined as the part of the niche where most *Abies* spp. occurrences are found, exhibited two peaks in mean annual temperature. The niche center predicted for today was ~2.2 °C and ~3 °C colder and had more precipitation annually (by ~40 mm) than in the Early Holocene. At the same time, the niche center for the Late Holocene was ~2.2 °C and ~3.3 °C colder and had ~120 mm more precipitation annually than in the Early Holocene (Fig. 5B–E).

Finally, under future climate scenarios, the RelAb of *Abies* species was predicted to increase in the mountainous regions of Europe, decrease slightly in the Czech Republic, Slovakia, and Hungary, and remain stable in most other areas (Fig. 6). When PoO was included (obtained with or without climate), it had a major influence on the spatial distribution of RelAb, causing it to decline more rapidly the further away from sites recording *Abies*. This result is consistent with past interpolations. Additionally, the model forced only with PoO predicted a general increase of the species in the area between the sites. Among the different climate scenarios, VarSel_STRK predicted a dramatic rise in RelAb and suggests very high RelAb even in areas not naturally colonized (e.g., Scotland) (Fig. 6). On the contrary, the Varsel_Exp results were more conservative and did not predict significant changes. Additionally, two of the climate variable dataset scenarios (VarSel_Exp and VarSel_STRK) exhibited a narrower RelAb distribution, while VarSel_Surv was found to be platykurtic (Fig. 6). Finally, the warmer future climate scenario yielded a smaller increase in RelAb, while the decline did not exhibit a consistent pattern.

For a full review of *Abies* spp. demographic history see Supplementary review.

## 4. Discussion

### 4.1. Strengths and limitations of the workflow

In this study, we developed and tested a new workflow to reconstruct past demographic changes of forest tree species using pollen, simulated paleoclimate data, and geostatistical modeling. The combination of two regression kriging approaches (Fig. 1) enabled us to predict the species’ PoO and to interpolate their RelAB, ultimately providing a spatially and temporally explicit demographic history at the regional level (~32 × 32 km) on a continental scale.

Our cross-validation suggests that reconstruction at this spatial scale can be achieved with a Brier’s score of ~0.1 when sufficient data are provided. For the *Abies* genus in Europe, this threshold was reached only during the Holocene. Oriani et al. (2024) likewise found that pollen-based interpolations with kriging are robust if a sufficiently large dataset is provided. In our dataset of *Abies* spp., 50 sites with records
<figure><img src="image_10.png" /><figcaption>**Fig. 5.** Evolution of the climate niche of *Abies* spp. across time. Gray points indicate the available climate space across the continent but not occupied by *Abies* spp., while colored points indicate the interpolated relative abundance (RelAb) using probability of occurrence (calculated including climate, VarSel_IK) and climate (VarSel_STRK).</figcaption></figure>

of the species (out of 200 evenly distributed sites) appeared sufficient, which limits the applicability of the workflow to non-refugial periods. Additionally, when only a few observations are available, the climate–species relationships cannot be captured accurately; thus, climate data cannot be used to fill the gaps (Fig. 4C). This result is more evident when using shorter time bins, as the data fragmentation increases signal isolation (Fig. S7, S8, S9). This is a well-known limitation in paleobiogeography, as fossil data cannot be interpreted as the realized distribution of the species because they operate at multilevel scales that are not always fine enough (taxonomic, temporal, and spatial scales; Varela et al., 2011 and references therein). A possible solution to overcome this limitation is to apply a looser limit to determine the local presence of the species, such as the absolute limit (i.e., the first occurrence of the species vs. the continuous presence of the empirical limit used in our workflow), or to interpolate the presences (or RelAb) when only a small number of past occurrences are available (e.g., GRIWM method; Barnosky et al., 2015). This, however, would increase inaccuracy and may introduce instability in the predictions.

Previously proposed methods of past vegetation reconstruction often rely on the stationarity assumption, i.e., that the relationship between environmental conditions (such as climate) and a species' response (presence or abundance) remains constant over time (Yates et al., 2018; Chevalier et al., 2020). However, this is potentially an unrealistic simplification in ecological modeling and pollen-based paleoclimate reconstructions (Chevalier et al., 2020). Species responses to environmental conditions are not static; they can be influenced by factors such as resource availability, biotic interactions, population density, and human disturbances (Yates et al., 2018; Hou et al., 2025; Serra-Maluquer et al., 2025). The case of European *Abies* spp. (Fig. 5) serves as an example of the potential consequences of ignoring non-stationarity: modern data would not explain the observations for the LGM, Late Glacial and Early Holocene, since the largest populations of this period were present under different climate conditions than today. Thus, relying on a stationary model across such major environmental shifts could hinder or bias ecological processes or mask adaptation to

a changing environment. We addressed this by defining different eras characterized by distinct demographic patterns and using the data from the eras only for training (Fig. 1). This allowed the regression model to capture potentially different relationships between species RelAb and climate variables during distinct demographic periods (e.g., refugia or expansion), and eventually to capture different demographic processes not explained by climate in the residuals. Additionally, we trained a separate regression model for each time frame when predicting RelAb (before adding the interpolated residuals) (Fig. 1). This ensured that the model could capture a potentially unique, non-stationary relationship between *Abies* RelAb and climate variables within each specific training time frame, leading to a more dynamic and potentially accurate representation of the species' response through time.

While the presented workflow offers flexibility in implementing climate and its dynamic relationship with species, this may not lead to improvements in predicting species distributions or RelAb. Generally speaking, including climate increases the spatial heterogeneity of the interpolated PoO and RelAb, with minor regional differences across the variable sets; however, it does not necessarily increase the accuracy across the entire map (Fig. 3). This highlights a crucial point regarding migration lag and climate suitability maps. For example, in predictions of RelAb at unknown locations after the end of the Late Glacial, the highest correlation was consistently achieved when PoO (without climate) was the only explanatory variable. Conversely, model performance suffered if PoO was excluded from the model or if PoO was calculated by incorporating climate data (Fig. 4C). This suggests that, during the post-glacial recolonization by *Abies* spp., the presence of existing populations that served as sources (PoO without climate) was the main predictor of RelAb, and including climate information in the PoO may decrease model performance by confusing environmental suitability with dispersal and legacy limits. However, it is essential to consider the limitations of the reconstructions before 8k CalYrBP, most importantly, the limited number of sites available for this period. Additionally, the low temporal resolution of most time series and our applied temporal binning may have masked rapid responses, as
<figure><img src="image_11.png" /><figcaption>**Fig. 6.** Prediction of the relative abundance (RelAb) of *Abies* spp. under different future climate scenarios. (A, D, G) Density plots of *Abies* spp. RelAb across Europe. Values smaller than 0.01 were excluded. (B, C, E, F, H, I) Difference in predicted RelAb under the climate scenarios ssp126 and ssp585 (Karger et al., 2017) and the last time bin of the interpolation of the past (500–0 CalYrBP). Note that the maps are not cropped. Distribution of *Abies* spp. from Caudullo et al. (2017).</figcaption></figure>

observed in very-high-resolution records when *Abies* population dynamics are directly compared with climate at (ecologically meaningful) decadal scales (e.g., Tinner and Lotter, 2001; Lang et al., 2023). Finally, the simulated climate might be inadequate for this period, as the lack of paleo proxies for some climate variables, e.g., precipitation patterns, prevents validation of the climate model results. Future research on proxies for these climate patterns will enable more reliable vegetation-response reconstruction at the continental scale.

Conversely, in predictions of RelAb at unknown locations for the Late Holocene, the difference in results between including PoO (with or without climate) and including only climate in the interpolation was slight, with PoO alone showing the larger error and lower correlation. This suggests that climate drove RelAb in this period, rather than migration constraints or suitability. Additionally, comparing the estimated demography in the youngest time bin (i.e. the last 500 years) obtained with models including climate or not, indicates that climate increased RelAb in the mountain chains and leading edge regions, and decreased RelAb in areas not covered today by the species, such as the Pannonian Plain, the Po Plain, and the Mediterranean French coast (Fig. 4D). Thus climate has an essential role in predicting RelAb in the last time bin and it is crucial to provide additional information. We speculate that the greater role of climate relative to source population, distribution, and suitability stems from the relatively static spatial structure of the sites after millennia of establishment (Fig. 4C, S5). This prolonged period after establishment should, in theory, allow the species to spread and colonize suitable areas nearby. This contrasts with the discrepancy

between the available suitable habitats and those colonized over the last two millennia (Fig. 5), which is consistent with previous studies. Some authors attributed this significantly lower abundance, as well as a species distribution more restricted than expected based on the species niche, to human impact (e.g., Henne et al., 2011; Ruosch et al., 2016; Mauri et al., 2022; Lang et al., 2023). However, others attributed this finding to a migration lag due to competition with the existing vegetation (Sandel et al., 2025) or to the location of refugia and the slow species migration rate (Bennett et al., 1991; Svenning and Skov, 2004, 2007; Normand et al., 2011; Cheddadi et al., 2013). While our results show the expected reduction in RelAb and a halt in species expansion, we cannot disentangle the roles of climate, human impact, and sample size, given the nature of our analysis. Indeed, the small prediction error in the Late Holocene might be due to data availability.

Our workflow also allows assessment of a species' climate niche across different eras. At our rather rough spatial and temporal resolution, we found that European *Abies* spp., particularly *Abies alba*, had a large niche and was able to colonize a wide range of habitats, spanning regional mean annual temperatures of -1.4 to 16 °C and regional annual precipitation sums of 444 to 1655 mm in the raw pollen data, and from -3 to 21.2 °C and from 145 to 2200 mm in the simulated data (Fig. 5). However, these values deviate from the known niche of *Abies*, being too extreme in temperature and precipitation. This deviation is primarily caused by the low spatial and temporal resolution of pollen and, consequently, climate data, as well as by an artifact of the regression, which may over- or underpredict PoO and RelAb when predicting
outside the training frame. Indeed, *Abies* spp. are only found in mountainous environments today and are considered mountain species by the forestry community (Mauri et al., 2016). Yet our findings, in agreement with previous regional analyses of paleoecological data (Tinner et al., 2013), suggest that in the Middle Holocene, *Abies alba* was dominant in warmer environments than today, provided conditions were sufficiently humid. Although we cannot fully rule out that the observed large niche size results from disequilibrium conditions (Sandel et al., 2025), our dynamic workflow accounts for population demographic trends, thereby limiting bias arising from mismatches between the realized niche and distribution/occurrence.

In our workflow, we used correction factors to account for differences in pollen productivity among taxa (Table S2); however, we did not account for transport and deposition dynamics, as in models such as REVEALS (Sugita, 2007a,b) or STEPPS (Paciorek and McLachlan, 2009; Dawson et al., 2016). We decided to ignore these processes because we believe they are challenging to capture at our spatial scale (grid cell size: ~32 × 32 km). While general patterns of differences among taxa in pollen productivity remain robust if pollen productivity is calculated as relative to a reference species (Liu et al., 2022), which in our case is *Fagus sylvatica*, using the correction factors from Lang et al. (2023), the modeling of pollen dispersal mechanisms in the landscape requires additional information, such as wind regimes and landscape heterogeneity, which are difficult to obtain and do not necessarily hold across space and time (Liu et al., 2022; Lang et al., 2023).

The flexibility of our workflow allows several possible improvements that may be addressed in future implementations. First, including more sites could certainly improve the prediction performance and will enable an increase in spatial and temporal resolution. Finer spatial resolution could, for example, allow the estimation of non-stationary responses of forests to past, present, and future extreme events, such as droughts (Serra-Maluquer et al., 2025). Second, future research could explore the inclusion of additional covariates, such as topography, soil properties, CO₂ concentrations, a proxy for human impact, or competition with other species. Additionally, depositional sites may be associated with specific lithologies or microclimates, such as volcanic lakes, which can have unique vegetation associations. Here, we attempted to include the presence of permafrost and semi-permafrost (following Levavasseur et al., 2011) in the variable selection procedure; however, they were not significant (unpublished results). Third, non-linear models may be employed to model the pollen-climate relationship, as applied in pollen-based climate reconstructions, since species may exhibit a bimodal or multimodal response to climate variables (Chevalier et al., 2020). Finally, one could explore the vegetation dynamics of an ecosystem inferred from pollen data over long periods, rather than focusing on a single species. This would allow the assessment of climate change risks and resilience of those ecosystems, although a careful consideration of nonlinear responses and spatial heterogeneity is pivotal (Sun et al., 2022; Hou et al., 2025).

## 5. Conclusions

Fossil records offer crucial insight into how species have responded to environmental changes over time, providing richer and more quantitative information than the static “snapshot” of a species’ current range by incorporating long-term demographic dynamics and expanding the species’ known niche (Varela et al., 2009). Fossil data and quantitative reconstructions at the regional to continental scale are fundamental for understanding the processes that shape today’s biogeographical patterns, such as the distribution of genetic diversity and the impact of past climate changes on a species’ evolutionary response. Our workflow applies to any species included in a pollen record, provided sufficient sites are available and the species-specific ecology is considered. Interpolated species abundance maps can be used in several downstream analyses, such as in spatially explicit population genetic models to test different demographic scenarios (Szép et al., 2022) or as a biologically

meaningful null hypothesis for detecting adaptation to environmental changes from genetic data (He et al., 2013). Additionally, identifying refugial areas and secondary contact zones between lineages could inform the selection of provenances for adaptive forest management practices, as these areas may harbor high genetic diversity and adaptive traits.

## CRediT authorship contribution statement

**Azzurra Pistone:** Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Formal analysis, Data curation, Conceptualization. **Denis Allard:** Writing – review & editing, Methodology, Conceptualization. **Christoph Schwörer:** Writing – review & editing, Data curation. **César Morales-Molino:** Writing – review & editing, Data curation. **Willy Tinner:** Writing – review & editing, Funding acquisition, Data curation. **Katalin Csilléry:** Writing – review & editing, Writing – original draft, Supervision, Methodology, Funding acquisition, Conceptualization.

## Funding sources

An ERC CoG grant to KC (MyGardenOfTrees, 101003296) funded this project and covered the salaries of AP and KC. The University of Bern also supported AP through funding by the Canton of Bern to WT. CMM was supported by an Atracción de Talento Modalidad 1 grant (2022-T1/AMB-23788) funded by the Comunidad de Madrid region.

## Declaration of competing interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: The authors declare no conflict of interest.

## Acknowledgments

Data were obtained from the European Pollen Database (EPD) and the Alpine Palynological Database (ALPADABA) via the Neotoma Palaeoecology Database (http://www.neotomadb.org). The work of data contributors, data stewards, and the EPD, ALPADABA and Neotoma communities is gratefully acknowledged. We thank Melissa Dawes for proofreading, Moreno Bevilacqua for his suggestions on GeoModels and W.O. (Pim) van der Knaap for his valuable contribution in Lang et al. (2023) which inspired this work.

## Appendix A. Supplementary data

Supplementary figures, tables, review, and results are available on GitHub (https://github.com/azzpis/pollenkriging).

Supplementary material related to this article can be found online at https://doi.org/10.1016/j.ecoinf.2025.103561.

## Data availability

Data and code related to our findings are available on GitHub (https://github.com/azzpis/pollenkriging). Lists of data sources for raw pollen data are provided in Supplementary Table 1 and in Lang et al. (2023). Current and future CHELSA v2.1 climate data are from Karger et al. (2021). Past climate data (CHELSA-TraCE21k) are from Karger et al. (2020). The distribution of *Abies* spp. is from Caudullo et al. (2017).
## References

Aitken, S.N., Bemmels, J.B., 2015. Time to get moving: Assisted gene flow of forest trees. Evol. Appl. (ISSN: 1752-4571) 9 (1), 271–290. http://dx.doi.org/10.1111/eva.12293.

Allard, D., Clarotto, L., Emery, X., 2022. Fully nonseparable gneiting covariance functions for multivariate space-time data. Stat. Stat. (ISSN: 2211-6753) 52, 100706. http://dx.doi.org/10.1016/j.spasta.2022.100706.

Allard, D., Comunian, A., Renard, P., 2012. Probability aggregation methods in geoscience. Math. Geosci. (ISSN: 1874-8953) 44 (5), 545–581. http://dx.doi.org/10.1007/s11004-012-9396-3.

Allen, C.D., 2009. Climate-induced forest dieback: An escalating global phenomenon. Unasylva 231 (232), 60.

Ammann, B., van Raden, U.J., Schwander, J., Eicher, U., Gilli, A., Bernasconi, S.M., van Leeuwen, J.F., Lischke, H., Brooks, S.J., Heiri, O., Nováková, K., van Hardenbroek, M., von Grafenstein, U., Belmecheri, S., van der Knaap, W., Magny, M., Eugster, W., Colombaroli, D., Nielsen, E., Tinner, W., Wright, H.E., 2013. Responses to rapid warming at termination 1a at gerzensee (central Europe): Primary succession, albedo, soils, lake development, and ecological interactions. Palaeogeogr. Palaeoclimatol. Palaeoecol. (ISSN: 0031-0182) 391, 111–131. http://dx.doi.org/10.1016/j.palaeo.2013.11.009.

Barnosky, A.D., Lindsey, E.L., Villavicencio, N.A., Bostelmann, E., Hadly, E.A., Wanket, J., Marshall, C.R., 2015. Variable impact of late-quaternary megafunal extinction in causing ecological state shifts in North and South America. Proc. Natl. Acad. Sci. (ISSN: 1091-6490) 113 (4), 856–861. http://dx.doi.org/10.1073/pnas.1505295112.

Benito Garzón, M., Sánchez de Dios, R., Sáinz Ollero, H., 2007. Predictive modelling of tree species distributions on the iberian peninsula during the last glacial maximum and Mid-Holocene. Ecography (ISSN: 1600-0587) 30 (1), 120–134. http://dx.doi.org/10.1111/j.0906-7590.2007.04813.x.

Bennett, A.C., McDowell, N.G., Allen, C.D., Anderson-Teixeira, K.J., 2015. Larger trees suffer most during drought in forests worldwide. Nat. Plants (ISSN: 2055-0278) 1 (10), http://dx.doi.org/10.1038/nplants.2015.139.

Bennett, K.D., Tzedakis, P.C., Willis, J., 1991. Quaternary Refugia of north European trees. J. Biogeogr. (ISSN: 0305-0270) 18 (1), 103. http://dx.doi.org/10.2307/2845248.

Bevilacqua, M., Morales-Oñate, V., Caamaño-Carrillo, C., 2024. GeoModels: Procedures for Gaussian and non Gaussian geostatistical (large) data analysis. R package version 2.1.0.

Birks, H.H., Ammann, B., 2000. Two terrestrial records of rapid climatic change during the Glacial–Holocene transition (14,00020139,000 calendar years B.P.) from Europe. Proc. Natl. Acad. Sci. (ISSN: 1091-6490) 97 (4), 1390–1394. http://dx.doi.org/10.1073/pnas.97.4.1390.

Bledý, M., Vacek, S., Brabec, M., Vacek, Z., Zukor, J., Černý, J., Ševčík, R., Brynchová, K., 2024. Silver fir (Abies alba Mill.): Review of ecological insights, forest management strategies, and climate change's impact on European forests. Forests (ISSN: 1999-4907) 15 (6), 998. http://dx.doi.org/10.3390/f15060998.

Caudullo, G., Welk, E., San-Miguel-Ayanz, J., 2017. Chorological maps for the main European woody species. Data Brief (ISSN: 2352-3409) 12, 662–666. http://dx.doi.org/10.1016/j.dib.2017.05.007.

Cheddadi, R., Birks, H.J.B., Tarroso, P., Liepelt, S., Gömöry, D., Dullinger, S., Meier, E.S., Hülber, K., Maiorano, L., Laborde, H., 2013. Revisiting tree-migration rates: Abies alba (Mill.), a case study. Veg. Hist. Archaeobotany (ISSN: 1617-6278) 23 (2), 113–122. http://dx.doi.org/10.1007/s00334-013-0404-4.

Chevalier, M., Davis, B.A., Heiri, O., Seppä, H., Chase, B.M., Gajewski, K., Lacourse, T., Telford, R.J., Finsinger, W., Guiot, J., Kühl, N., Maezumi, S.Y., Tipton, J.R., Carter, V.A., Brussel, T., Phelps, L.N., Dawson, A., Zanon, M., Vallé, F., Nolan, C., Mauri, A., de Vernal, A., Izumi, K., Holmström, L., Marsicek, J., Goring, S., Sommer, P.S., Chaput, M., Kupriyanov, D., 2020. Pollen-based climate reconstruction techniques for late quaternary studies. Earth-Sci. Rev. (ISSN: 0012-8252) 210, 103384. http://dx.doi.org/10.1016/j.earscirev.2020.103384.

Chiles, J.-P., Delfiner, P., 2012. Geostatistics: Modeling Spatial Uncertainty. John Wiley & Sons.

Davis, M.B., Shaw, R.G., Etterson, J.R., 2005. Evolutionary responses to changing climate. Ecology (ISSN: 0012-9658) 86 (7), 1704–1714. http://dx.doi.org/10.1890/03-0788.

Dawson, A., Paciorek, C.J., McLachlan, J.S., Goring, S., Williams, J.W., Jackson, S.T., 2016. Quantifying pollen-vegetation relationships to reconstruct ancient forests using 19th-century forest composition and pollen data. Quat. Sci. Rev. (ISSN: 0277-3791) 137, 156–175. http://dx.doi.org/10.1016/j.quascirev.2016.01.012.

Deza-Araujo, M., Morales-Molino, C., Conedera, A., Henne, P.D., Krebs, P., Hinz, M., Heitz, C., Hafner, A., Tinner, W., 2022. A new indicator approach to reconstruct agricultural land use in Europe from sedimentary pollen assemblages. Palaeogeogr. Palaeoclimatol. Palaeoecol. (ISSN: 0031-0182) 599, 111051. http://dx.doi.org/10.1016/j.palaeo.2022.111051.

Dubrule, O., 2017. Indicator variogram models: Do we have much choice? Math. Geosci. (ISSN: 1874-8953) 49 (4), 441–465. http://dx.doi.org/10.1007/s11004-017-9678-x.

Elith, J., Kearney, M., Phillips, S., 2010. The art of modelling range-shifting species. Methods Ecol. Evol. (ISSN: 2041-210X) 1 (4), 330–342. http://dx.doi.org/10.1111/j.2041-210x.2010.00036.x.

Erdtman, G., 1969. Handbook of Palynology: Morphology - Taxonomy - Ecology: An Introduction to the Study of Pollen Grains and Spores. Hafner Publishing Co.

Finsinger, W., Giesecke, T., Brewer, S., Leydet, M., 2017. Emergence patterns of novelty in European vegetation assemblages over the past 15,000 years. In: Williams, J. (Ed.), Ecol. Lett. (ISSN: 1461-0248) 20 (3), 336–346. http://dx.doi.org/10.1111/ele.12731.

Fyfe, R.M., de Beaulieu, J.-L., Binney, H., Bradshaw, R.H.W., Brewer, S., Le Flao, A., Finsinger, W., Gaillard, M.-J., Giesecke, T., Gil-Romera, G., Grimm, E.C., Huntley, B., Kunes, P., Kühl, N., Leydet, M., Lotter, A.F., Tarasov, P.E., Tonkov, S., 2009. The European pollen database: Past efforts and current activities. Veg. Hist. Archaeobotany (ISSN: 1617-6278) 18 (5), 417–424. http://dx.doi.org/10.1007/s00334-009-0215-9.

Genuer, R., Poggi, J.-M., Tuleau-Malot, C., 2013. VSURF: Variable selection using random forests. Title of the publication associated with this dataset: CRAN: Contributed Packages.

Giesecke, T., Davis, B., Brewer, S., Finsinger, W., Wolters, S., Blaauw, M., de Beaulieu, J.-L., Binney, H., Fyfe, R.M., Gaillard, M.-J., Gil-Romera, G., van der Knaap, W.O., Kuneš, P., Kühl, N., van Leeuwen, J.F.N., Leydet, M., Lotter, A.F., Ortu, E., Semmler, M., Bradshaw, R.H.W., 2013. Towards mapping the late quaternary vegetation change of europe. Veg. Hist. Archaeobotany (ISSN: 1617-6278) 23 (1), 75–86. http://dx.doi.org/10.1007/s00334-012-0390-y.

Giesecke, T., Wolters, S., van Leeuwen, J.F.N., van der Knaap, P.W.O., Leydet, M., Brewer, S., 2019. Postglacial change of the floristic diversity gradient in Europe. Nat. Commun. (ISSN: 2041-1723) 10 (1), http://dx.doi.org/10.1038/s41467-019-13233-y.

Githumbi, E., Fyfe, R., Gaillard, M.-J., Trondman, A.-K., Mazier, F., Nielsen, A.-B., Poska, A., Sugita, S., Woodbridge, J., Azua, J., Feurdean, A., Grindean, R., Lebreton, V., Marquer, L., Nebout-Combourieu, N., Stančikaitė, M., Tančău, I., Tonkov, S., Shumilovskikh, L., LandclimII data contributors. 2022. European pollen-based REVEALS land-cover reconstructions for the holocene: Methodology, mapping and potentials. Earth Syst. Sci. Data (ISSN: 1866-3516) 14 (4), 1581–1619. http://dx.doi.org/10.5194/essd-14-1581-2022.

Hammond, W.M., Williams, A.P., Abatzoglou, J.T., Adams, H.D., Klein, T., López, R., Sáenz-Romero, C., Hartmann, H., Breshears, D.D., Allen, C.D., 2022. Global field observations of tree die-off reveal hotter-drought fingerprint for earth's forests. Nat. Commun. (ISSN: 2041-1723) 13 (1), http://dx.doi.org/10.1038/s41467-022-29289-2.

Hartmann, H., Bastos, A., Das, A.J., Esquivel-Muelbert, A., Hammond, W.M., Martínez-Vilalta, J., McDowell, N.G., Powers, J.S., Pugh, T.A., Ruthrof, K.X., Allen, C.D., 2022. Climate change risks to global forest health: Emergence of unexpected events of elevated tree mortality worldwide. Annu. Rev. Plant Biol. (ISSN: 1545-2123) 73 (1), 673–702. http://dx.doi.org/10.1146/annurev-arplant-102820-021820.

He, Q., Edwards, D.L., Knowles, L.L., 2013. Integrative testing of how environments from the past to the present shape genetic structure across landscapes. Evolution (ISSN: 0014-3820) 67 (12), 3386–3402. http://dx.doi.org/10.1111/evo.12159.

Hengl, T., Heuvelink, G.B., Stein, A., 2004. A generic framework for spatial prediction of soil variables based on regression-kriging. Geoderma (ISSN: 0016-7061) 120 (1–2), 75–93. http://dx.doi.org/10.1016/j.geoderma.2003.08.018.

Henne, P.D., Elkin, C.M., Reineking, B., Bugmann, H., Tinner, W., 2011. Did soil development limit spruce (Picea Abies) expansion in the central alps during the Holocene? Testing a palaeobotanical hypothesis with a dynamic landscape model. J. Biogeogr. (ISSN: 0305-0270) 38 (5), 933–949. http://dx.doi.org/10.1111/j.1365-2699.2010.02460.x.

Hijmans, R.J., 2020. terra: Spatial data analysis. http://dx.doi.org/10.32614/cran package.terra.

Hou, L.-F., Li, L., Chen, R., Wu, Y.-P., Feng, G.-L., Sun, G.-Q., 2025. Vegetation dynamics: Modeling, mechanisms, and emergent properties. Phys. Rep. (ISSN: 0370-1573) 1145, 1–87. http://dx.doi.org/10.1016/j.physrep.2025.09.003.

Hu, D., Shu, H., Hu, H., Xu, J., 2017. Spatiotemporal regression kriging to predict precipitation using time-series MODIS data. Clust. Comput. (ISSN: 1573-7543) 20 (1), 347–357. http://dx.doi.org/10.1007/s10586-016-0708-0.

Kaltenrieder, P., Belis, C.A., Hofstetter, S., Ammann, B., Ravazzi, C., Tinner, W., 2009. Environmental and climatic conditions at a potential glacial refugial site of tree species near the Southern Alpine glaciers. New insights from multiproxy sedimentary studies at Lago della Costa (Euganean hills, Northeastern Italy). Quat. Sci. Rev. (ISSN: 0277-3791) 28 (25–26), 2647–2662. http://dx.doi.org/10.1016/j.quascirev.2009.05.025.

Karger, D.N., Conrad, O., Böhner, J., Kawohl, T., Kreft, H., Soria-Auza, R.W., Zimmermann, N.E., Linder, H.P., Kessler, M., 2017. Climatologies at high resolution for the earth's land surface areas. Sci. Data (ISSN: 2052-4463) 4 (1), http://dx.doi.org/10.1038/sdata.2017.122.

Karger, D.N., Conrad, O., Böhner, J., Kawohl, T., Kreft, H., Soria-Auza, R.W., Zimmermann, N.E., Linder, H.P., Kessler, M., 2021. Climatologies at high resolution for the earth's land surface areas. http://dx.doi.org/10.16904/envidat.228, EnviDat.

Karger, D.N., Nobis, M.P., Normand, S., Graham, C.H., Zimmermann, N.E., 2020. CHELSA-TraCE21k: downscaled transient temperature and precipitation data since the last glacial maximum. http://dx.doi.org/10.16904/envidat.211, EnviDat.

Karger, D.N., Nobis, M.P., Normand, S., Graham, C.H., Zimmermann, N.E., 2023. CHELSA-TraCE21k – high-resolution (1 km) downscaled transient temperature and precipitation data since the last glacial maximum. Clim. Past (ISSN: 1814-9332) 19 (2), 439–456. http://dx.doi.org/10.5194/cp-19-439-2023.
Kremer, A., Ronce, O., Robledo-Arnuncio, J.J., Guillaume, F., Bohrer, G., Nathan, R., Bridle, J.R., Gomulkiewicz, R., Klein, E.K., Ritland, K., Kuparinen, A., Gerber, S., Schueler, S., 2012. Long-distance gene flow and adaptation of forest trees to rapid climate change. Ecol. Lett. (ISSN: 1461-0248) 15 (4), 378–392. http://dx.doi.org/10.1111/j.1461-0248.2012.01746.x.

Kuhn, Max, 2008. Building predictive models in R using the caret package. J. Stat. Softw. 28, 1–26. http://dx.doi.org/10.18637/jss.v028.i05.

Kursa, M.B., Rudnicki, W.R., 2010. Feature selection with the Boruta package. J. Stat. Softw. 36 (11), 1–13.

Lang, G., Ammann, B., Behre, K.-E., Tinner, W., 2023. Quaternary Vegetation Dynamics of Europe. Haupt Verlag.

Lebourgeois, F., Gomez, N., Pinto, P., Mérian, P., 2013. Mixed stands reduce Abies alba tree-ring sensitivity to summer drought in the Vosges mountains, western Europe. Forest Ecol. Manag. (ISSN: 0378-1127) 303, 61–71. http://dx.doi.org/10.1016/j.foreco.2013.04.003.

Levavasseur, G., Vrac, M., Roche, D.M., Paillard, D., Martin, A., Vandenberghe, J., 2011. Present and LGM permafrost from climate simulations: contribution of statistical downscaling. Clim. Past (ISSN: 1814-9332) 7 (4), 1225–1246. http://dx.doi.org/10.5194/cp-7-1225-2011.

Li, J., Heap, A.D., 2011. A review of comparative studies of spatial interpolation methods in environmental sciences: Performance and impact factors. Ecol. Inform. (ISSN: 1574-9541) 6 (3–4), 228–241. http://dx.doi.org/10.1016/j.ecoinf.2010.12.003.

Liaw, A., Wiener, M., 2002. Classification and regression by randomforest. R News 2 (3), 18–22.

Liepelt, S., Cheddadi, R., de Beaulieu, J.-L., Fady, B., Gömöry, D., Hussendörfer, E., Konnert, M., Litt, T., Longauer, R., Terhürne-Berson, R., Ziegenhagen, B., 2009. Postglacial range expansion and its genetic imprints in Abies alba (Mill.) 2014 a synthesis from palaeobotanic and genetic data. Rev. Palaeobot. Palynol. (ISSN: 0034-6667) 153 (1–2), 139–149. http://dx.doi.org/10.1016/j.revpalbo.2008.07.002.

Lindner, M., Maroschek, M., Netherer, S., Kremer, A., Barbati, A., García-Gonzalo, J., Seidl, R., Delzon, S., Corona, P., Kolström, M., Leker, M.J., Marchetti, M., 2010. Climate change impacts, adaptive capacity, and vulnerability of European forest ecosystems. Forest Ecol. Manag. (ISSN: 0378-1127) 259 (4), 698–709. http://dx.doi.org/10.1016/j.foreco.2009.09.023.

Liu, Y., Ogle, K., Lichstein, J.W., Jackson, S.T., 2022. Estimation of pollen productivity and dispersal: How pollen assemblages in small lakes represent vegetation. Ecol. Monograph. (ISSN: 1557-7015) 92 (3), http://dx.doi.org/10.1002/ecm.1513.

Magri, D., 2007. Patterns of post-glacial spread and the extent of glacial refugia of European beech (Fagus sylvatica). J. Biogeogr. (ISSN: 1365-2699) 35 (3), 450–463. http://dx.doi.org/10.1111/j.1365-2699.2007.01803.x.

Masson-Delmotte, V., Zhai, P., Pirani, A., Connors, S.L., Péan, C., Berger, S., Caud, N., Chen, Y., Goldfarb, L., Gomis, M., et al., 2021. Climate change 2021: The physical science basis. p. 2391, Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change 2.1.

Mauri, A., De Rigo, D., Caudullo, G., 2016. Abies Alba in Europe: Distribution, Habitat, Usage and Threats. Publication Office of the European Union Luxembourg, pp. 48–49.

Mauri, A., Girardello, M., Strona, G., Beck, P.S.A., Forzieri, G., Caudullo, G., Manca, F., Cescatti, A., 2022. EU-Trees4F, a dataset on the future distribution of European tree species. Sci. Data (ISSN: 2052-4463) 9 (1), http://dx.doi.org/10.1038/s41597-022-01128-5.

Milborrow, S., Hastie, T., Tibshirani, R., 2007. Earth: Multivariate adaptive regression splines. http://dx.doi.org/10.32614/cran/package.earth.

Mottl, O., Flantua, S.G.A., Bhatta, K.P., Felde, V.A., Giesecke, T., Goring, S., Grimm, E.C., Haberle, S., Hooghiemstra, H., Ivory, S., Kuneš, P., Wolters, S., Seddon, A.W.R., Williams, J.W., 2021. Global acceleration in rates of vegetation change over the past 18,000 years. Science (ISSN: 1095-9203) 372 (6544), 860–864. http://dx.doi.org/10.1126/science.abg1685.

Normand, S., Ricklefs, R.E., Skov, F., Bladt, J., Tackenberg, O., Svenning, J.-C., 2011. Postglacial migration supplements climate in determining plant species ranges in Europe. Proc. R. Soc. B (ISSN: 1471-2954) 278 (1725), 3644–3653. http://dx.doi.org/10.1098/rspb.2010.2769.

Oriani, F., Mariethoz, G., Chevalier, M., 2024. EUPollMap: The European atlas of contemporary pollen distribution maps derived from an integrated kriging interpolation approach. Earth Syst. Sci. Data (ISSN: 1866-3516) 16 (1), 731–742. http://dx.doi.org/10.5194/essd-16-731-2024.

Paciorek, C.J., McLachlan, J.S., 2009. Mapping ancient forests: Bayesian inference for spatio-temporal trends in forest composition using the fossil pollen proxy record. J. Amer. Statist. Assoc. (ISSN: 1537-274X) 104 (486), 608–622. http://dx.doi.org/10.1198/jasa.2009.0026.

Pecchi, M., Marchi, M., Burton, V., Giannetti, F., Moriondo, M., Bernetti, I., Bindi, M., Chirici, G., 2019. Species distribution modelling to support forest management. A literature review. Ecol. Model. (ISSN: 0304-3800) 411, 108817. http://dx.doi.org/10.1016/j.ecolmodel.2019.108817.

Prentice, C., 1988. Records of vegetation in time and space: the principles of pollen analysis. In: Vegetation History. Springer, pp. 17–42.

Probst, P., 2018. VarImp: RF variable importance for arbitrary measures. Title of the publication associated with this dataset: CRAN: Contributed Packages.

R Core Team, 2024. R: A Language and Environment for Statistical Computing.

Rufibach, K., 2010. Use of brier score to assess binary predictions. J. Clin. Epidemiol. (ISSN: 0895-4356) 63 (8), 938–939. http://dx.doi.org/10.1016/j.jclinepi.2009.11.009.

Ruosch, M., Spahni, R., Joos, F., Henne, P.D., van der Knaap, W.O., Tinner, W., 2016. Past and future evolution of Abies alba forests in Europe 2013 comparison of a dynamic vegetation model with palaeo data and observations. Global Change Biol. (ISSN: 1365-2486) 22 (2), 727–740. http://dx.doi.org/10.1111/gcb.13075.

Sandel, B., Merow, C., Serra-Diaz, J.M., Svenning, J.-C., 2025. Disequilibrium in plant distributions: Challenges and approaches for species distribution models. J. Ecol. (ISSN: 1365-2745) 113 (4), 782–794. http://dx.doi.org/10.1111/1365-2745.70009.

Schild, L., Ewald, P., Li, C., Hébert, R., Laepple, T., Herzschuh, U., 2025. Legacy Vegetation: Northern hemisphere reconstruction of past plant cover and total tree cover from pollen archives of the last 14 kyr. Earth Syst. Sci. Data (ISSN: 1866-3516) 17 (5), 2007–2033. http://dx.doi.org/10.5194/essd-17-2007-2025.

Sekulić, A., Kilibarda, M., Protić, D., Tadić, M.P., Bajat, B., 2019. Spatio-temporal regression kriging model of mean daily temperature for Croatia. Theor. Appl. Climatol. (ISSN: 1434-4483) 140 (1–2), 101–114. http://dx.doi.org/10.1007/s00704-019-03077-3.

Serge, M., Mazier, F., Fyfe, R., Gaillard, M.-J., Klein, T., Lagnoux, A., Galop, D., Githumbi, E., Mindrescu, M., Nielsen, A., Trondman, A.-K., Poska, A., Sugita, S., Woodbridge, J., Abel-Schaad, D., Åkesson, C., Alenius, T., Ammann, B., Andersen, S., Anderson, R., Andrič, M., Balakauskas, L., Barnekow, L., Batalova, V., Bergman, J., Birks, H., Björkman, L., Bjune, A., Borisova, O., Broothaerts, N., Carrion, J., Caseldine, C., Christiansen, J., Cui, Q., Currás, A., Czerwiński, S., David, R., Davies, A., De Jong, R., Di Rita, F., Dietre, B., Dörfler, W., Doyen, E., Edwards, K., Ejarque, A., Endtmann, E., Etienne, D., Faure, E., Feeser, I., Feurdean, A., Fischer, E., Fletcher, W., Franco-Mágica, F., Fredh, E., Froyd, C., Garcés-Pastor, S., García-Moreiras, I., Gauthier, E., Gil-Romera, G., González-Sampériz, P., Grant, M., Grindean, R., Haas, J., Hannon, G., Heather, A.-J., Heikkilä, M., Hjelle, K., Jahns, S., Jasiunas, N., Jiménez-Moreno, G., Jouffroy-Bapicot, I., Kabailiené, M., Kamerling, I., Kangur, M., Karpínska-Kolaczek, M., Kasianova, A., Kolaczek, P., Lagerås, P., Latalowa, M., Lechterbeck, J., Leroyer, C., Leydet, M., Lindbladh, M., Lisitsyna, O., López-Sáez, J.-A., Lowe, J., Luelmo-Lautenschlaeger, R., Lukanina, E., Macijauskaitė, L., Magri, D., Marguerie, D., Marquer, J., Martínez-Cortizas, M., Mehl, I., Mesa-Fernández, J., Mighall, T., Miola, A., Miras, Y., Morales-Molino, C., Mrotzek, A., Sobrino, C., Odgaard, B., Ozola, I., Pérez-Díaz, S., Pérez-Oblol, R., Poggi, C., Rego, P., Ramos-Román, M., Rasmussen, P., Reille, M., Rösch, M., Ruffaldi, P., Goni, M., Savukynienė, N., Schröder, T., Schult, M., Segerström, U., Seppä, H., Vives, G., Shumilovskikh, L., Smettan, H., Stancikaitė, M., Stevenson, A., Stivrins, N., Tantau, I., Theuerkauf, M., Tonkov, S., van der Knaap, W., van Leeuwen, J., Vecmane, E., Verstraeten, G., Veski, S., Voigt, R., Von Stedingk, H., Waller, M., Wiethold, J., Willis, K., Wolters, S., Zernitskaya, V., 2023. Testing the effect of relative pollen productivity on the reveals model: A validated reconstruction of europe-wide holocene vegetation. Land (ISSN: 2073-445X) 12 (5), 986. http://dx.doi.org/10.3390/land12050986.

Serra-Maluquer, X., Astigarraga, J., Morales-Molino, C., Ruiz-Benito, P., 2025. Non-stationary forest responses to hotter droughts: a temporal perspective considering the role of past legacies. Ecography (ISSN: 1600-0587) http://dx.doi.org/10.1002/ecog.07740.

Sturrock, R.N., Frankel, S.J., Brown, A.V., Hennon, P.E., Kliejunas, J.T., Lewis, K.J., Worrall, J.J., Woods, A.J., 2011. Climate change and forest diseases. Plant Pathol. (ISSN: 1365-3059) 60 (1), 133–149. http://dx.doi.org/10.1111/j.1365-3059.2010.02406.x.

Sugita, S., 2007a. Theory of quantitative reconstruction of vegetation I: Pollen from large sites REVEALS regional vegetation composition. Holocene (ISSN: 1477-0911) 17 (2), 229–241. http://dx.doi.org/10.1177/0959683607075837.

Sugita, S., 2007b. Theory of quantitative reconstruction of vegetation II: All you need is LOVE. Holocene (ISSN: 1477-0911) 17 (2), 243–257. http://dx.doi.org/10.1177/0959683607075838.

Sun, G.-Q., Li, L., Li, J., Liu, C., Wu, Y.-P., Gao, S., Wang, Z., Feng, G.-L., 2022. Impacts of climate change on vegetation pattern: Mathematical modeling and data analysis. Phys. Life Rev. (ISSN: 1571-0645) 43, 239–270. http://dx.doi.org/10.1016/j.prev.2022.09.005.

Svenning, J.-C., Fleiggaard, C., Marske, K.A., Nógues-Bravo, D., Normand, S., 2011. Applications of species distribution modeling to paleobiology. Quat. Sci. Rev. (ISSN: 0277-3791) 30 (21–22), 2930–2947. http://dx.doi.org/10.1016/j.quascirev.2011.06.012.

Svenning, J.-C., Skov, F., 2004. Limited filling of the potential range in European tree species. Ecol. Lett. (ISSN: 1461-0248) 7 (7), 565–573. http://dx.doi.org/10.1111/j.1461-0248.2004.00614.x.

Svenning, J.-C., Skov, F., 2007. Could the tree diversity pattern in Europe be generated by postglacial dispersal limitation? Ecol. Lett. (ISSN: 1461-0248) 10 (6), 453–460. http://dx.doi.org/10.1111/j.1461-0248.2007.01038.x.

Szép, E., Trubenová, B., Csilléry, K., 2022. Using <scp>gridcoal</scp> to assess whether standard population genetic theory holds in the presence of spatio-temporal heterogeneity in population size. Mol. Ecol. Resour. (ISSN: 1755-0998) 22 (8), 2941–2955. http://dx.doi.org/10.1111/1755-0998.13676.

Theuerkauf, M., Couwenberg, J., 2018. ROPES reveals past land cover and PPEs from single pollen records. Front. Earth Sci. (ISSN: 2296-6463) 6, http://dx.doi.org/10.3389/feart.2018.00014.
Tinner, W., Colombaroli, D., Heiri, O., Henne, P.D., Steinacher, M., Untenecker, J., Vescovi, E., Allen, J.R.M., Carraro, G., Conedera, M., Joos, F., Lotter, A.F., Luterbacher, J., Samartin, S., Valsecchi, V., 2013. The past ecology of Abies alba provides new perspectives on future responses of silver fir forests to global warming. Ecol. Monograph. (ISSN: 1557-7015) 83 (4), 419–439. http://dx.doi.org/10.1890/12-2231.1.

Tinner, W., Lotter, A.F., 2001. Central European vegetation response to abrupt climate change at 8.2 ka. Geology 29 (6), 551–554.

Trachsel, M., Dawson, A., Paciorek, C.J., Williams, J.W., McLachlan, J.S., Cogbill, C.V., Foster, D.R., Goring, S.J., Jackson, S.T., Oswald, W.W., Shuman, B.N., 2020. Comparison of settlement-era vegetation reconstructions for STEPPS and REVEALS pollen 2013vegetation models in the northeastern United States. Quat. Res. (ISSN: 1096-0287) 95, 23–42. http://dx.doi.org/10.1017/qua.2019.81.

van Zoest, V., Osei, F.B., Hoek, G., Stein, A., 2019. Spatio-temporal regression kriging for modelling urban NO2 concentrations. Int. J. Geogr. Inf. Sci. (ISSN: 1362-3087) 34 (5), 851–865. http://dx.doi.org/10.1080/13658816.2019.1667501.

Varela, S., Lobo, J.M., Hortal, J., 2011. Using species distribution models in paleobio-
geography: A matter of data, predictors and concepts. Palaeogeogr. Palaeoclimatol.
Palaeoecol. (ISSN: 0031-0182) 310 (3–4), 451–463. http://dx.doi.org/10.1016/j.
palaeo.2011.07.021.

Varela, S., Rodríguez, J., Lobo, J.M., 2009. Is current climatic equilibrium a guarantee for the transferability of distribution model predictions? A case study of the spotted hyena. J. Biogeogr. (ISSN: 1365-2699) 36 (9), 1645–1655. http://dx.doi.org/10.1111/j.1365-2699.2009.02125.x.

Vitali, V., Büntgen, U., Bauhus, J., 2017. Silver fir and Douglas fir are more tolerant to extreme droughts than Norway spruce in south-western Germany. Global Change Biol. (ISSN: 1365-2486) 23 (12), 5108–5119. http://dx.doi.org/10.1111/gcb.13774.

von Post, L., 1916. Mötet den 2 november 1916. Geol. Föreningen I Stock. Förhan-
dlingar 38 (6), 383–394. http://dx.doi.org/10.1080/11035891609444403, eprint:
https://doi.org/10.1080/11035891609444403.

Wermelinger, B., 2004. Ecology and management of the spruce bark beetle ips
typographus—a review of recent research. Forest Ecol. Manag. (ISSN: 0378-1127)
202 (1–3), 67–82. http://dx.doi.org/10.1016/j.foreco.2004.07.018.

Williams, J.W., Grimm, E.C., Blois, J.L., Charles, D.F., Davis, E.B., Goring, S.J., Graham, R.W., Smith, A.J., Anderson, M., Arroyo-Cabrales, J., Ashworth, A.C., Betancourt, J.L., Bills, B.W., Booth, R.K., Buckland, P.I., Curry, B.B., Giesecke, T., Jackson, S.T., Latorre, C., Nichols, J., Purdum, T., Roth, R.E., Stryker, M., Takahara, H., 2018. The neotoma paleoecology database, a multiproxy, international, community-curated data resource. Quat. Res. (ISSN: 1096-0287) 89 (1), 156–177. http://dx.doi.org/10.1017/qua.2017.105.

Yates, K.L., Bouchet, P.J., Caley, M.J., Mengersen, K., Randin, C.F., Parnell, S., Fielding, A.H., Bamford, A.J., Ban, S., Barbosa, A.M., Dormann, C.F., Elith, J., Embling, C.B., Ervin, G.N., Fisher, R., Gould, S., Graf, R.F., Gregr, E.J., Halpin, P.N., Heikkinen, R.K., Heinänen, S., Jones, A.R., Krishnakumar, P.K., Lauria, V., Lozano-Montes, H., Mannocci, L., Mellin, C., Mesgaran, M.B., Moreno-Amat, E., Mormede, S., Novaczek, E., Oppel, S., Ortuño Crespo, G., Peterson, A.T., Rapaciuolo, G., Roberts, J.J., Ross, R.E., Scales, K.L., Schoeman, D., Snelgrove, P., Sundblad, G., Thuiller, W., Torres, L.G., Verbruggen, H., Wang, L., Wenger, S., Whittingham, M.J., Zharikov, Y., Zurell, D., Sequeira, A.M., 2018. Outstanding challenges in the transferability of ecological models. Trends Ecol. Evolut. (ISSN: 0169-5347) 33 (10), 790–802. http://dx.doi.org/10.1016/j.tree.2018.08.001.