Multi-decadal trend analysis and forest disturbance assessment of European tree species: concerning signs of a subtle shift

Climate change poses a significant threat to the distribution and composition of forest tree species worldwide. European forest tree species ’ range is expected to shift to cope with the increasing frequency and intensity of extreme weather events, pests and diseases caused by climate change. Despite numerous regional studies, a continental scale assessment of current changes in species distributions in Europe is missing due to the difficult task of modeling a species realized distribution and to quantify the influence of forest disturbances on each species. In this study we conducted a trend analysis on the realized distribution of 6 main European forest tree species ( Abies alba Mill., Fagus sylvatica L., Picea abies L. H. Karst., Pinus nigra J. F. Arnold, Pinus sylvestris L. and Quercus robur L.) to capture and map the prevalent trends in probability of occurrence for the period 2000 – 2020. We also analyzed the impact of forest disturbances on each species ’ range and identified the dominant disturbance drivers. Our results revealed an overall trend of stability in species ’ distributions (85% of the pixels are considered stable by 2020 for all species) but we also identified some hot spots characterized by negative trends in probability of occurrence, mostly at the edges of each species ’ latitudinal range. Additionally, we identified a steady increase in disturbance events in each species ’ range by disturbance (affected range doubled by 2020, from 3.5% to 7% on average) and highlighted species-specific responses to forest disturbance drivers such as wind and fire. Overall, our study provides insights into distribution trends and disturbance patterns for the main European forest tree species. The identification of range shifts and the intensifying impacts of disturbances call for proactive conservation efforts and long-term planning to ensure the resilience and sustainability of European forests.


Introduction
Tree species have repeatedly demonstrated remarked adaptability over time, changing their geographical distributions in response to large scale fluctuations in the climate (Cheddadi et al., 2016;San-Miguel-Ayanz et al., 2016;Svenning et al., 2010;Svenning and Skov, 2007).Over centuries up to millennia, these adjustments have highlighted the delicate balance between species survival and the environmental conditions that either support or negate their existence.Owing to their long generation time, trees are however particularly sensitive to climate change.The current, unprecedented rapid climate change has the potential to drive major shifts in tree species distributions (Hanberry and Hansen, 2015;IPCC, 2022) to the point that it could alter the global location of entire biomes (Berner and Goetz, 2022;Bonannella et al., 2023;Gonzalez et al., 2010).
Two decades of studies have significantly improved the detection and attribution of biome shifts to climate change (Higgins et al., 2023;Lindner et al., 2010).Numerous studies address the expected shifts in distributions of tree species using Species Distribution Models (SDMs) (Franklin, 2010), which map the ecological niche of the species in space under current conditions, and then include different climatic scenarios projections to predict the geographical boundaries of the niche in the future (Dyderski et al., 2017;Mauri et al., 2022;Thuiller et al., 2008;Zhang et al., 2017).There is general consensus regarding the response of species' ranges to climate change, pointing towards poleward (Berner and Goetz, 2022;Zhu et al., 2014Zhu et al., , 2012) ) and upward (Feeley et al., 2011) shifts in distribution ranges.Studies using tools like SDMs have mainly focused on predicting potential species distributions.In contrast, modeling the realized distribution is an inherently complex task and extrapolating this model into the future is even more challenging.Despite its complexity, modeling the realized distribution in current conditions may still be feasible, for example using Earth Observation (EO) data, but there are only few examples of such exercises (Bonannella et al., 2022;Gelfand and Shirota, 2021;Hefley and Hooten, 2016).Compared to previous studies focusing on long term future predictions using climate models, using EO data based models allows for timely and consistent assessments of ongoing changes (Strickland et al., 2020): EO data does not only provide higher temporal resolution compared to, for example, periodic surveys such as National Forest Inventories (NFIs) programs, but also offers a broader spatial coverage, including remote or inaccessible regions which are difficult to monitor from the ground.EO data also offers a compelling advantage to extensive field surveys in terms of cost efficiency: accessing and processing EO data is nowadays remarkably streamlined due to the free availability of large collections of data, like the Landsat mission (Zhu et al., 2019), or tools able to analyze such large collections like the cloud computing platform Google Earth Engine (Gorelick et al., 2017).While field campaigns like NFI programs have been used to detail current ongoing changes in species distributions at regional or national scale (Ewald, 2012;Rigling et al., 2013;Scherrer et al., 2022), the decadal temporal resolution, the limited geographic extent (national at most) and access (only few users can actually harvest all the useful information these collections of plot data offer) of these surveys makes their use for tree species shifts assessment challenging (Nabuurs et al., 2022).Due to these reasons, a comprehensive understanding of changes in tree species distributions trends at large (i.e., continental) scale in the face of climate change is still lacking.
Particularly European forests have recently faced a barrage of severe threats as an effect of previous forest management decisions and climate change (Senf and Seidl, 2021a).These disturbances have increased consistently over the last century (Patacca et al., 2023;Schelhaas et al., 2003;Seidl et al., 2011) with the most pronounced intensification occurring in the last two decades (Senf et al., 2018).This intensification has caused extensive tree loss and drastic changes in the structure and dynamics of European forests (Maes et al., 2023;Popkin, 2021;Senf and Seidl, 2021a).Several studies identified positive correlations between tree species diversity and abundance and the provision of ecosystem services (Brockerhoff et al., 2017;Gamfeldt et al., 2013;Himes et al., 2020).Hence, their decline not only negatively affects the environment but human societies as well (FAO, 2022;FOREST EUROPE, 2020).Given the importance of the topic, numerous studies have analyzed European forest disturbances qualitatively (Forzieri et al., 2021;Seidl et al., 2017;Sommerfeld et al., 2018), quantitatively (Patacca et al., 2023) and spatially (Senf and Seidl, 2021a;b).While some researchers focused on dissecting specific disturbance events, seeking to unravel the causes or the impacts of individual storms (Chirici et al., 2019;Kronauer, 2007), fires (Ganteaume et al., 2013;San-Miguel-Ayanz et al., 2012) or pests (Hlásny et al., 2021a), others analyzed the species-specific effects of these disturbances.Schmidt et al. (2010) tried to estimate the risk of storm damage for different groups of species, spruce and Scots pine in particular in Southwestern Germany; Hlásny et al. (2021b) and Kautz et al. (2023) analyzed drivers and symptoms of European spruce bark beetle outbreaks on the spruce in Central Europe and Moris et al. (2023) tested the resilience of beech stands to fire in Northern Italy.These localized or context-specific studies offer valuable insights into the species' behavior in response to disturbances, especially if the study area is at the edges of the species range: in this case, they can provide valuable knowledge in understanding what the species range response could be in the future, if it will advance or recede.While insightful, these studies primarily concerned localized (regional or national) contexts; in contrast, the pivotal aspect of our study is to explore disturbances at a continental level, matching the scale of European forest tree species' ranges.Our understanding of disturbances on species' ranges at this scale remains currently uncharted territory.
Considering this knowledge gap and the relevance of European forests, this paper characterizes changes in the distribution of key European forest tree species and the contribution of forest disturbance events to species' ranges for the period 2000-2020.Our research addressed two specific research questions: (a) are there observable range shifts in European forest tree species?(b) To what extent do forest disturbance events impact different tree species?To answer these we conducted a trend analysis on high spatial resolution (30 m) time series maps of six main European forest tree species realized distribution.We next aggregated the results of the trend analysis to coarser resolution (1 km) maps to capture the prevalent trends.Furthermore, we used available high resolution forest disturbance maps detailing both the year and type of disturbance to quantify how much of each species' range was affected by disturbances and identify the dominant disturbance types.By following this methodology, our analysis provides spatial evidence of tree species distributions dynamics, uncovers potential shifting patterns and provides species-specific insights on disturbance regimes.

Tree species probability of occurrence datasets
We used the realized distribution maps produced by Bonannella et al. (2022) to quantify the probability of occurrence over the twenty years period of analysis.In their study, they used a combination of EO data and Machine Learning (ML) to produce a time series for the period 2000-2020 covering the whole European continent for 16 forest tree species at 30 m resolution.To account for the multiple conditions that affect the realized distribution, the predictor variables used in the model include well known static (i.e., not changing over time) variables such as long term climatic and topographical data, but also dynamic (i.e.time series) variables such as multispectral data coming from Landsat.A novelty of their approach was to use chorological maps of different tree species as well to account for mutualism/competition effects.The presence-absence point dataset used for model training was published separately and is available on Zenodo (https://doi.org/10.5281/zenodo.6516590): presence data comes from a harmonized and preprocessed version of existing tree species data from the Global Biodiversity Information Facility (GBIF), the EU-Forest project (Mauri et al., 2017) and the LUCAS survey (EUROSTAT, 2017); absence data comes instead exclusively from the LUCAS survey.The authors then tested and combined the predictions of several ML algorithms using an ensemble approach known as stacked generalization (Wolpert, 1992) to deliver the final predictions.The probability of occurrence (0-100%) is provided for each species individually, irrespective of co-occurrence, meaning that the sum of the probabilities of different species in the same pixel can exceed 100%.The 2000-2020 time period is split in time windows of different extent: 2000-2002, 2002-2006, 2006-2010, 2010-2014, 2014-2018 and 2018-2020, for a total of six observations in time for each species for the period of analysis.The reason of these uneven time windows is due to design decision (i.e. a time window large enough to see changes in realized distribution over time and short enough to have multiple points in the period 2000-2020) and data constraints and irregularities (i.e. the SLC failure from Landsat 7 (Zhang et al., 2007)).The data is available on Zenodo, with a different entry for each species (example for Abies alba Mill.at: https://doi.org/10.5281/zenodo.5818021),and in European Environmental Data Cube (Witjes et al., 2023) at https://ecodatacube.eu.For the purposes of this paper we focused our analysis on only six of the tree species from Bonannella et al. (2022): Abies alba Mill.(silver fir), Fagus sylvatica L. (European beech), Picea abies L. H. Karst.(Norway spruce), Pinus nigra J. F. Arnold (black pine), Pinus sylvestris L. (Scots pine) and Quercus robur L. (Common oak).The species were chosen because of their common occurrence (Barbati et al., 2014) in European forest, their economic importance (Hanewinkel et al., 2013) and the fact that they represent more than 70% of the growing stock volume in the region (FOREST EUROPE, 2020).

Forest disturbance datasets
We used a combination of two datasets to collect data on the year of disturbance (Senf and Seidl, 2021a) and type of disturbance (Senf and Seidl, 2021b).Both datasets have a pan-European extent, 30 m spatial resolution and cover the time period 1986-2020.However, they do not cover all European countries (i.e., Iceland is missing); our analysis was therefore limited to the countries covered by the disturbance datasets.Firstly, the forest disturbances map only lists the year of the greatest disturbance over the whole 1986-2020 time series, regardless of the disturbance agent.We used the latest version of this product, available on Zenodo at https://zenodo.org/record/7080016 (v. 1.1.4).A combination of reference data from Landsat satellite images of visually interpreted disturbed and non disturbed patches (available on Zenodo at: https://zenodo.org/doi/10.5281/zenodo.3561924)and outputs from the well known LandTrendr algorithm (Kennedy et al., 2010) is used to train a random forest classifier; the model classifies each pixel in non-forest, undisturbed forest or disturbed forest.Information on the year of disturbance is provided from the trend analysis of spectral bands and indices.Secondly, the forest disturbances type map attributes to each disturbed pixel of the forest disturbances map a code based on type of event detected, as described in Senf and Seidl (2021b).We used the latest version of this product, available on Zenodo at https://zenodo.org/record/8202241 (v.1.2).Information on forest disturbance agents was gathered from the FORWIND database for storms (Forzieri et al., 2020) and the European forest fire information system (EFFIS) for fire related disturbances (San-Miguel-Ayanz et al., 2012); areas classified as disturbed but with low probability for wind or fire disturbances were assigned to the "other" class.A set of ten predictors (four spectral, three spatial and three landscape-based ones) was used to model the probability of a disturbed pixel to be affected by storm, fire or neither through a random forest model.Probabilities were then converted into hard classes in the final version of the disturbance maps, with the following legend: "0" for no disturbances, "1" for other disturbances (i.e.logging, drought, biotic disturbances such as bark beetles etc), "2" for wind-related disturbances and "3" for fire-related disturbances.In both cases, maps are provided on a country scale: for our analysis we first downloaded all the countries available and then merged them in a single file with European-wide coverage matching the same spatial properties (extent, pixel size, coordinate reference system) used by Bonannella et al. (2022) to produce the tree species probability of occurrence maps.We disregarded disturbances outside the time period of the analysis (prior year 2000).

Trend analysis and forest disturbances attribution
To analyze the change in probability of occurrence over time, we fitted simple OLS regression models with the probability of occurrence as the dependent variable and time as the independent variable.We applied this procedure at pixel level and for all six species individually.We then calculated the t-test statistic to determine the presence of an increasing or decreasing trend according to the OLS model.By combining the regression slope (β) of the linear model and the p-value from the t-statistics, we assigned each pixel to one of three categories: The chosen threshold of ± 0.25 for the regression coefficient denotes a total change of 5% in the probability of occurrence over the 20 years period of the analysis, which is considered a substantial change.Pixels in the stable class have relatively constant values of probability of occurrence over the time period analyzed.Pixels in the positive class have an increase in probability of occurrence over the time of the analysis while the opposite is true for pixels in the negative class; stable pixels could also experience fluctuations in the probability of occurrence, alternating favorable and unfavorable conditions for the species, but ultimately the net effect of these events would balance out by the end of the time series.The focus of the analysis was on identifying the direction of change in the probability of occurrence for each pixel and species, thereby revealing overall patterns in species distribution changes rather than quantifying the magnitude of the trend.
To better identify and analyze these patterns, we prepared speciesspecific 1 km resolution layers displaying the proportion (0-100%) of positive, stable and negative pixels within each cell: this aggregation of the 30 m results at 1 km allowed us to filter out local mature stand changes or other local variations, distinguishing genuine signals from isolated ones or false positives.Each species was modeled independently in Bonannella et al. (2022), so model accuracy varies between different species; the aggregation at 1 km helps in removing false signals caused by this variability.In addition, in areas undergoing intense forest management, due to clearcuts or relatively aggressive thinning operations, 30 m pixels may display very diverse spatial patterns, with a pixel having a strong positive slope (β > 1) next to a pixel with strong negative slope (β < − 1) (Fig. 1).While from a trend analysis perspective those isolated pixels may provide information on, respectively, the species advancing or retreating, the spatial context is instrumental in discerning weak signals from true niche shifts.The aggregation at 1 km filters those weak signals out and ensures us that the identified patterns are robust, providing a clearer understanding of the ecological dynamics at play.
As a final step, following the trend analysis and the spatial quantification, we overlayed slope values of each species with the forest disturbances maps.This allowed us to differentiate between disturbed and non-disturbed pixels and, for the disturbed pixels only, assign the year of disturbance and disturbance type.

Changes in probability of occurrence
For each species we visually inspected a histogram of the slope values.The modes of all histograms corresponds to zero slope (i.e., stable), as shown in Fig. 2.However, for four of the six species investigated, the modal histogram bin is very different from its direct neighbors.For these four species, the zero slope bin holds 40% to 50% of the total number of pixels, with beech having the greatest proportion, closely followed by spruce.Between 50% to 70% of the pixels are distributed in the "− 0.1", "0", and "0.1" bins, indicating that the majority of the pixels is classified as stable and not having a clear increasing or decreasing trend.For the two pine species instead, the "0" bin holds less than 20% of the pixels and the frequency distributions are remarkably flatter, with the majority distributed over the bins left of the mode.
Fig. 3 depicts the spatial patterns of positive, negative or no trends (stable in the figure): the column representing pixels with stable probability of occurrence over time is the one that has a more evenly distributed presence throughout the study area.The mapped distributions (Fig. 3) confirm the general trends in Fig. 2 that despite a large fraction of stable situations, there are much more negative than positive trends.Furthermore, the pixels in the stable column that have high values form large clusters that cover the extent of one or more.
countries rather than being concentrated in specific localized regions.This suggests that the areas considered stable cover a large portion of the total realized distribution of each species compared to areas with majority of positive or negative trends, which exhibit more distinct and localized hot spots.While this is true for almost all the species, the spatial patterns for black pine match what was highlighted previously in the distribution of the slope values: when compared with the other species, black pine is the sole species that lacks hot spots for positive trends and has a much more dispersed distribution of pixels with prevalent negative trends than pixels with a prevalent stable probability of occurrence over time.Except for black pine, the hot spots for negative trends consistently exhibit a notable pattern characterized by their localization at the latitudinal edges of each species distribution.
Even though five of the six species share this pattern, there are some differences in the negative hot spots' spatial locations; with most of them being located in Central and Northern Europe.For the silver fir, beech and oak, the hot spots seem to coincide with the northeastern edge of their distribution in the study area.The proportion values vary between 25% and 60%.Spruce presents a consistent negative hot spot towards the northernmost edge of its distribution like the other species, but also exhibits additional hot spots in the southern edge and in other areas.Notably, as the range extends further north, these hot spots are mostly located in regions with relatively low (500-1000 m.a.s.l.) elevation for the species (i.e.Carpathians, Ardennes) or even at sea level (i.e.Poland coastline).Similar considerations can be made for Scots pine, where the negative hot spots are located in the north (Finland, Baltic states) or at the lowest elevations of places located in the southern edge of the range (Provence in France, hilly areas surrounding the Pyrenees, the Cantabrian and Leon mountains in Spain).Areas with prevalent negative trends exhibit in general moderate values (40-50%) but certain portions of the negative hot spots for spruce, Scots pine, and black pine stand out with very high values (95-100%).Notably, these three species show elevated values primarily in the southern regions such as Spain, France, Slovenia, and parts of Romania, while the northern regions' hot spots remain in the moderate range.
Regarding positive trends, no clear hot spots can be discerned: pixel values vary between very low (5-10%) to low (20-30%) proportions and there are only few exceptional cases where they reach moderate (40-50%) proportions.Most of the pixels with low to moderate values correspond to specific mountain ranges (Alps, Norwegian Alps and Carpathians).Other areas with pixels in that value range are more species-specific: the Dinaric Alps for beech and oak, Scotland and Ireland for spruce and Scots pine, with Scots pine having some additional concentrations of positive trends in Galicia and Normandy.

Effects of forest disturbances
In addition to analyzing the slope values distributions, we flagged pixels affected by disturbances and used color-coded bars in the histograms to have a visual representation of the disturbance status of the pixels (see Fig. 2).Pixels affected by disturbances fall into the "0", "− 0.1" and "− 0.2" histogram bins, with a prevalence of disturbed pixels in the "0" bin.The similarities stop there, since each slope distribution seems to behave differently and no consistent trend can be picked up.Beech and spruce have the greatest proportion of pixels in the "0" histogram bin.Both black pine and Scots pine have an equal proportion of disturbed pixels across all histogram bins, mostly left of the mode (until bin "− 0.9"), while on the right it stops much earlier (bin "0.3").
The distribution of disturbed pixels over the years shows some common patterns among all the species (see Fig. 4 top); at least 2% of the total range of all the species is affected by disturbances every year, with most species having peaks in the year 2000, 2005 and 2020, which corresponds to a sudden increase of wind-related disturbance events.It is especially evident for the silver fir, where the percentage of pixels disturbed by wind goes from a value of 10-15% every year, up to more than 50% in the year 2020.There is an increasing trend of disturbed pixels over time, with all the species having on average 2% more of their range affected by disturbance events by the end of the time series.The other class is the most prevalent across all the years, reaching peaks of 90% in some years, with the exception of years affected by intense wind events where the percentage drops to 40-50% for some species.Among the coniferous species, spruce and Scots pine are the ones most affected by disturbances, with 5% of the range on average every year and peaks around 6.5%; while the silver fir has a lower average score.Silver fir is the coniferous species with the highest percentage of its range affected over the time period of analysis (8%).Black pine is the only species which is consistently affected by fire events (10% of the disturbed pixels on average every year, with peaks around 30%), while the other coniferous species are affected by fire much more sporadically.Both broadleaved species follow almost the same pattern in both proportion of affected range and distribution of disturbance types, with beech slightly more affected by fire events than oak.

Species' range assessment
In this study we conducted a trend analysis of the time series of probability of occurrence for different forest tree species on a European scale to identify potential shifts in tree species distributions in the last two decades.Our results show that, on average, the probability of occurrence among the pixels for all the species is rather stable for significant areas in Europe (no significant positive or negative trend), suggesting a predominant status quo in the species distributions over the study area.However, despite this high stability, there are some ongoing changes in species distribution patterns shared across the species and for different regions.While among these subtle variations we have identified both positive and negative trends, the latter are more prevalent than the former and that Northern and Central Europe shows more negative hot spots than other parts.Furthermore, the various patterns among the species analyzed suggest a different response to the wide range of environmental stressors each species is subjected to.These results still need to be approached with caution as for some of the analyzed species from Bonannella et al. (2022), the model performance is not very high.This could lead to underestimation or overestimation of the species' realized distribution and incorrect conclusions about increasing or receding species ranges.While the authors provide model uncertainty maps together with the probability maps, their analysis did not involve plot data such as NFI plots with precise locations.They used a filtering procedure to achieve a high spatial confidence (i.e.enough to be used for 30 m resolution maps) for datasets whose locational precision is either difficult to assess (i.e.GBIF) or is degraded to coarse resolution by definition (EU-Forest).For this reason, the precision of the distribution maps is expected to be somewhere between the species level and forest type level.In the current study we aggregated the results from 30 m to 1 km to filter out false signals (i.e.areas with intense forest management, local mature stands changes etc.); this process may also have canceled out the effect of potentially misclassified areas in the distribution maps.
In our analysis, it is interesting to note that areas with the highest values (90-100%) of pixels with no increasing or decreasing trend in probability of occurrence across all the species roughly align with what is reported in literature as their native range (Bonannella et al., 2022;San-Miguel-Ayanz et al., 2016): the places where species have historically persisted and adapted over time despite external perturbations are associated with a high degree of stability.While it is true that the current state of European forests is far from being natural (Strona et al., 2016) and that they have been influenced by humans for centuries (Kaplan et al., 2009), the identification of these stable regions within the native range of the species highlights the importance of preserving and protecting these areas to safeguard the diversity of European forests.For some species such areas can be considered as ecological refugia (Ashcroft, 2010), potentially providing valuable insights for conservation strategies.For instance, areas with stable probability of occurrence for the silver fir are localized in the Alpine regions of Central Europe (Tinner et al., 2013), well known as the species native range; for beech and oak vast parts of temperate Europe are classified as stable, which are areas where the two species have thriven for centuries (Fang and Lechowicz, 2006;Jones, 1959); for spruce in the central and northeastern Alps, the Carpathians and Scandinavia (Aarrestad et al., 2014;Latałowa and van der Knaap, 2006) and for Scots pine the stable areas are mostly located in Scandinavia, with some clusters on the border between Germany and Poland (Eckenwalder, 2009;Krakau et al., 2013).
While this overall trend is in agreement with the results of Maes et al. (2023), who showed that forest conditions across Europe have been even improving from 2000 to 2018, the areas with a prevalence of negative trends can't be ignored.Multiple studies have shown that climate change, by the end of the century, will considerably alter the distribution of European tree species (Mauri et al., 2022) and will be especially threatening for coniferous species (Dyderski et al., 2017), so it becomes crucial to scrutinize any signs of possible range contractions to differentiate between localized phenomena and large-scale shifts.The areas with a prevalence of negative trends can represent hot spots where this phenomenon could spread or, alternatively, they could be isolated episodes caused by extremely particular events occurring in those locations.
Our results reveal that the hot spots of negative trends generally lie at the edges of each species latitudinal range, suggesting these areas might be more exposed to unfavorable factors (climate change impacts, natural disturbances or lack of thereof, habitat fragmentation and speciesspecific disturbances).These factors could push the species beyond their optimal conditions, leading to range contractions or reduced fitness.Examples of the effects of these stressors on a European scale are numerous: Neumann et al. (2017) localized hot spots of tree mortality in both Northern (Scandinavia, mainly Finland) and Southern (Cantabrian and Leon mountain ranges in Spain) Europe for the period 2000-2012 due to extreme climatic variability, hot spots that align with our results for some tree species, mainly spruce, black and Scots pine; Senf et al. (2018) later proved how tree mortality rate, mainly due to changes in climate and land-use, has doubled in the last decade compared to the previous two.
Quantifying climate change impacts on tree species and their range may be difficult, as each species may respond differently to these stressors (Fei et al., 2017).Previous studies on the subject have shown a variety of responses, the most common being upward (Feeley et al., 2011;Lenoir et al., 2008;Morin et al., 2008) and poleward (Berner and Goetz, 2022;Hanberry and Hansen, 2015) shifts due to temperature increases, but downward shifts are also not uncommon (Crimmins et al., 2011).While in some cases a slight increase in temperature outside the species optimal range has proven to be beneficial for tree growth due to the CO2 fertilization effect (Martínez-Vilalta et al., 2008), a steady increase has proven to be one of the most important drivers of tree mortality, especially in Mediterranean Europe (Archambeau et al., 2020).Because of their economical importance, most of the species analyzed in this study have been historically planted outside of their native range, as is the case for spruce and Scots pine.
It is no coincidence that several of the negative trends hot spots that also have the highest proportion (90-100%) values for these two species are found at low elevations in Central (i.e.Slovakia, Czech Republic, Poland) and Southern (all Spain except the Pyrenees) Europe that are outside of their native range, where the species have been naturalized through continuous forest management practices (Caudullo et al., 2017;San-Miguel-Ayanz et al., 2016).These results align with a long list of studies assessing either stress or decline for the two species in those areas (Aguadé et al., 2015;Kolář et al., 2017;Rigling et al., 2013;Sedmáková et al., 2019).In contrast, the negative trends hot spots in the moderate (40-50%) range found in the northernmost part of Scandinavia for the same two species align with the results of Maes et al. (2023), which indicate a general reduction in ecosystem productivity (NDVI) and soil organic carbon (SOC) in the area rather than a species-specific issue.Hot spots for Scots pine extend further in the northern part of Finland, while spruce stays stable in the area: both trends have been confirmed by a tree-ring analysis conducted by Mäkinen et al. (2022) on the data from the last Finnish NFI, which found a declining trend for Scots pine especially in sites with mineral soils and in peatlands.Hot spots for very stable species like beech in south-eastern Sweden are also well known: beech does well especially in the south-west, where water availability is high, and gets progressively worse moving north, as shown by Martinez del Castillo et al. (2022).
Our results indicate that the proportion of positive trends is notably smaller compared to stable and negative trends for all the tree species analyzed, reaching moderate (40-50%) values at best and without consistent hot spots.This was also expected: colonizing new suitable areas through upward or poleward shifts to compensate for range erosion is a difficult task for most of the European tree species due to their low dispersal ability (Mauri et al., 2022;Svenning and Skov, 2004).Evidence of latitudinal (i.e.poleward or southward) shifts in communities of adult trees are still rare (Hanberry and Hansen, 2015) and trees natural dispersion capabilities are deemed to be insufficient to face the future risks posed by climate change (Jump and Peñuelas, 2005;Zhu et al., 2012), with assisted migration being a possible solution to the issue (Mauri et al., 2023).Our results for the species analyzed in this study suggest a higher degree of stability in the species' ranges, with signs of localized negative trends in the last two decades.A combination of lean and crash phenomena, as described by Lenoir and Svenning (2015), characterizes the situation observed in the species analyzed.The presence of overall stable range edges and shifting optimum, alongside localized declines, within the existing range, suggests significant dispersal limitations and the early stages of a shifting process that lags behind climate change.This phenomenon is often referred to as climatic debt in the scientific literature (Dullinger et al., 2012;Tilman et al., 1994).As we see in our results, the species can persist within their current range and even their native range, as indicated by the areas with prevalent stable probability of occurrence, despite environmental conditions being sub-optimal, for instance, for sexual reproduction and expansion.However, the species may eventually confront this climatic debt through abrupt range shifts, potentially leading to widespread forest dieback, particularly when critical thresholds are crossed due to repeated extreme climatic events or climate-induced pest outbreaks (Bałazy, 2020;Menezes-Silva et al., 2019;Pietrzykowski and Woś, 2021).

Species disturbances patterns
Our results show a consistent increase over time in disturbances over the species realized distribution, which is in line with the reported intensification of disturbance regimes in the last two decades for European forests (Hlásny et al., 2021a;Patacca et al., 2023;Senf and Seidl, 2021b).While our study didn't focus on the magnitude of these disturbance events on each species, the slope distributions show that most of the disturbed pixels have values of β coefficient close to 0 for the time period analyzed.This could be interpreted as disturbances having a negligible effect on the probability of occurrence during this time period; however, it might also be a result of the majority of the disturbance events being concentrated at the beginning or the end of the time series.The OLS model is fitted over only six points in time, with each point capturing the average situation of 3-4 years, and the information on the year of.
disturbance is limited to the greatest disturbance only.Information on recurring events over the same area is therefore not available, even though the probability of occurrence could experience fluctuations due to other disturbance events; this limitation in the data could lead to a severe underestimation of the effect of disturbances on species probability of occurrence.Furthermore, as for the species' ranges assessment, some areas of Europe the distribution maps may depict forest types rather than forest species.The disturbance analysis was conducted at 30 m and is not affected by the 1 km aggregation we applied to the species range analysis; any potential misclassification in the distribution maps directly influences species-specific conclusions of the disturbance analysis, with once again possible under-or overestimation of disturbance events per species.This is an additional reason on why we decided to focus on qualitative trends only regarding disturbances and we advice caution when reporting these results for future research.
Results for both the affected range and disturbance types for spruce and Scots pine align with what in general is reported in literature: well known for their high susceptibility to wind-induced disturbances (Mitchell, 2013;Peltola et al., 1997), information on the effects on these species of the major wind storms (Kronauer, 2007(Kronauer, , 2000;;Seidl and Blennow, 2012;Valinger et al., 2014) or bark beetles outbreaks (Jaime et al., 2019;Mezei et al., 2017;Wermelinger et al., 2021) are also well documented.Despite this, our findings indicate that for other species, like silver fir or even the broadleaved species, the percentage of pixels affected by wind disturbances in the years of severe storms is higher (≥ 50%) than for spruce and Scots pine.While it has been proved that silver fir is more sensitive to wind than Scots pine and less than the spruce (Schmidt et al., 2010), broadleaved species, and in particular the beechoak group, are in general less sensitive to wind than coniferous species (Klaus et al., 2011;Seidl et al., 2017;Wohlgemuth et al., 2022).This discrepancy can be attributed to the significant impact of the other disturbance class, which includes logging, drought and bark beetles, for spruce and Scots pine: not only these two species have a relatively larger distribution range compared to other species, but also higher harvesting rates (FOREST EUROPE, 2020).Coniferous species have also been more affected by drought (Pardos et al., 2021) due to forest management practices (i.e.monospecific stands) than broadleaved species and by bark beetles (Hlásny et al., 2021b), with spruce in particular being affected by multiple severe Ips typographus (L.) outbreaks in the last decade (Patacca et al., 2023;Romeiro et al., 2022).A combination of all these factors may overshadow the prominence of wind-induced disturbances for spruce and Scots pine.This highlights once again the importance of having accurate and up to date spatiotemporal information on forest disturbance drivers.
Black pine stands out among the other species not only for keeping an overall stable trend in affected range but also for being affected by fire more frequently and intensely: compared to the other species analyzed in this study, black pine native range includes regions where wildfires have historically been a regular occurrence (Isajev et al., 2004).When compared to other Mediterranean pines (i.e., Pinus halepensis Mill.or Pinus pinaster Ait.) however, black pine lacks fire-adaptation mechanisms (Rodrigo et al., 2004); it can tolerate low intensity surface fires, but crown fires are usually either stand replacing or lead to the development of non-forest ecosystems (i.e.grasslands or shrublands) (Ganatsas, 2010;Lucas-Borja et al., 2021).This behavior is consistent with the negative trends observed for the species in Catalonia, which, according to the disturbance type map we used from Senf and Seidl (2021b), is also one of the two hot spots mostly affected by wildfires in Mediterranean Europe.There is a chance that the majority of the pixels in the area classified as undisturbed but with negative trends were previously affected by wildfires in the years not covered by our analysis.Large fires in the area are common and the last one before 2000 was in 1998, with more than 24,000 ha of forest cover completely burned (Rodrigo et al., 2004).While fire can explain a significant portion of the disturbances and its higher proportion compared to the other species is significant, it is still limited in comparison to the prevalence of the other class.As we know that black pine is also fairly harvested (FOREST EUROPE, 2020), there is a significant component of logging in the other class.However, due to the lack of clear distinction within this class, further interpretation on the effects of the other class on black pine is challenging, making it a limitation in our ability to clearly infer the individual impact of the different disturbance types within the other class.
Broadleaved species show distinct patterns in their response to disturbances compared to coniferous species.While they exhibit similarities with other conifers in terms of the percentage of affected range and the proportion of wind and other disturbances, the distribution of disturbed pixels for broadleaved species is concentrated either in the mode of the slope distribution or evenly spread to the right of the mode (reaching "1.2" -"1.3" bins for the oak), indicating their ability to cope with disturbances.Other studies have highlighted how disturbances seem to be of minor concern for broadleaved species, to the point of being advantageous to them: Moris et al. (2023) showed that beech has high resilience to fires due to its fast resprouting ability and Scherrer et al. (2022) observed how beech outlived and even replaced spruce and Scots pine in those areas affected by wind or other disturbances.Salvage logging, a common practice for windthrown trees (Lindenmayer et al., 2012), has also contributed in the past to conversion from coniferous to broadleaved stands through natural regeneration only (Götmark and Kiffer, 2014;Jonášová et al., 2010).However, recent studies have pointed out that despite their resilience, broadleaved species may face challenges under increased frequency and intensity of disturbances, especially during prolonged drought spells.Mathes et al. (2023) analyzed the effect of three consecutive years of drought (2018-2020) on beech dominated forests in Germany and their study revealed a significant reduction in radial growth, indicating that severe droughts can override the effects of competition on tree growth.Overall, while disturbances might pose a danger by themselves for coniferous species, an increase in their frequency and intensity could potentially act as a tipping point for broadleaved species.

Conclusion
In this study, we conducted a comprehensive trend analysis of the probability of occurrence for the realized distribution of different forest tree species across Europe to investigate potential shifts in their distributions over the last two decades.Our results demonstrate an overall trend of stability in the distributions for many regions in Europe, with stable regions largely aligning with the native ranges of each species.These stable areas hold ecological significance and should be prioritized for conservation efforts to safeguard the diversity of European forests.However, we also identified areas with a high proportion of negative trends mostly situated at the edges of each species latitudinal range and exhibiting, specifically located in Central and Northern Europe for the species investigated.These areas may serve as potential hot spots of range contractions or reduced fitness for the species due to exposure to unfavorable factors, so they should be monitored closely in the coming years.Lastly, we identified areas with localized negative trends that were more species-specific, indicating that the response patterns to environmental stressors vary among different species.While our results align with several related studies, the potential limitations of the distribution maps used in the current assessment calls for caution and also highlights the need for more precise distribution maps that include useful information provided by detailed field surveys such as National Forest Inventories.Our analysis on forest disturbances focused on quantifying affected range and most prominent disturbance drivers for each species.Disturbance patterns have shown an increase over time in species affected range for all species, with some disturbance agents being very species-specific.The other disturbance class, remains the most prevalent disturbance type for all the species: this highlights the difficulty in precisely discerning the individual impact of different disturbances due to the absence of information on each disturbance agent other than wind and fire.
In conclusion, our study provides insights into distribution trends and disturbance patterns for some of the most important tree species in European forests in the recent past.Notably, this research introduces a novel approach by making use of Species Distribution Models mostly based on Earth Observation data, which enables consistent assessments of the ongoing changes in tree species distributions.This methodology contributes to understanding the observed trends, particularly the emerging changes in certain areas, which appear to be more pronounced in the recent years.Furthermore, the observed trends emphasize the significance of species-specific responses to environmental stressors and the need for more accurate spatiotemporal information on forest disturbances to monitor and understand these responses in the long term.It is imperative to proactively plan for the future species composition of European forests, alongside with preserving what we already have or adopting assisted migration practices.Careful evaluation and consideration of which species to plant where are essential to strike a balance between fostering adaptability and safeguarding the integrity of natural ecosystems.Our results for the last two decades, however, underscore the urgency of taking prompt actions, especially considering the long time period forest will require to adapt to climate change.

Fig. 1 .
Fig. 1.Example of the 30 m workflow for a production forest in Southern Sweden.Production forests are areas intensively managed: this creates local variations in probability of occurrence over time.Results of the trend analysis would then reflect this through chessboard-like spatial patterns in the slope values.The aggregation of the results at 1 km resolution isthen necessary to correctly interpret the results of the analysis.

Fig. 2 .
Fig. 2. Distributions of slope pixels per species; given the significant difference in absolute number of pixels across the species analyzed, the values are normalized.The color of the stacked histogram bins represents the frequency of pixels affected or not by a disturbance event.

Fig. 3 .
Fig. 3. Proportions of pixels with, from left to right, increasing (positive class) trend in probability of occurrence, no trend (stable class), and decreasing (negative class) trend in probability of occurrence over 1 × 1 km blocks for all species.The dark gray background marks the countries included in the study area.Note the legend employs unequal interval widths to facilitate distinguishing varying proportions and comparing spatial patterns within or across species.

Fig. 4 .
Fig. 4. Species specific distribution of disturbances over time.Top: percentage of 30 m pixels affected by disturbances relative to the species total realized range.Bottom: disturbance types normalized to the total number of disturbed pixels.