Is New Always Better? Frontiers in Global Climate Datasets for Modeling Treeline Species in the Himalayas

Comparing and evaluating global climate datasets and their effect on model performance in regions with limited data availability has received little attention in ecological modeling studies so far. In this study, we aim at comparing the interpolated climate dataset Worldclim 1.4, which is the most widely used in ecological modeling studies, and the quasi-mechanistical downscaled climate dataset Chelsa, as well as their latest versions Worldclim 2.1 and Chelsa 1.2, with regard to their suitability for modeling studies. To evaluate the effect of these global climate datasets at the meso-scale, the ecological niche of Betula utilis in Nepal is modeled under current and future climate conditions. We underline differences regarding methodology and bias correction between Chelsa and Worldclim versions and highlight potential drawbacks for ecological models in remote high mountain regions. Regarding model performance and prediction plausibility under current climatic conditions, Chelsa-based models significantly outperformed Worldclim-based models, however, the latest version of Chelsa contains partially inherent distorted precipitation amounts. This study emphasizes that unmindful usage of climate data may have severe consequences for modeling treeline species in high-altitude regions as well as for future projections, if based on flawed current model predictions. The results illustrate the inevitable need for interdisciplinary investigations and collaboration between climate scientists and ecologists to enhance climate-based ecological model quality at meso- to local-scales by accounting for local-scale physical features at high temporal and spatial resolution.


Introduction
The Himalayas are the largest mountain system in the world with a distinct threedimensional geoecological differentiation, with a high variation of altitude, temperature, precipitation, soils, and land use [1][2][3]. The small-scale heterogeneity of habitats and site conditions provides a multitude of ecological niches associated with high diversity of species and communities [4][5][6]. Compared to other regions of the world, anthropogenic influence is less pronounced, but even in the Himalayas, near-natural treeline sites can hardly be found [7]. In the course of past climate changes, refugia at higher elevations facilitated the preservation of biodiversity, and even greater importance could be attached to these refugia in the future [8].
In recent decades, interest in treeline research in the world's mountains has considerably increased. Research on treeline causation, spatial patterns, and treeline dynamics has particularly intensified from the 1990s to today [9]. Inspired by easy-to-use software and freely accessible data, analyzing underlying factors of treeline position and dynamics as well as forecasting future changes of treeline species composition, altered species phenology, and future species ranges represents an ongoing trend in recently published

Methodology of Climate Datasets
Worldclim 1.4 [31] The WORLDCLIM 1.4 climate dataset is based on interpolation of point observations from weather stations, which regionalizes monthly observations of precipitation and temperature with a thin-plate smoothing spline algorithm, using latitude, longitude, and elevation as statistical predictor variables. Thus, climatic variations from areas with high station densities are fairly well resolved, whereas areas with limited weather stations yield higher insecurity. To account for this uncertainty in input data, weather station density, elevation bias in the weather stations, and elevation variation within grid cells, data partitioning and cross-validation were applied. Despite the high spatial resolution of 30 arc seconds, WORLDCLIM 1.4 neglects local-scale atmospheric processes which are crucial for the formation of site-specific topo-climatic conditions in high mountain environments. The authors clearly emphasized that uncertainty in WORLDCLIM 1.4 climate estimates resulting from high variation within grid cells is highest in mountainous and poorly sampled regions. Irregular distribution and low density of weather stations is a common feature of almost all mountain areas outside of Europe and North America, and especially of the Himalayas. At the time of publication, WORLDCLIM 1.4 had higher spatial resolution, higher number of weather station records, improved elevation data, and more information on spatial distribution of uncertainties compared to other hitherto available products. For WORLDCLIM 1.4, monthly minimum, maximum, and mean temperature, monthly precipitation sums, and bioclimatic variables are freely available for download.
Worldclim 2.1 [32] For WORLDCLIM 2.1, weather station data were interpolated using likewise thin-plate splines with covariates including elevation, distance to the coast, and satellite-derived covariates (i.e., maximum and minimum land surface temperature and MODIS cloud cover). Data from weather stations (9000-60,000) between 1970 and 2000 included monthly tem- peratures (minimum, maximum, and mean), precipitation sums, incoming solar radiation, vapor pressure, and wind speed.
Compared to the former version, prediction accuracy for temperature variables increased by 5-15% (0.07-0.17 • C) due to remote sensing data. However, prediction error remained high for all climate variables in areas with low weather station density (i.e., remote, mountainous regions in less developed countries).
For the final product, the best performing model for each region and variable was selected. Finally, global cross-validation correlations revealed r ≥ 0.99 for temperature and humidity, 0.86 for precipitation, and 0.76 for wind speed. The authors highlight the necessity to expand weather station data networks, since remote sensing covariates only marginally improved most WORLDCLIM climate estimates in version 2.1. Besides monthly minimum, maximum, and mean temperature, monthly precipitation sums, and bioclimatic variables, WORLDCLIM 2.1 offers solar radiation, wind speed, water vapor pressure, and elevation. Chelsa 1.1 [33] CHELSA 1.1 climate data are based on a quasi-mechanistical statistical downscaling of climatic reanalysis interim (ERA-Interim) global circulation model with Global Precipitation Climatology Centre (GPCC) and Global Historical Climate Network (GHCN) bias correction. The algorithm for temperatures in CHELSA 1.1 relies on a statistical downscaling of monthly vertical atmospheric means of daily mean temperatures at ERA-Interim pressure levels from 1000 to 300 hPa. The free atmospheric temperature lapse rates were calculated using linear regression for each ERA-Interim grid cell. After the reduction to sea level, it was interpolated between grid cells and afterwards projected back on the elevational surface of the DEM GMTED2010 [47]. For the precipitation algorithm, orographic predictors such as wind fields, valley exposition, and boundary layer height with a subsequent bias correction were incorporated. The latter includes bias correction at different scales (large scale at~150 × 150 km and small scale at~0.75 • × 0.75 • latitude/longitude) and temporal aspects (backward bias correction) for aggregation of monthly precipitation sums to annual rainfall estimates. The resulting data consist of a monthly temperature and precipitation climatology for the years 1979-2013 available at a resolution of 30 arc seconds. Especially, precipitation amounts showed higher accuracy compared to precipitation data from other climate datasets (e.g., WORLDCLIM). The authors emphasized the appropriate performance of CHELSA 1.1 in mountainous regions since orographic precipitation effects were considered. For CHELSA 1.1, monthly minimum, maximum, and mean temperature, precipitation, and bioclimatic derivates are available. Chelsa 1.2 [34] For the latest version of CHELSA climate data, the original quasi-mechanistical statistical downscaling approach was enhanced to compute more accurate temperature and precipitation estimates.
As in version 1.1, free atmospheric temperature lapse rates (1000 to 300 hPa) were used to estimate near-surface air temperatures at sea level, which were then interpolated between grid cells using B-spline interpolation, and afterwards projected back on the elevational surface of the DEM GMTED2010 [47]. Maximum and minimum temperatures were calculated using climatological-aided interpolation. Mean monthly temperature values were used and added or subtracted to the maximum or minimum daily temperature derived from the three-hourly data of previous post-processing data fields in ERA-Interim of minimum or maximum temperature.
The downscaling of precipitation is in accordance with version 1.1, applying the bias correction with GHCN and GPCC data. Advances in the latest version refer to the calculation of annual precipitation sums for the climatologies and the handling of bias correction, which now includes GPCC data at a higher resolution of 0.25 • × 0.25 • latitude/longitude.
The authors aim at providing climatologies for precipitation and temperature at high spatial resolution for environmental and ecological studies, which might prove valuable in varied scientific applications that rely on a good representation of small-scale precipitation and temperature patterns. For CHELSA 1.2, the following additional predictors are available: aridity, hydrothermic coefficient, potential evapotranspiration, growing season length, temperature and precipitation, growing degree days, frost change frequency, glacier height, and solar radiation.

Modeling Procedure
To evaluate the effect of different climate datasets when modeling the ecological niche of B. utilis in Nepal under current and future climate conditions, we followed the approach described by Bobrowski et al. [36]. Out of 360 occurrences of B. utilis, 13 georeferenced records  were accessed via the Global Biodiversity Information Facility (gbif.org), 31 records were extracted from a database compiled from a literature survey ( [15]; unpublished data), and 316 occurrences were obtained via satellite images of Google Earth (GoogleEarthTM, ver. 7.1.1.1888, Google, 2015) ( Figure 1). These occurrence localities were validated through expert knowledge, obtained from numerous field visits in the Himalayan region. This method has been proven to be useful in global treeline research [48,49].
The downscaling of precipitation is in accordance with version 1.1, applying the bias correction with GHCN and GPCC data. Advances in the latest version refer to the calculation of annual precipitation sums for the climatologies and the handling of bias correction, which now includes GPCC data at a higher resolution of 0.25° × 0.25° latitude/longitude.
The authors aim at providing climatologies for precipitation and temperature at high spatial resolution for environmental and ecological studies, which might prove valuable in varied scientific applications that rely on a good representation of small-scale precipitation and temperature patterns. For CHELSA 1.2, the following additional predictors are available: aridity, hydrothermic coefficient, potential evapotranspiration, growing season length, temperature and precipitation, growing degree days, frost change frequency, glacier height, and solar radiation.

Modeling Procedure
To evaluate the effect of different climate datasets when modeling the ecological niche of B. utilis in Nepal under current and future climate conditions, we followed the approach described by Bobrowski et al. [36]. Out of 360 occurrences of B. utilis, 13 georeferenced records  were accessed via the Global Biodiversity Information Facility (gbif.org), 31 records were extracted from a database compiled from a literature survey ( [15]; unpublished data), and 316 occurrences were obtained via satellite images of Google Earth (GoogleEarthTM, ver. 7.1.1.1888, Google, 2015) ( Figure 1). These occurrence localities were validated through expert knowledge, obtained from numerous field visits in the Himalayan region. This method has been proven to be useful in global treeline research [48,49]. To limit the space for pseudo-absence selection (i.e., the space where the species is known not to occur, e.g., the Tibetan Plateau and the Indian lowlands), a convex hull with a buffer of 5° was drawn around the occurrences utilizing the R-package 'rgeos' [50]. Additionally, a buffer of 5 km was placed around the known occurrences to ensure that To limit the space for pseudo-absence selection (i.e., the space where the species is known not to occur, e.g., the Tibetan Plateau and the Indian lowlands), a convex hull with a buffer of 5 • was drawn around the occurrences utilizing the R-package 'rgeos' [50]. Additionally, a buffer of 5 km was placed around the known occurrences to ensure that pseudo-absences where not placed directly next to occurrences. Following the procedure described in [25], a total of 4000 pseudo-absences were generated as randomly stratified points (i.e., one pseudo-absence per grid cell), as implemented in the R-package 'sp' [51].
The data were split into 80% training and 20% testing data utilizing the 'caret' package [52]. We applied generalized linear models (GLM, [53]) and followed the procedure described in Bobrowski et al. [36], with binomial distribution, logit-link function, and polynomial terms of second order were used to build the model. We calculated the root mean square error (RMSE) in order to account for overfitting of the data (i.e., RMSE should be very similar between training and testing datasets), indicating a good fit of all models. Five model runs were averaged, and final models were evaluated using Pseudo R 2 [54] to determine the explained variance of the test data.
To compare the different CHELSA and WORLDCLIM climate datasets (Table 1), we downloaded the bioclimatic variables of the respective datasets directly in R using the package 'ClimDatDownload' [55], which is a recently published package designed to download, process, and compare different climate datasets.
In accordance with the species-specific ecology and to allow for a direct comparison, the same climatic variables were used (after Spearman's rank correlation coefficient of r s ≤ |0.7| with the 'psych' package [56], which was applied to check for autocorrelation between predictor variables) viz. Mean Temperature of Wettest Quarter (i.e., mean temperature of the growing season), and Precipitation of the Coldest Quarter (i.e., proxy for melt water availability in the pre-monsoon season). Additionally, we calculated the average mean precipitation sums of March, April, and May, and included that variable into the models to account for pre-monsoon drought stress [7,57]. All utilized bioclimatic variables have been shown to reflect ecological requirements of B. utilis at treeline sites, accounting for growth limitations due to low temperature and insufficient moisture availability, and to describe the ecological niche (for detailed interpretation, see [35,36]).
To demonstrate the significance of current predictions, we projected the ecological niche of B. utilis under future climate change conditions. Only for CHELSA 1.2 and WORLDCLIM 1.4 future scenarios were available at 30 arc seconds resolution and subsequently downloaded, cropped, and snapped directly in R with the package 'ClimDatDownload' [55]. We used future scenarios of the American National Center for Atmospheric Research (NCAR) with their Community Climate System Model (CCSM4, [58]) with two RCPs, namely 4.5 as medium-low and 8.5 as high scenario for the timeframe 2080. In detail, RCP 4.5 assumes an increase in temperature of around 2 • C relative to the period between 1950 and 2000, whereas in RCP 8.5, an increase of around 4 • C is predicted for most parts of the study area. Furthermore, a slight increase in precipitation amounts is projected for central Nepal, which is highest under RCP 8.5.
All statistical analyses were performed using the programming language R [59] and all maps were created in ArcMap 10.7.1 [60].

Comparison of CHELSA and WORLDCLIM Climate Datasets
Prior to the modeling approach, we inspected selected predictor variables and found striking differences between the climate datasets and respective versions. As mentioned above, differences between temperature-related variables were smaller compared to precipitation-related variables.
For the Mean Temperature of the Wettest Quarter, all climate datasets show highly significant positive correlations (Table 2), and similar value ranges with consistent distribution ( Figure 2). The negative correlation between increasing altitude and decreasing temperature becomes clearly visible. WORLDCLIM 2.1 reproduces comparatively low temperature values on the Tibetan plateau ( Figure 2).  The comparison of precipitation-related variables between the climate datasets revealed major differences, and even among versions, wide disparities could be found (Figures 3 and 4). Correlation analyses showed high coefficients among versions and lower coefficients between climate datasets, whereas coefficients were higher for CHELSA versions and WORLDCLIM 2.1 compared to 1.4 (Table 3). These differences are mainly attributed to different origin of input data and calculation methods between the versions. Precipitation-related variables of CHELSA versions are based on precipitation data from ERA-Interim data, which are downscaled and include a three-fold bias correction [34]. CHELSA precipitation estimates are corrected using GPCC and GHCN data to execute monthly bias correction, inclusion of orographic effects, and weather station bias correction (for more details, see [34]). WORLDCLIM-related variables are solely based on station data, which are point data that are interpolated with latitude, longitude, and elevation as independent variables [31]. The very low density of weather stations on the Tibetan plateau is one of the main challenges, leading to an unequal distribution of measurements for the climate dataset calculations and to the need to apply different correction methods. The comparison of precipitation-related variables between the climate datasets revealed major differences, and even among versions, wide disparities could be found (Figures 3 and 4). Correlation analyses showed high coefficients among versions and lower coefficients between climate datasets, whereas coefficients were higher for CHELSA versions and WORLDCLIM 2.1 compared to 1.4 (Table 3). These differences are mainly attributed to different origin of input data and calculation methods between the versions. Precipitation-related variables of CHELSA versions are based on precipitation data from ERA-Interim data, which are downscaled and include a three-fold bias correction [34]. CHELSA precipitation estimates are corrected using GPCC and GHCN data to execute monthly bias correction, inclusion of orographic effects, and weather station bias correction (for more details, see [34]). WORLDCLIM-related variables are solely based on station data, which are point data that are interpolated with latitude, longitude, and elevation as independent variables [31]. The very low density of weather stations on the Tibetan plateau is one of the main challenges, leading to an unequal distribution of measurements for the climate dataset calculations and to the need to apply different correction methods. Although more consistent between versions than WORLDCLIM, CHELSA versions show remarkable differences. One major difference is the bubble-shaped structure of the northern edge of the values in version 1.2 compared to a smoother curve shape in version 1.1 (Figures 3 and 4). This difference is attributed to weighting and smoothing factors applied to predictions when correcting precipitation values by including orographic effects. In version 1.1, the monthly bias-corrected precipitation grids are proportionally distributed at ERA resolution and boundary layer of corrected wind effect surface, whereas in version 1.2, this relationship is assumed to be linear to achieve data scaled to 0.75 • resolution [33,34]. For both approaches, precipitation values are erroneously overpredicted in topographically caused rain-shadow regions behind the Himalayan main range and on the Tibetan Plateau.
predictions when correcting precipitation values by including orographic effects. In ve sion 1.1, the monthly bias-corrected precipitation grids are proportionally distributed ERA resolution and boundary layer of corrected wind effect surface, whereas in versio 1.2, this relationship is assumed to be linear to achieve data scaled to 0.75° resolutio [33,34]. For both approaches, precipitation values are erroneously overpredicted in top graphically caused rain-shadow regions behind the Himalayan main range and on th Tibetan Plateau.    The two WORLDCLIM versions show severe artefacts on the Tibetan plateau, wi these artefacts appearing ball-shaped in version 1.4 and more flattened and stretching t wards the East in version 2.1 (Figure 3). Differences between the two versions can be pa ticularly attributed to the number of weather stations used for interpolations. Fick an Hijmans [32] stated higher prediction error with increase in station elevation and distan  The two WORLDCLIM versions show severe artefacts on the Tibetan plateau, with these artefacts appearing ball-shaped in version 1.4 and more flattened and stretching towards the East in version 2.1 (Figure 3). Differences between the two versions can be particularly attributed to the number of weather stations used for interpolations. Fick and Hijmans [32] stated higher prediction error with increase in station elevation and distance to nearest neighboring station, and further reported higher relative error in areas with complex topography (i.e., mountains) and low station density. Compared to WORLDCLIM version 1.4, a higher number of weather station data points (version 1.4: 47,554, version 2.1: 60,419) were incorporated in version 2.1 to interpolate precipitation values. Furthermore, they included cloud cover to increase accuracy in regions with low weather station density and sharply contrasting microclimates to emphasize topographic effects and avoid oversmoothing relationships between elevation and precipitation with stations at equivalent elevations across rain-shadows [61].
Correlation analyses revealed consistent high coefficients among versions and rather low coefficients between climate datasets, whereas the lowest coefficients were found for CHELSA 1.1 and WORLDCLIM 1.4 as well as 2.1 (Table 4). In general, WORLDCLIM shows lower precipitation values compared to CHELSA for Average Precipitation in March, April, and May ( Figure 4). Comparing CHELSA 1.1 and 1.2, the circular arrangement of precipitation amounts remains, with highest values being predicted in the eastern parts of the study area and north of the Himalayan main range. This is suspicious, since the main range constitutes a barrier, leading to very dry climate on the Tibetan plateau. This artefact in CHELSA 1.2 is attributed to the interpolation of station-derived biases. The effect of remote sensing products to decrease bias of precipitation variables can be observed for Average Precipitation in March, April, and May, with artefacts in the eastern part of the study area on the Tibetan Plateau being diminished in WORLDCLIM 2.1 compared to 1.4.

Effects of Climate Datasets on Model Performance and Prediction
Although the modeled potential distribution of Betula utilis stretches along the Himalayan mountain range, we found significant differences in the explained variance between the models using the different climate datasets (Friedmann X 2 = 9.96, p ≤ 0.05), but not among respective versions ( Figure 5). Explained variance was highest for CHELSA 1.2 with 73%, followed by CHELSA 1.1 (72%), whereas WORLDCLIM 1.4 and WORLDCLIM 2.1 received considerably lower explained variances, with 63% and 62%, respectively. CHELSA-related model predictions generally show higher predicted values compared to WORLDCLIM-related predictions, resulting in significant differences between CHELSA 1.   We found conspicuous differences between variable importance for models based on climate datasets, whereas Mean Temperature of the Wettest Quarter gained the highest values for all models ( Figure 6). Precipitation of the Coldest Quarter reached relative importance between 5.4% (CHELSA 1.2) and 11% (CHELSA 1.1). For models based on WORLDCLIM 1.4, variable importance for Mean Temperature of Wettest Quarter and Average Precipitation in March, April, and May was considerably lower (<3%). Most balanced variable importance was found for models based on WORLDCLIM 2.1. The relative importance of Average Precipitation in March, April, and May was always higher for WORLDCLIM-based models than for CHELSA-based models.

Comparison of Future Projections between CHELSA and WORLDCLIM
For demonstration purposes, we projected model predictions under future climate conditions in 2070 using the Community Climate System Model (CCSM4, [58]) of the American National Center for Atmospheric Research (NCAR) for RCPs 4.5 and 8.5 for each climate dataset where future scenarios were available at 30 arc seconds resolution. Differences between the model projections of CHELSA 1.2 and WORLDCLIM 1.4 could hardly be greater (Figure 7). CHELSA future model projection with RCP 4.5 shows a shift of suitable habitat of B. utilis to higher elevations, whereas suitable habitat is no longer projected as a narrow vegetation belt but becomes blurred towards the upper distributional limits. This trend is reinforced under RCP 8.5, where suitable habitat is projected at even higher elevations and on the Tibetan plateau. Based on CHELSA, an increase in potentially suitable habitat for B. utilis is projected under both RCPs. Based on WORLDCLIM future climate projections, a totally different picture emerges with the potential distribution of B. utilis being mainly predicted in west Nepal, and future projections being overruled by the artefact-adhesive data. This leads, for both RCPs, to a projected future distribution on the Tibetan plateau, which is skewed to the western part of the study area. Under RCP 4.5, the potential habitat is projected as more extensive compared to RCP 8.5.

How Does the Methodology of Climate Datasets Influence Model Performance and Predictions?
Freely downloadable climate data with global extent and high resolution (~1 km) have led to a tremendous increase in ecological modeling studies predicting local to global distributional ranges of species under current, future, and past climatic conditions. Such modeling studies are often conducted without evaluation and comparison of climate data

How Does the Methodology of Climate Datasets Influence Model Performance and Predictions?
Freely downloadable climate data with global extent and high resolution (~1 km) have led to a tremendous increase in ecological modeling studies predicting local to global distributional ranges of species under current, future, and past climatic conditions. Such modeling studies are often conducted without evaluation and comparison of climate data quality. The most profound effect on climate data quality can be attributed to solid databases of input parameters (i.e., long-term meteorological records and high density of weather stations). This constitutes the major challenge for global climate datasets since weather stations' distribution is highly uneven due to accessibility and funding, resulting in a poor data basis for regions with high climatic variability (i.e., mountainous regions, coastal areas, and lowland rainforests in the tropics) [62][63][64]. Although both climate datasets use in parts the same raw data (i.e., weather station data) to produce the bioclimatic raster-layers, we found striking differences, especially for precipitationrelated variables. For temperature-related variables, differences were not paramount due to the negative correlation between altitude and temperature, allowing only little room for variation. Further, smaller annual and diurnal ranges (e.g., −14 to 22 • C) of variations in temperatures as well as more bell-shaped distributions of near-surface air temperatures compared to monthly precipitation sums (e.g., 31 to 1620 mm) with highly skewed gamma distributions result in much higher disparities for spatial rainfall climatologies.
The quantity and quality of raw data, applied methodology, and subsequent bias correction are key to understanding differences between the climate datasets. For regions like Europe and Northern America with high densities of weather stations, WORLDCLIM performs well for stations which were used in modeling but fails to predict between stations [34]. Although in WORLDCLIM 2.1 more weather station data were utilized, the interpolation approach can be regarded as inferior to statistical downscaling approaches or climatologies based on remote sensing for ecological approaches [34,63]. Contrary to WORLDCLIM, the challenge of the heterogeneity of spatial and temporal distribution of weather station data incorporated in CHELSA was addressed with a Model Output Statistics algorithm for ERA-Interim reanalysis data, which was corrected using gaugederived products from the GPCC and the GHCN datasets [34].
For remote mountainous regions with low weather station density and complex topography, climate modeling constitutes a challenging endeavor, since local-scale atmospheric processes, which cause site-specific topo-climatic conditions, have to be included in the final product. Generally, local-scale atmospheric conditions are mainly influenced by the underlying terrain, leading to a complex temperature pattern due to anisotropic heating at different slope positions as well as cold air drainage and pooling in mountain valleys during autochthonous weather conditions [35]. Likewise, the spatial pattern of precipitation is influenced by windward and leeward slope positions, leading to hyper-humid climate conditions at the southern declivity of the Himalayan range and semi-arid to arid conditions in the Trans-Himalayan valleys.
For Europe, Marchi et al. [65] assessed good representativeness for mean annual temperature, but high discrepancies for mean annual precipitation between WORLDCLIM 1.4 and climate observations. Generally, they stress that the main shortcoming of WORLD-CLIM 1.4 is the lack of reliable precipitation values due to its low spatial and temporal autocorrelation and missing statistical relationships with physiographic parameters like elevation. The lacking capability of WORLDCLIM 1.4 to capture precipitation in mountainous regions has been debated (e.g., in References [61,[65][66][67][68]). To overcome this challenge, the CHELSA climate dataset includes orographic wind effects and boundary layers in the validation of the main correction step to increase validity of precipitation at stations before downscaling to 1 km resolution [34]. Overall, CHELSA shows more plausible precipitation values compared to WORLDCLIM for our study region, confirming better precipitation patterns compared to WORLDCLIM versions, as stated in Karger et al. [34]. Nevertheless, some limitations to precipitation variables of CHELSA 1.2 exist due to the interpolation of station-derived biases, especially for complex terrain with low densities in either windward and leeward areas (e.g., the Himalayas and the Tibetan Plateau).

Modeling Treeline Species-Oversimplifying Treeline Dynamics?
In general, global climate datasets are rarely robustly evaluated for ecological applications and ecologists often do not account for afflicted uncertainties in their studies [64]. Comparisons and investigations on shortcomings and uncertainties in climate datasets remain largely out of focus in ecological modeling studies (e.g., in References [62][63][64][69][70][71][72][73][74]) and the number of ecological modeling studies in the Himalayas incorporating robust evaluations of climate datasets is very limited (e.g., [42][43][44]). Bobrowski and Schickhoff [44] and Suwal et al. [43] compared CHELSA 1.1 and WORLDCLIM 1.4 under current climate conditions. The former focused on B. utilis, a principal treeline species along the Himalayan arc, while the latter modeled the ecological niche of two Macaca assamensis subspecies in the Hindu-Kush-Himalaya region. Bobrowski and Schickhoff [44] found significantly better results for CHELSA climate data compared to WORLDCLIM. In contrast, Suwal et al. [43] found marginally better results for WORLDCLIM 1.4, however, differences between their model performance criteria were neither statistically validated nor comprehensively discussed in terms of climatological plausibility. Additionally, they modeled potential future range shifts using WORLDCLIM climate data since future CHELSA climate data were not available at that time. Datta et al. [42] compared CHELSA 1.2 and WORLDCLIM 2 under the current climate for modeling Ageratina adenophora in South Asia and concluded marginally better results for WORLDCLIM 2, which were likewise not statistically validated. Although they provided several reasons for their findings, profound climatological explanations in relation to the methodology of the used climate datasets were missing. In order to gauge the potential and challenges of global climate datasets for meso-and macro-scale studies, a critical examination of climate data generation approaches and their ability to capture the underlying climatological processes is mandatory. Our findings expand and update former research of Bobrowski and Schickhoff [44], as most prominent climate datasets and their versions were evaluated for modeling the ecological niche of B. utilis in Nepal with regard to underlying databases, methodological approaches, and their application for ecological models in remote high-altitude regions. In recently published modeling studies on B. utilis in the Himalayas, a critical examination of uncertainties and limitations of global climate datasets has not been carried out (for WORLDCLIM, e.g., References [75][76][77][78][79][80][81]; for CHELSA, e.g., [37]).
We demonstrated the effect of different climate datasets on model performance and prediction plausibility and found major differences between model predictions and serious shortcomings that can be attributed to flawed climate predictors of the respective climate datasets. Model predictions based on CHELSA climate datasets performed significantly better, regarding the explained variance in the test data, than models based on WORLDCLIM (Friedmann X 2 = 9.96, p ≤ 0.05). A comparison among versions of CHELSA revealed no significant differences, while CHELSA 1.2 received the highest explained variance. However, the latest version should not be favored without carefully verifying afflicted distortions of precipitation amounts. To conclude, more plausible precipitation estimates for CHELSA 1.1 compared to CHELSA 1.2 can be found. Thus, CHELSA 1.1 should be preferred when modeling the ecological niche of B. utilis in Nepal. Although Fick and Hijmans [32] highlight significant improvements of WORLDCLIM 2.1 compared to the former version 1.4, we found no significant differences when modeling the ecological niche of B. utilis in Nepal. Especially, the variable Precipitation of Coldest Quarter reflects biased precipitation values similarly skewed towards the Tibetan plateau. We confirmed shortcomings of both WORLDCLIM versions regarding model performance for our case study detected for WORLDCLIM 1.4 in a previous study [44].
Additionally, we found differences in contributions to variable importance between the climate datasets, with all models showing highest importance for temperature during the growing season (i.e., Temperature of Wettest Quarter). These findings are in accordance with the primary dependency of alpine treelines on low air and soil temperatures during growing season [82][83][84]. In general, the position of treelines is controlled by heat deficiency, leading to high susceptibility to climate change impacts [9,85]. For precipitation-related variables, which obtained overall lower importance, Precipitation of the Coldest Quarter was higher for models based on CHELSA climate datasets. The role of thickness and snow cover providing sufficient moisture at the beginning of the growing season for Betula forests has been emphasized [1,[86][87][88]. However, WORLDCLIM climate datasets inherit high biases in precipitation values for the coldest quarter. The climate at treeline sites ranges from more continental in the NW to more oceanic in the SE of the Himalayan mountain system. The amount of annual precipitation increases with increasing monsoonal influence in the same direction along the southern front of the range [15,89]. Average Precipitation of March, April, and May showed higher importance in models using WORLDCLIM, although precipitation amounts are highly patchily distributed, neglecting the degree of continentality towards the eastern part of the Himalayan mountains.
Furthermore, we demonstrated how flawed predictions under current climate conditions impact future projections. The striking differences between CHELSA-and WORLD-CLIM-based projections emerge from incoherent values of the reference period (i.e., current climate conditions). Although CHELSA projects a more conservative distribution in 2070 for both scenarios, extensive care in interpretation should be taken since only one GCM has been projected. Robustness of projections can be increased by comparing multiple GCMs [90]. Wrong assumptions could lead to severe consequences, when these predictions are used to implement conservation strategies without quantification of afflicted uncertainties [91][92][93].

Implications for Future Research
Global climate datasets often neither capture local-scale physical features nor satisfy requirements for climate impact studies [62,64,94]. Furthermore, they have been found to be inferior to local fine-scale climate datasets [69,95,96], whereas climate data with higher resolution (<1 km), capturing microclimates and spatial heterogeneity, shows increased prediction accuracy when modeling species in mountainous regions [97][98][99][100]. Moreover, local, fine-grained climate datasets (250 m) intended for ecological applications have been shown to improve model performance by incorporating climate variability and average extremes [101]. Especially in the course of climate change, the role of microrefugia in topographically complex regions may be emphasized depending on spatiotemporal extent and resolution of the studied area [102,103]. Data on microclimate in mountainous regions is needed to assess changes in species distributions (e.g., extinction, upward shift, range contraction) [104][105][106], but is often not covered in coarse grain, global climate datasets (i.e., CHELSA and WORLDCLIM).
Due to high resolution and spatial and temporal continuity without geographical biases, variables derived from remote sensing products present a valuable option and have been successfully applied to improve ecological models [63,107]. Most notably, information on topo-climate conditions (e.g., elevation, aspect, slope, solar insolation) are vital in mountainous regions, and can be obtained via high-resolution Global Digital Elevation Models (30 m) [108]. In addition, availability and resolution of remote sensing products for ecological applications are advancing, offering great potential for future studies. Currently, a purely remotely sensed bioclimatic predictor set at global extent is available at 2.5 arc minutes [109]. The "classical" modeling approach of bioclimatic and topographical predictors can be enhanced by remotely sensed products of cloud cover [110], snow cover [111], and soil moisture [112] to derive tailored predictors for ecological models. Remotely sensed land surface temperature (LST) and precipitation derivates such as cloud cover were found to provide fair substitutes for climate data derived from CHELSA, when modeling the ecological niche of B. utilis in the Himalayas [36,113]. Moreover, the potential of remotely sensed data as proxies for interpolated weather station data could be confirmed [114]. To generate fine spatiotemporal microclimate consisting of current, future, and past predictors of sufficient extent, Lembrechts et al. [115] proposed a combination of in situ microclimate measurements, long-term weather station data, and high-resolution remote sensing data (e.g., airborne LiDAR and hyperspectral images). Different approaches to reduce errors in climate datasets have been proposed, such as improvement of physics in models, quantification of uncertainty by applying multi-model ensembles, and climate data post-processing to remove model biases [116,117]. However, for the first approach, high research effort is needed, whereas the latter requires tremendous computational power.
Climate data sources and subsequent predictor variables' selection are dependent on the target species and on temporal and spatial extent of the studied region, and may vary between taxa and growth forms, especially in mountainous regions [118]. Incorporation of ecophysiological predictors, other than climate, are needed for complete niche quantification and thus will increase predictive power of the ecological model [119][120][121]. In other words, although treeline elevations can be attributed to prevailing thermal conditions, climate is not the exclusive factor and information on further abiotic and biotic factors driving species distribution are needed to account for species-environment relationships, interactions, and dynamics. To date, an array of additional, standardized predictor variable sets have been published, e.g., CliMond [122,123], EcoClimate [124], ENVIREM [125], MERRAclim [109], and CMCC-BioClimInd [126], varying in spatio-temporal resolution (from~1 to~340 km at the equator). The spatial resolution of ≥1 km may be too coarse to account for differences between slopes in topographically heterogeneous regions like the Himalayas, yet their comparison and evaluation in ecological modeling studies remains an open question.

Conclusions
Climate modeling and thus ecological modeling using climatic predictors face numerous challenges in high-altitude regions. The most important constraint is data availability, which is often sparse due to poor accessibility of the terrain. In many ecological studies, the investigation and evaluation of climate datasets are neglected. Global climate datasets often do not account for local climatic conditions in topographically heterogeneous landscapes due to their quantitatively limited weather station data, computation method, and relatively coarse resolution (≥1 km). Depending on spatial and temporal scales, these climate datasets lack precision of local temperature and precipitation patterns and thus may lead to flawed climate grids used in ecological modeling studies.
Although climate is not exclusively delimiting species niches in mountainous regions, it is an indispensable prerequisite as species occurrences and their dynamics at treelines are governed by low temperatures and local precipitation regimes. Incorporation of further local abiotic and biotic important factors in ecological models are mandatory to account for species' ecophysiological needs, dispersal abilities, biotic interactions, site-specific processes, and their interrelations.
We showed that unmindful usage of global climate datasets may have severe consequences for modeling treeline species in remote high-altitude regions under current climatic conditions, and even more so if flawed current model predictions provide the basis for future projections. These results deepen our understanding of model performance and prediction plausibility related to climate data quality for treeline species in remote high-altitude regions and stress the need for interdisciplinary investigations to advance climate and ecological model quality. Especially in remote high-altitude regions, the underlying data basis of climate models needs urgent improvements, such as quality and quantity of weather station data, predictors, computation methodology, and resolution. However, coherent investigations of afflicted shortcomings and uncertainties in climate data and limitations of ecological model prediction accuracy will remain an essential basis for future modeling studies in high-altitude ecosystems with limited data availability. We advocate a closer collaboration between climate scientists and ecologists to develop high temporal and spatial resolution climate datasets intended for ecological applications, and thus also for conservation strategies.