The study area covers the entire Basilicata administrative region in southern Italy (41°00'N – 16°45'E; 36°54'N – 15°21'E) which extends over an altitudinal range of 0 to 2,267 m. Most of the region is characterized by the Lucanian Apennines, a limestone mountain range encompassing approximately 170,000 hectares. This area includes 64 NATURA 2000 sites of the 92/43/EEC “Habitat” and “Bird” Directives (NATURA 2000 Viewer), two National Parks and seven Regional Parks. The vegetation landscape is characterized by evergreen and thermophilic deciduous oak forests (Fraxino orni-Quercion ilicis Biondi et al. 2003 and Carpinion orientalis Horvat 1958) within the coastal and hilly bioclimatic belt, as well as mixed woods and mesophilic oak forests (Fraxino orni-Ostryon carpinifoliae Tomazic 1940 and Melittio-Quercion frainetto Barbero et al. in Bonin et Gamisans 1976) within the submontane and lower montane belt. Beech woods (Geranio striati-Fagion Gentile 1970) are the dominant forest type throughout the upper montane belt of the whole region, with the co-dominance of Abies alba only in restricted areas of the Lucanian Apennines (Ruoti, Laurenzana) and Pollino massif. Only in the latter does a Pinus heldreichii open forest belt occur in the subalpine belt (Di Pietro and Fascetti 2005; Di Pietro et al. 2010; Fascetti et al. 2010; Biondi et al. 2014; Mucina et al. 2016; Di Pietro et al. 2021).

According to the checklist of the Italian vascular Flora (Bartolucci et al. 2024) 11 native oak species occur in Italy, whereas according to Flora of Italy (Pignatti et al. 2017-2019) this number rises to 19 species. This discrepancy can be attributed to the divergent systematic interpretations provided by various authors particularly concerning certain critical taxa belonging to the collective group of downy oaks (subgen. Quercus, sect. Quercus). Some of these taxa (e.g. Q. amplifolia Guss., Q. congesta C. Presl., Q. leptobalana Guss., etc.) are endemic to southern Italy, whereas others (e.g. Q. virgiliana (Ten.) Ten., Q. dalechampii Ten.) while maintaining their locus classicus in southern Italy they are also considered in the national floras of other European countries (Jordanov 1936–1979; Matyas 1970; Christensen 1997; Kaplan et al. 2022). Unfortunately, the analytical keys currently available do not yet allow for unambiguous distinction between the typical Q. pubescens and its taxonomically critical “sister species” within the collective group of Q. pubescens (Di Pietro et al. 2016, 2020a, 2020b). In view of the considerable morphological variability exhibited by this collective group, and the overlap in leaf and acorn morphological traits, when comparing different taxa within it, it was deemed most appropriate to make a generic reference to Q. pubescens s.l. in this paper.

To evaluate time variations, satellite images from the period 2015–2022 were used based on the four Sentinel-2 bands at 10 m × 10 m resolution: Blue, Green, Red, and NIR. The data were updated at intervals ranging from 8 to 16 days, and were acquired from 1 March 2015 to 31 August 2022. This data forms part of the European Space Agency’s Copernicus monitoring system (ECMWF 2018; Conte et al. 2019). The Normalized Difference Vegetation Index (NDVI) was downloaded from TerraClimate’s online database (Climate Hazard Sub-matrix) based on (Huntington et al. 2017) and available through the Climate Engine (2022–2024) web service (https://app.climateengine.com/climateEngine).

The presence and distribution of stands affected by oak decline syndrome, and their coenological identification, were assessed over the years between summer 2015 and summer 2022 (see Fig. 2). This was achieved by combining vegetation response data (NDVI), visible satellite image analysis, and field sampling based on the phytosociological approach (Braun-Blanquet 1964). Due to pragmatic considerations, each stand encompasses an area of 225 m² (corresponding to a plot size of 15 m × 15 m), constituting a slight excess approximation compared to the 200 m² proposed by Chytrý and Otýpková (2003) as a conventional standard plot-size for phytosociological relevés within European deciduous woodlands. The stands affected by oak decline are forest areas that clearly show symptoms of i) crown damage by die back (i.e. the death of twigs and other branches from the tips, starting from the apex of the branch and continuing down to the base) (Fig. 1); ii) leaf changes (i.e. the leaves may turn from yellow-green to yellow-orange or fall prematurely in September or October); iii) branch dieback and canopy thinning (i.e. the crown thins due to the loss of leaves, twigs and branches); and iv) cracks in the bark (i.e. cracks may appear in the bark, from which a watery fluid may run down). Furthermore, sprouts may also develop on the main branches and stems (epicormic branches) (Manion 1981; Ragazzi et al. 2000; Thomas et al. 2002; Choat et al. 2012; Keča et al. 2016; Colangelo et al. 2017; Conte et al. 2019).

In the initial phase of the study, the oak decline syndrome was assessed using the Normalized Difference Vegetation Index (NDVI) (Tucker 1979; Pettorelli 2013; Huang et al. 2021). It was utilized to identify plots showing a decrease in greenery and to monitor the spatio-temporal evolution of oak decline within the study area. In a subsequent step, visible satellite images were analyzed to rectify any issues identified by the vegetation response provided by the NDVI. The NDVI data were obtained from the Copernicus Sentinel-2 program (Climate Engine 2023). Field sampling was carried out in areas recognized as being affected by decline, and a Decline Severity Index (DES) was assessed in accordance with the protocol developed by Conte et al. (2019). The DES scale ranges from a value of 0 to 6, indicating the degree of damage observed in the stand due to oak decline syndrome. This observation was also utilized to elucidate the variability exhibited by the other variables.

The data set included 22 (18 quantitative and 4 qualitative) variables, collected over the period 2015–2022 (Table 1). The climatic variables were selected based on their potential role in triggering oak decline, as highlighted in numerous publications on forest decline and ecological stresses (Manion 1981; Lindner et al. 2010; Marques et al. 2025). These publications considered these variables individually (see Costa and Rodriguez 2017) or in association with other factors (see Brasier and Scott 1994; Colantoni et al. 2015; Kasim et al. 2025). The selection of the variables also considered certain of their attributes, such as the possibility of free use, a large spatial and temporal resolution and the accuracy and reliability of the data sources. The slope aspect (S and E in Table 1) was the only variable whose values were recorded directly in the field, while all the other data were obtained from online data sets (Funk et al. 2015; Climate Engine 2022–2024; Rodell et al. 2004; Rodell and Beaudoing 2007; MASE 2024). A comprehensive account of the forest stands, including their code, geographic coordinates, altitude, aspect, and dominant tree species is provided in the Suppl. material 1 (Suppl. material 1: tables S3, S4).

The following variables collected from the online data set: rainfall of the first semester: PL; surface incoming shortwave flux: Rad; average of the minimum temperature of the meteorological winter months: TW; average of the maximum temperature of the meteorological summer months: TS; surface wind speed: Wind; rain anomaly index: RAI; root zone soil moisture: RZsm; heatwave amplitude, as daily mean (2 m) temperature on hottest day satisfying the heatwave criteria of at least three consecutive days above the 90^th percentile: HWh; heatwave number, the count of events satisfying the heatwave criteria of at least three consecutive days above the 90^th percentile: HWc; longest number of consecutive days satisfying the heatwave criteria of at least three consecutive days above the 90^th percentile: HWL; elevation from digital terrain model: ELE; topographic ruggedness index: TRI; topographic position index: TPI; nitrogen dioxide: NO2; sulfur dioxide: SO2; formaldehyde: HCHO) (Acker and Leptoukh 2007; Climate Engine 2023). The distribution of quantitative variables has been conducted within four major groups according to the following parameters: climatic, ecological stress, topographic, and atmospheric pollution indicators.

Finally, four qualitative variables, namely: the type of oak wood and the dominant species: WoDs; type of forest management: ForM; degree of naturalness: NATUR; Ecopedological category: EcP, were extrapolated from the Regional Forestry Map (Corona 2006) and the Ecopedological Map of Italy (Rusco et al. 2003). The complete descriptions of the 22 variables can be found in the Suppl. material 1 (Suppl. material 1: table S1).

A data matrix comprising values recorded between 2015 and 2022 (eight years) for each stand was subjected to statistical analysis (Suppl. material 1: table S2). This analysis was conducted using JASP software v. 0.19 (JASP Team 2024) and PAST software (Hammer et al. 2001).

The explanatory variables were transformed into a smaller set of principal components (PC), which account for most of the data variance. Principal component analysis (PCA) was performed on 18 explanatory variables using correlation matrix option in Past software (Hammer et al. 2001). In accordance with the scores expressed by the variables on the first three components, it was determined that only those with scores greater than 0.6 would be selected for further processing using descriptive and inferential statistics (Stevens 1992; Field 2005).

Descriptive statistics are a means of summarizing and graphically representing data, including box plots and summary tables. The central tendency (mean, median and mode) and dispersion (variance and standard deviation) of the data were assessed using descriptive statistics. The Shapiro-Wilk test was performed to ascertain whether the assumptions required for an inferential parametric test were satisfied. Based on the results of the Shapiro-Wilk test, the non-parametric Kruskal-Wallis test was employed in conjunction with a post-hoc Dunn’s Test to analyze the impact of DES on the number of years, for each designated variable.

A vulnerability map was created by combining environmental and climatic variables most strongly associated with oak decline. The identification of these variables was achieved through a multi-step process involving remote sensing, field data, and statistical modelling. Variables with scores greater than 0.6 in the first three PCA components were selected. Given the non-normal distribution of the data set, the median values were utilized as a metric of centrality and regarded as thresholds which were subsequently spatialized.

Statistical raster data (based on the median threshold values) were subjected to a process of normalization, with the aim of rendering them comparable to a range between 0 and 1. This process involved adjusting the values that lay above and below the threshold (Gemitzi et al. 2010). In this approach, all parameters were given equal weight. No estimation was made of the relative suitability of one site over another based on geostatistical scores. This approach enabled the standardization of the values associated with oak decline as they emerged from statistical analysis. To standardize the spatial resolution, it was necessary to scale the raster grid to a minimum horizontal pixel resolution of 20 m × 20 m. The NDVI was calculated using band 8A as opposed to band 8, on the basis that band 8A is more sensitive to changes in vegetation status. Given that band 8A has a resolution of 20 m × 20 m, it was necessary to scale band 4 to match this dimension. Consequently, an NDVI with a resolution of 20 m × 20 m was utilized to discern variations in forests with greater efficiency (Mandanici and Bitelli 2016; Zhang et al. 2017).

All operations were executed using the open-source software QGIS 3.16 (QGIS Development Team 2023).

Thirty stands affected by oak decline were identified within the Basilicata region between summer 2015 and summer 2022 based on NDVI satellite images and field surveys where different DES values were assigned (Fig. 2). The comprehensive description of the forest stands, including their codes, geographic coordinates, altitude, aspect, and dominant tree species is available in the Suppl. material 1 (Suppl. material 1: tables S3, S4).

According to the physiognomic characteristics of the studied forest types (Suppl. material 1: table S4, column WoDs), the dominant tree species in 60% stands affected by oak decline was Quercus pubescens Willd. This was followed by Q. cerris L. (30%), Q. frainetto Ten. (6.66%), and Q. ilex L. (3.33%).

Based on the naturalness class reported in the Regional Forestry Map (Suppl. material 1: table S4, Column NATUR), approximately 76.66% of the stands were found to be managed as standard coppice with short cutting cycles and were classified as included in the low naturalness category. The “young high forests” stands Rv11 and Cs16 and the “mature high forests” stands Ab10, Gg28 and Sc29 were classified as being in the medium naturalness category. Finally, the “young forests” stands Va07 and Tc27 were classified as “high naturalness category stands”.

The eco-pedological category “Apennine reliefs composed of undifferentiated tertiary sedimentary rocks and developed in Mediterranean-montane climate” (Suppl. material 1: table S4, column EcP), was found to be the most common substrate within the study area, accounting for 40% of the total analyzed stands.

Eleven, out of the thirty stands analyzed, fall within various protected areas, including two National Parks, seven Regional Parks and six NATURA 2000 sites, SPAs, or SACs (Suppl. material 1: table S4, column Protected area).

Oak decline-related damage to forest vegetation was observed during the specified time interval in 2017, 2019, 2020, 2021, and 2022 with a total of 57 events distributed unevenly across the 30 stands considered (Suppl. material 1: tables S5, S6). No occurrence of oak decline damage was recorded in the years 2015, 2016 or 2018. The highest number of oak decline episodes occurred in 2017, affecting all 30 stands. This was followed by 2020 (11 stands), 2019 (7 stands), 2021 (6 stands), and 2022 (3 stands). The calculation of the total number of oak decline episodes for each stand over the entire time interval (2015–2022) revealed that Pi05 (Quercus pubescens) experienced the most episodes (five), followed by Sm08 (Q. pubescens), Rv18 (Q. cerris), and Sc29 (Q. cerris) with four episodes each.

The data of 18 quantitative variables for 30 stands recorded over 8 years was combined to obtain a data matrix of 18 × 240 (columns × rows) (Suppl. material 1: table S2). PCA analysis of this matrix yielded an explanation of 49.26% of the total variance along the first three axes.

As illustrated in Figure 3, Component 1 of the PCA clearly differentiates between the observation years 2015, 2016, 2018, 2019, 2020, and 2021, all of which fall entirely or largely within the negative section of the diagram, and the observation years 2017 and 2022, which fall enterely within the positive section. Component 2 reveals the extent of the dispersion of the stands within the same year. As shown in Table 2, the variables most significantly correlated with Component 1 are primarily responsible for the spatial separation of the different years in the PCA diagram. It is evident that all these variables serve as indicators of climatic and ecological stress. Specifically, RAI, RZsm, and PL exhibit a negative relationship with Component 1, whereas Rad, TS, HWh, HWc, and HWL exhibit a positive one. Component 2 is more closely associated with topographical and atmospheric pollution indicators. The data indicates a positive correlation with TW, wind speed, E, NO2, and TPI. Nevertheless, Component 2 exhibits a negative correlation with SO2, S, ELE, and TRI. In general, Component 1 demonstrated a stronger relationship with inciting variables, while Component 2 showed a stronger relationship with predisposition variables (Manion and Lachance 1992; Brown et al. 2018). The biplot of PC1 and PC3 highlights the role of NO2 and TW (which are negatively correlated) in differentiating between stands along the third component. Conversely, a minimal variability was observed in the biplot PC2 and PC3 (Suppl. material 1: Sheet S11). The spatial separation of the ellipses relating to different years confirms the episodic nature of climatic and ecological stress. The PCA also reveals that precipitation and heat-related indices are the variables that most negatively affect oak decline patterns. According to the loadings presented in Table 2, eight variables (PL, Rad, TS, RAI, RZsm, HWh, HWC, and HWL) exhibited loadings that exceeded the statistically significant threshold of 0.6 (Stevens 1992; Field 2005). Consequently, these variables were subjected to further analysis using inferential statistics.

The Shapiro-Wilk test indicated that, except for TS, the variables were not normally distributed. Consequently, the median value was used instead of the mean value in all the descriptive statistical analyses (Suppl. material 1: table S7).

In terms of the trend in climatic variables, the average precipitation (PL) was recorded during the first semester of the 2015–2022 time series. At 243.2 mm, PL in 2020 was the lowest on record for this period. This was followed by 2022 (246.7 mm) and 2017 (255.1 mm) (Fig. 4). For all the other years in the time series, PL was greater than 300 mm (in 2019 and 2021) and greater than 400 mm (in 2015, 2016 and 2018) (Suppl. material 1: table S7). Rad values exceeded 280 W/m² in 2021 (286.3 W/m²) and 2017 (282.6 W/m²) (Fig. 4 and Suppl. material 1: table S7). In all other years, including in 2015, 2016, 2018, 2020, and 2022, Rad values were consistently below 280 W/m². The highest TS value was recorded in 2017 (37.4 °C), followed by 2015 (35.2 °C). In all the other years, the temperatures were below 35 °C (Fig. 4).

Regarding the trend of ecological stress variables, the RAI was found to be negative for the years 2022 (-2.1), 2020 (-1.5), 2017 (-1.3) and 2021 (-0.2). In contrast, positive values were recorded in 2015 (0.8), 2016 (0.9), 2018 (0.5), 2019 (0.3) and, most notably, in 2023 (2.5) (Fig. 5 and Suppl. material 1: table S7). From 2017 to 2022, RZsm showed negative values, with the following measurements: 2017 (-41.1 mm), 2020 (-26.5 mm), 2022 (-25.9 mm), 2021 (-9.5 mm), 2018 (-5.3 mm) and 2019 (3.9 mm). The only years to record positive values were 2015 (2.752 mm), 2016 (11.7 mm), and 2023 (46.1 mm) (Suppl. material 1: table S7). HWc showed its highest median value in 2022 with 9 events, followed by 6 events in 2021, 4 events in 2017 and 2019, 3 events in 2015 and one event in 2018 and 2020. No heat wave events were registered in 2016 (Fig. 5).

The median value recorded for HWh in 2017 was the highest in the series, at 31.7 °C. The lowest medians were recorded in 2018 and 2020 (24.3 and 24.4 °C respectively). A clear increase occurred between 2020 and 2021, rising from 24.4 to 30.3 °C. The median values for 2019, 2021 and 2022 were all high (> 29 °C; 29.2 °C, 30.3 °C and 29.4 °C respectively), suggesting an upward stabilization in recent years following the decline between 2018 and 2020. The HWL for 2017 showed a median value of 8.5 days. This was followed by 2022 (6.8 days) and 2019 (6.0 days), 2021 (5.2 days) and 2020 (4.6 days) (Fig. 5). No HWL was recorded in 2016, whereas it was 5.6 days in 2015 and 3.0 days in 2018 (Suppl. material 1: table S7).

The Kruskal-Wallis test, performed using either the years of observation or the DES as discriminatory factors, revealed non-causal variability in the dataset (Suppl. material 1: table S8), except for the comparison between HWc and DES, where the test was not significant. The results of the post-hoc test (Dunn’s test) demonstrate the disparities between the various years from 2015 to 2022 as determined by the eight variables selected by PCA.

PL underwent significant changes over the eight-year period, with the most profound changes occurring after 2017, particularly between 2020 and 2022. This suggests a temporal threshold effect, meaning a key event or trend that began around 2017 and significantly altered precipitation amounts, particularly in 2020 and 2022.

Rad shows robust and sustained change over time, with nearly all comparisons between pairs of years being statistically significant. Significant differences were observed between 2015 and 2016, 2017, 2018, 2019, and 2021. In contrast, the difference between 2015, 2020, and 2022, is not significant.

TS showed significant variation between 2015 and 2016, as well as between 2015 and 2018. 2019 marked a turning point, with no significant differences observed compared to subsequent years. From 2019 onwards, the results suggest that the TS has reached a state of equilibrium.

Significant variations in RAI occurred between 2015 and 2017, 2020, 2021, and 2022, respectively. The last three years of the series, 2020, 2021, and 2022, showed some convergence in rainfall anomalies.

RZsm has changed significantly over time, particularly during years when oak decline occurred in the different stands (2017, 2020, 2021, and 2022, except 2019).

HWh shows a significant difference between 2015, 2017, and 2021. Data for 2016 are unavailable because no heat waves were recorded in that year.

The analysis reveals a significant increase in HWc between 2015 and 2016, 2018, 2020, 2021, and 2022. This trend suggests an intensification of extreme weather events in the final years of the series.

HWL changed significantly in the series between 2015 and 2018. The years 2019, 2020, and 2021 exhibit more similar trends than almost all other years.

Three vulnerability classes were identified for all the oak forests occurring in the study area based on the values exhibited by the most active variables in triggering oak decline. The median values were regarded as critical thresholds beyond which oak decline was deemed more probable (Suppl. material 1: table S9). The normalized raster layers were then combined by calculating the mean score for each spatial unit across all eight variables. This process yielded a composite vulnerability index for the forest areas (i.e., a single numerical index with different values for each pixel of the raster). The final value of the index, which ranged from 0 to 1.833 in the study area, was divided into three discrete values of ranges by statistical algorithms (from 0 to 0.701; from 0.702 to 1.430 and from 1.431 to 1.833). The aforementioned factors gave rise to three distinct vulnerability classes: low vulnerability (below 0.701), medium vulnerability (ranging between 0.702 and 1.430) and high vulnerability (above 1.430 up to 1.833). The resulting map is a raster file linked to a weighted linear combination of the theoretical binary variables. As illustrated in Table 3, approximately 58% of the vulnerable forest areas are classified as being “medium” class, 38.71% into the “high” class and 2.93% into the “low” class. A total of 31.7% of the oak forests in the Basilicata region (approximately 116,000 hectares out of 355,409 hectares of forest) have been classified as medium or high vulnerability (Fig. 6).

A total of 57 episodes of oak decline has been documented over the past eight years. The highest number of episodes was observed in 2017, followed by 2019, 2020, 2021, and 2022. The majority of oak decline episodes within the study area were found to be concentrated in the hilly bioclimatic belt, and, to a lesser extent, in the sub-montane belt, i.e. within an altitudinal range between 400 and 800 m a.s.l. This range is characterized by a preponderance of thermo-mesophilic and xero-thermic oak woodlands. These results corroborate the findings of previous studies carried out in the Italian Peninsula (Sicoli et al. 1992; Ragazzi et al. 2000; Gentilesca et al. 2015; Gentilesca et al. 2017; Conte et al. 2019; Coluzzi et al. 2020).

Throughout southern Italy, the hilly bioclimatic belt is characterized by the co-dominance of evergreen sclerophyllous oak forests (which are predominantly represented by Quercus ilex) and thermophilic deciduous oak forests, which are characterized by a significant presence of Quercus pubescens s.l. Q. frainetto, and, to a lesser extent, Q. cerris. In fact, the Q. cerris forests become dominant in the submontane belt and in the first fringe of the lower montane belt (see Bonin and Gamisans 1976; Aita et al. 1977; Zanotti et al. 1995; Di Pietro et al. 2014).

About the distribution of the oak decline episodes within the investigated forest vegetation, it was found that 96% of affected populations were in deciduous woods.

It is hypothesized that an increase in summer aridity, average temperatures and the number of heatwaves will have a greater impact on deciduous forests than on sclerophyllous evergreen forests. Nevertheless, the low incidence of oak decline in Quercus ilex evergreen forests is noteworthy, particularly since such forests have been severely impacted by oak decline due to climate change and pathogen outbreaks (Frisullo et al. 2018; Encinas-Valero et al. 2022) in other areas of Italy and the Mediterranean basin (e.g. Spain and North Africa). It is evident that even a forest type that is notoriously well adapted to the Mediterranean climate is increasingly affected by damage caused by global warming. This is an unequivocal sign of the extent of the ongoing climate change. In this context, the unexpected emergence of Q. pubescens forests as the deciduous oak forest type characterized by the highest number of oak decline episodes in our study area could also be interpreted as a slightly alarming finding. In fact, it is widely accepted that Q. pubescens is one of the most xerothermic European deciduous oak species. It is known to be well adapted to the Mediterranean climate, showing a broad ecological spectrum, as well as drought tolerance strategies and effective protection against high solar radiation (Damesin and Rambal 1995; De Rigo et al. 2016; Tognetti et al. 2019).

However, according to some authors (Brullo et al. 1999; Bacchetta et al. 2004; Biondi et al. 2004) the “typical” Quercus pubescens Willd is best suited to the continental or steppic woodlands in Central and Eastern Europe. It would also occur in northern Italy and in some relict areas throghout the Apennines. Other pubescent oak species, such as Q. amplifolia, Q. congesta, Q. dalechampii, Q. leptobalana, and Q. virgiliana, would largely substitute it in the warmer Mediterranean areas of the Italian peninsula and its major islands. Based on this theory, we could hypothesize that most of the pubescent oak forest stands that we investigated were, in fact, relict communities dominated by the typical Q. pubescens with, eventually, only a minimal participation from the pubescent oak “sister” species from southern Italy mentioned above. From this perspective, i.e. considering Q. pubescens as a relict species in the Mediterranean context that is presumably unprepared for the xeric climate change, the high incidence of oak decline that we detected among the Q. pubescens stands would find a plausible explanation. However, several biosystematics studies published in the last decade and aimed at establishing the taxonomic consistency of all the southern Italian taxa currently included in the Quercus pubescens complex seem to exclude the possibility of more than one pubescent oak species being present in the study area and in Italy as a whole (Di Pietro et al. 2016, 2020a, 2020b, 2020c; Proietti et al. 2021; Fortini et al. 2022).

Conversely, the hypothesis that the Q. pubescens forests would be the most damaged, given that they were the most widespread forest type in the study area is unfounded. This is because Q. cerris is the dominant species in most oak forests in the Basilicata Region (Di Pietro et al. 2010; Fascetti et al. 2010; Borghetti et al. 2024). We would have expected the latter to be particularly susceptible to oak decline. However, despite being much more demanding in terms of soil moisture than Q. pubescens (Gellini and Grossoni 1997), Q. cerris appears (at least in this case) to well tolerate the negative effects of global warming. One possible explanation for the ability of Q. cerris to adapt to climate change in the Basilicata region is simply topographical: Q. cerris forests in southern Italy are mainly found in the submontane and lower montane bioclimatic belts (600–1200 m) and especially on north-facing slopes (Aita et al. 1977; Blasi et al. 2018). Within these altitudinal range and northern slopes preference, Q. cerris forests receive greater amounts of precipitation, air humidity, decreased values of direct insolation (consequently decreased rates of evapo-transpiration) than Q. pubescenss.l. forests in the hilly belt, where the latter species ecological optimum is found. The location of Q. cerris forests at this “high” altitude enables them to mitigate the negative effects of summer droughts, heatwaves and other climatic anomalies to some extent.

A more detailed analysis of the role played by the individual variables considered in this study, revealed that several of these were only slightly significant in triggering oak decline. The topographic variables (ELE, TRI, TPI, E, S) were found to have a negligible effect, accounting for a mere 12.4% of the expressed variability, as demonstrated in Component 2 of PCA. As Manion (1981) previously observed, these variables function as “predisposition variables”. They amplify the destabilizing effect of the anomalous values expressed by the climatic variables, thereby rendering forest stands more susceptible to oak decline. The most logical explanation for the absence of significance of these variables is that the majority of the oak decline episodes occurred within the same altitudinal belt (the hilly belt) and under similar environmental conditions. A much larger sampling area and/or a significantly greater number of sampling sites might have provided additional data for finer-scale comparisons.

Variables linked to atmospheric pollution (i.e., NO2, SO2, HCHO) were found to be scarcely significant too. It is likely that the spatialization of data produced by a limited number of isolated air quality control units might not be adequate to record variations in the response chemical parameters on a finer scale (see Brown et al. 2018). This is particularly evident in the context of sampling sites that exhibit significant natural features, where it is presumed that disparities in the number and incidence of pollutants between contaminated and uncontaminated sites are minimal.

However, statistically significant variables were identified as triggers for oak decline. These were primarily climatic (PL, Rad, and TS) and ecological variables (RAI, RZsm, HWh, HWc, and HWL). As evidenced in previous studies (Brown et al. 2018; Colangelo et al. 2018; Senf et al. 2020; Bose et al. 2021; Bussotti et al. 2023; Gosling et al. 2024), this study found that variables related to drought and temperature anomalies (e.g. warm summers and high annual and seasonal drought) were the most significant in triggering oak decline episodes. In particular, the PCA revealed that the variables PL, Rad, TS, RAI, RZsm, HWh, HWc, and HWL displayed the highest correlation. The highest DES values were frequently detected in forest stands managed for closely spaced coppicing cycles (10–15 years) with a low number of standards left in place. In some cases, these standards were not reflected in the most mature trees (Fig. 7). The high values of TS, HWh, HWL and HWc (with extremes occurring in years when oak decline episodes were most evident and frequent) are consistent with this type of forest management.

As far as RZsm is concerned, it can be hypothesized that forest stands characterized by a short coppicing cycle (especially when coupled with a low number of remaining standards) may not adequately maintain sufficient humidity in the uppermost soil layers over extended periods.

The absence of a dominant tree layer and shrubby undergrowth renders the soil particularly susceptible to physical erosion and chemical leaching, resulting in a rapid decrease in litter thickness. The absence of continuous tree canopy and the scarcity of litter result in accelerated soil desiccation, slowed pedogenic processes and diminished nutrient availability. A forest that has been subjected to such a process is undoubtedly vulnerable to extreme weather conditions and, consequently, is more susceptible to oak decline. A salient finding of the study is the pivotal function of inappropriate land use and forest management. This finding aligns with the conclusions of earlier studies (Rambal and Debussche 1995; Hernández-Lambraño et al. 2019; Tognetti et al. 2019). It is estimated that approximately 76% of the forest stands affected by oak decline as identified in this study have undergone standard coppice management (Fig. 7). This has had a significant effect in reducing the structural complexity of forest ecosystems. Conversely, only negligible negative effects were recorded in the so-called “undisturbed stands”, i.e. those characterized by a higher degree of naturalness (Suppl. material 1: table S4).

One of the objectives of this study was to ascertain the threshold values of each variable that trigger oak decline. Our results demonstrated that oak decline is not triggered by the intensity of a solitary negative factor, but rather by the cumulative and combined effect of multiple factors.

As previously mentioned, the years 2017 and 2022 are particularly pertinent for comprehending this concept. Despite the PCA indicating comparable ecological stresses experienced during these two years, a significant disparity was observed in the number of oak decline events and the extent of recorded damage (Table 4 and Suppl. material 1: table S10).

A comparison of the two years in question reveals that 2017 is particularly noteworthy due to the heightened intensity of specific variables (i.e., very high summer temperatures TS and solar radiation Rad). However, the RAI (precipitation anomalies) does not attain extreme values, in contrast to the 2022 phenomenon. Conversely, 2022 exhibited drier conditions in terms of precipitation (RAI) and the frequency of heat waves (HWc), albeit with slightly lower intensity thermal anomalies (TS, Rad, HWh). It can thus be concluded that the severity of the oak decline episodes that occurred in 2017 appears to be linked to a synergic combination of the following factors: i) high average summer temperature (TS), ii) intense solar radiation (Rad), iii) particularly low soil moisture (RZsm), and iv) intense temperatures and the duration of heat waves (HWh and HWL). In contrast, although there was an increase in the number of heat waves in 2022, the intensity of thermal stress (TS, Rad, and HWh) and soil moisture (RZsm) were both lower than their respective critical thresholds. This resulted in a lower overall impact and more localized episodes of oak decline. The year 2017 serves as a paradigmatic example of a scenario in which multiple critical variables simultaneously exceeded their threshold values, thereby amplifying each other’s effects. In contrast, there was no similar peak in summer thermal stress in 2022, resulting in less pronounced mutual amplification. This is particularly noteworthy given the very low precipitation levels in 2022. These findings confirm once again that oak decline is not attributable to the anomalies of a single variable, but rather, it is the result of more complex interactions among multiple variables.

Table 4 clearly shows that, out of the eight identified variables as triggering oak decline, only four of them (Rad, RZsm, HWh and HWL) recorded the most critical values in 2017. However, it should be noted that the other two variables (PL and TS) exhibited values close to the highest or lowest critical values recorded in other years within the time interval considered.

In essence, each variable has a negative effect with a proportional relationship to its absolute value, contributing to an increase in the negative effects expressed by the other critical variables in the same proportion. In addition to the incremental contribution of a single variable to the effects caused by the others, an opposite trend was also identified: the “compensatory effect” between variables (Table 5).

The existence of this compensatory effect can be demonstrated through an examination of the trend in the values expressed by the variables PL, RAI, RZsm, and TS, Rad, HWh, HWL, and HWc in 2015, in comparison with other years. Despite the presence of anomalously high temperatures throughout the summer months (Table 4), 2015 does not appear to be among the years characterized by oak decline, for any of the forest stands considered. In this case, the abundant rainfall experienced in the first semester of the year served to counterbalance the deleterious effects of the unusually high summer temperatures.

A vulnerability map which considers three oak decline risk categories (low, medium and high), was created for the entire Basilicata region (Fig. 6). It summarizes all the concepts and assessments that have been previously expressed. The situation appears quite alarming. A large proportion of the forest heritage in the study area is classified as medium or high risk. Unfortunately, the plots in the high-vulnerability category far outnumber those in the low-vulnerability category. Nevertheless, caution is required when drawing conclusions or hypothesizing future scenarios, as this constitutes merely the first proposal of its nature. As there are no vulnerability maps available for previous years, it is not possible to determine whether the forest situation of oak decline in the study area is worsening, remaining stable, or improving over time. Nevertheless, we can be certain of excluding at least the last one of these three hypothetical options. The constant increase in global warming over the last two decades, together with the ever-worsening desertification of the Mediterranean region (Colantoni et al. 2015; Carvalho et al. 2022), does not seem to support such an optimistic hypothesis. However, according to some recent estimates, a promising future seems to be opening for the Quercus genus. Some authors (Hanewinkel et al. 2013; Holland et al. 2014; Früchtenicht et al. 2018) suggest that the Quercus group could expand its potential distribution range in Europe (at the expense of other forest types) from 11% to 28–40% by 2100. While we don’t want to dampen enthusiasm, our vulnerability map indicates that prospective episodes of oak decline, on a regional scale, could affect 63.43% of mesophilic and meso-thermophilic oak forests. Given that future climate change scenarios indicate an increase in heatwaves and drought periods (ESOTC 2024), it is more likely that a significant proportion of the land acquired by oak forests, as postulated in the aforementioned studies, will be heavily taxed by oak decline.