**4. Discussion**

The quality of the reference baseline climate has a fundamental role in predictions of the potential impact of climate change on organisms and natural ecosystems. The stability and reliability of the estimated projections calculated from species distribution models [14,45,46], managemen<sup>t</sup> simulators [18,47], or the estimation of the geographical shift of climate zones [48] rely on the differences between current and future climate. While the representativeness of WorldClim is adequate concerning air temperatures, large differences were found in the precipitation surfaces. Our results demonstrate a systematic difference of 0.76 ◦C between observed and interpolated values. According to the most recent IPCC report [49], the observed increase in temperature has been around 0.2 ◦C per decade. As a consequence such difference might affect future projections of WorldClim dataset adding uncertainties on the further modelling efforts [22,50]. In this case a likelihood analysis should be more adequate than deterministic ones and in order to include a sensitivity analysis and evaluate the probability of success of empirical models based on phenotypic plasticity and applied future projections [51]. Precipitation, by contrast, proved to be the main weakness of WorldClim surfaces in Europe. Despite its relatively small spatial extent, the European environment is characterized by many different forest systems that reflect broad climate variability, spanning from the Mediterranean to the Arctic.

The 1961–1990 baseline period is a fundamental dataset for ecological modelling because records from earlier periods were often affected by different instrumentation or changes in observational practice [30,37]. Therefore, numerous studies from climatology to biology, ecology and forestry [36,48,52] have used this baseline period, and WorldClim has been used extensively. We can expect a further warming trend in the next two decades at a rate of about 0.1 ◦C per decade, due mainly to the slow response of the oceans. As a consequence, even though the linear regression analysis showed a good match between observed and interpolated data (adjusted R<sup>2</sup> = 0.856), the difference is higher than the expected rate of change, which could heavily affect model predictions, adding uncertainties on future projections and smoothing results (i.e., land suitability projections) in an uncontrolled way [14,53–55]. This issue is then amplified when analysing MAP, where higher differences were found in combination with a poor regression analysis result. As a consequence, important biases may be introduced when using WorldClim's precipitation dataset. This is particularly true when WorldClim is used as the reference line and climate projections are locally downscaled and added to the WorldClim surfaces, as in the "Delta method" [56]. As a result, the calculation of climate indices might be difficult. For example, many studies used reference evapotranspiration [57–59] as the main predictor in statistical models [3,60,61]. In this case, the mathematical combination of differences in MAT and MAP might introduce uncontrolled biases through the study area. These biases could

represent a critical issue, especially in the Mediterranean and anywhere else that moisture deficit is identified as the most relevant climate driver.

The main advantage of our compiled dataset might be its representativeness at small scale. The authors of WorldClim themselves warn that the high resolution of the climate surfaces does not imply high data quality in all places as this depends on local climate variability, quality and density of observations and the degree of the fitted spline [31]. In a similar study, when compared with PRISM and Daymet datasets for the continental United States, many concerns were expressed, especially regarding the quality of WorldClim's precipitation grids in mountainous areas [34,35,62,63]. For this reason, several studies at regional or national scale at higher resolution (e.g., 100–250 m) preferred the use of meteorological variables obtained at nearby observational sites [28,64–66]. Regardless of the distances of the investigation sites from the locations where meteorological datasets were gathered, orography and land use, and the surrounding area and variable characteristics, must be considered. At small scale their variability may be a strong driver of frequently overlooked heterogeneities, leading to significant discrepancies in transferred datasets used for otherwise appropriate processing methods [67,68].

The lack of any relationship between BIAS and the main physiographic parameters (i.e., latitude, longitude, and elevation) does not allow for any statistical adjustment (e.g., downscaling, locally calibrated lapse rate, etc.) for either temperature or precipitation. However precipitation regimes are very di fficult for meteorological stations to record properly and this issue has often been found in other databases [59,69]. Many more data are required, especially in the case of forest monitoring, as a result of the lack of temporal autocorrelation during the timeframe [70].

The need for a freely available and representative global climate dataset is large and growing, as evidenced by WorldClim's citation statistics. These goals can be achieved with local up-to-date monitoring networks, which could play a key role in evaluating global grids at small scale [71] as well as providing data for the construction of additional global climate datasets. Harmonization e fforts, as well as increased representativeness of the established networks, are paramount for construction of more accurate climate surfaces. Enhanced data recovery with regular spatial coverage may overcome the lack of dense environmental or climatological sampling [28,70,72,73]. Derived surfaces are fundamental in order to plan future managemen<sup>t</sup> strategies. For instance, and concerning forestry, additional strata, such as homogeneous climate zones, are needed as a fundamental tool to plan the transfer of genetic resources and reproductive materials across specific geographic areas [74,75]. WorldClim grids were interpolated with spline functions, a fast method known to yield results similar to polynomial functions but without mathematical instability. Such methods do not consider the spatial autocorrelation between observations, only partially achieved by more complex models where latitude and longitude are included as predictive variables [28,76]. Therefore, the exhibited spatial aggregation of the BIAS in the case of denser observations of our dataset (i.e., Sweden and Germany) may be relevant for research activities and improvements of the climate surfaces.
