Downscaled paleoclimate datasets are widely used in biogeographical research, aiding our understanding of past environmental shifts and species’ responses to climate change. Numerous datasets exist, varying in spatiotemporal resolution and underlying methodologies, resulting in variation in estimated temperature. Understanding this variability is important for accurately reconstructing past biogeographical dynamics, especially in complex regions like mountains. We compare the Mean Annual Temperature (MAT) at the Last Glacial Maximum (LGM) from five different downscaled paleoclimate datasets — Beyer, CHELSA-TraCE21k, EcoClimate, PALEO-PGEM-series, WorldClim — against MAT estimates from paleoenvironmental proxy records (fossil pollen and plant macrofossils) within and outside mountains. Additionally, we test the performance of a ‘global grid cooling’ method (i.e. lowering local temperatures by a global LGM estimate) against proxy records. Then, we evaluate the implications of inter-dataset variability for reconstructing temperature-delimited biomes in mountains by reconstructing LGM treeline elevations. We find that LGM temperature cooling and treeline reconstructions strongly vary amongst paleoclimate datasets and between datasets and proxy records. The temperature gradient with elevation is poorly captured by datasets with a coarser spatial resolution. Paleoclimate datasets generally suggest a warmer LGM than proxy records, especially in mountains, while the global grid cooling method more closely aligns with proxy records. Inter-dataset variability can strongly affect the outcome of temperature-delimited reconstructions of biomes and their boundaries, such as treelines. We call for greater awareness and more transparency about the limitations of downscaled paleoclimate datasets in mountainous areas and suggest further research to be aimed at capturing the small-scale heterogeneity of mountains in paleotemperature datasets.

Biogeographical reconstructions, global grid cooling, inter-dataset variability, Last Glacial Maximum, Mean Annual Temperature, mountain regions, paleo-proxy records, paleotemperature, treeline elevation

Paleoclimate models are essential in reconstructing past environmental conditions, providing crucial insights into the Earth’s climate history. Developed over decades, these models simulate past climates by integrating atmospheric, geological and biological data that serve a wide range of scientific disciplines. Paleoclimate simulations (Box 1) enable researchers to understand long-term climate trends, assess the natural variability of the climate system and the impact of climate events on biotic and abiotic systems (e.g. Haywood et al. (2019); Kageyama and Paillard (2021)). For instance, paleoclimate simulations help biogeographers understand historical changes in species distributions (e.g. Brown et al. (2018)), allow cryosphere researchers to track glacial extents and dynamics (e.g. Ludwig et al. (2019)), aid climatologists in studying the climate sensitivity of the Earth system to atmospheric CO₂ (e.g. Annan and Hargreaves (2015); Tierney et al. (2020a)), help geoscientists model the hydrological cycle (e.g. Haywood et al. (2004)) and archaeologists to understand the effect of past climate conditions on the spatial distribution of human populations (e.g. van Andel and Davies (2003)). By simulating climate conditions during key periods, such as the Last Glacial Maximum (LGM), these models offer valuable baselines against which current and future climate changes can be assessed, making them indispensable for our understanding of the Earth’s climatic past.

Since their introduction, paleoclimate models and outputs have evolved substantially and become increasingly accessible to a broad scientific community. The first climate simulations using a General Circulation Model (GCM, Box 1) were published in the late 1970s for the LGM (e.g. Gates (1976)). Proxy records such as ice cores, sediment records and fossil data aided in improving the simulations through testing and validation (e.g. Haywood et al. (2019); Tierney et al. (2020b)) and, more recently, through data-model assimilation (e.g. Tierney et al. (2020a)). Since the early 2000s, the outputs of both GCMs and Regional Circulation Models (RCM, Box 1) have been downscaled and specifically aimed to reconstruct past climate conditions. Through downscaling, the spatial and temporal resolution was enhanced to make it suitable for ecological and biogeographical applications. The availability of these high-resolution datasets facilitated widespread adoption, providing ready-to-use climatic reconstructions for the wider interdisciplinary community. Their global and spatially continuous coverage, provided in gridded format (e.g. raster), were an advantage over the spatially scattered proxy records (Braconnot et al. 2012; Haywood et al. 2019). Additionally, these datasets include multiple consistent variables, opposed to the few variables that can be derived from a proxy. As a result, the datasets of downscaled paleoclimate conditions became an attractive and essential source to assess environmental conditions on regional and continental scales over deeper timescales. In this study, we focus on these datasets and refer to them as “(downscaled) paleoclimate datasets” hereafter (Box 1).

In biogeographical research, the prevalence of downscaled paleoclimate datasets has played a crucial role in enabling researchers to reconstruct historical distributions of species and ecosystems and for understanding the effect of climatic changes on current diversity patterns. Examples of such study topics include, but are not limited to, range shifts (e.g. Martínez-Meyer and Peterson (2006)), migration pathways (e.g. Garzón et al. (2008)), niche stability (e.g. Worth et al. (2014); Kandziora et al. (2024)), refugia (e.g. Svenning et al. (2011); Gavin et al. (2014)), speciation and extinction patterns (e.g. Stigall (2013)) and vegetation or biome reconstructions (e.g. Werneck et al. (2011); Collevatti et al. (2013)). Amongst the most popular datasets are WorldClim (Hijmans et al. 2005) and CHELSA-TraCE21k (Karger et al. 2023), which provide high-resolution (0.04° and 0.008°, respectively) climate data for the period since the LGM. These datasets have been instrumental in numerous studies due to their open accessibility, comprehensive climate variables (“bioclimatic variables”), global and temporal coverage and easy integration with ecological models such as Species Distribution Models (SDMs) and Ecological Niche Models (ENMs) (Varela et al. 2015).

The LGM has been a period of especially high interest for climatic and ecological reconstructions. It represents a time of significant climatic contrast, providing a benchmark for understanding the impacts of extreme climate conditions on ecosystem and species distributions (Varela et al. 2015; Carotenuto et al. 2016). LGM estimates have been frequently used to study, quantify and test, for example, refugia (e.g. Tang et al. (2018)), phylogenetic hypotheses (e.g. Muellner-Riehl et al. (2019)) and species and biome distributions (e.g. Svenning et al. (2011); Kandziora et al. (2024)), often aided by SDMs and ENMs (e.g. Nogués-Bravo (2009); Varela et al. (2011)). Additionally, these estimates are used for studies based on climate derivatives over time, such as climatic stability (e.g. Harrison and Noss (2017); Brown et al. (2020)) and climate change velocity (e.g. Sandel et al. (2011)). Evidently, paleoclimatic data are essential to support numerous applications, from testing ecological and evolutionary hypotheses to understanding species and biome responses to extreme and dynamic climate conditions.

Though paleoclimate datasets all provide global coverage and paleotemperature estimations, they vary vastly in spatial resolution, temporal resolution, and temporal range (Barreto et al. 2023). Moreover, they are based on different GCMs that each have their parameterisations and can, therefore, differ substantially (Varela et al. 2015). Inter-dataset variation can impact findings derived from biogeographical studies (Collevatti et al. 2013), but underlying methodologies and uncertainties are often difficult to grasp for those outside the climatology community (Lima-Ribeiro et al. 2015; Varela et al. 2015). In addition, previous paleoclimatic studies have shown that, while simulations can generally capture the direction of climatic changes from the LGM to the pre-industrial period, they perform poorly in estimating the magnitude of change. In comparison to proxy records, simulations tend to overestimate LGM sea surface cooling in the Tropics (Waelbroeck et al. 2009) and underestimate LGM land cooling (Masson-Delmotte et al. 2005; Kageyama et al. 2006; Bartlein et al. 2011; Braconnot et al. 2012). Harrison et al. (2013) quantified these discrepancies for a range of different simulations, showing that all simulations derived from CMIP5 and PMIP2 consistently underestimate LGM temperatures over land. Mountain regions might be especially prone to potential discrepancies due to larger temperature changes at high elevation during the LGM (Harrison 2005; Loomis et al. 2017), but this has not been tested so far.

Paleoclimate simulation and downscaling in mountain regions present specific challenges (Badgley et al. 2018). Mountains exhibit stark differences in both biotic and abiotic factors over small spatial scales (Shafer et al. 2005). Complex environments and topography create elevational gradients, microclimates and diverse habitats within a close geographical range that are challenging to capture accurately in simulations (Lima-Ribeiro et al. 2015). The more complex the topography, the more temperature information is lost with decreasing spatial resolution of the climate simulation. In addition, mountain regions often have sparse weather stations and limited historical climate data (Hik and Williamson 2019; Pepin et al. 2022), hindering accurate calibration and validation of downscaled climate data in mountains and resulting in less reliable predictions (Thornton et al. 2022). Given the wide usage of paleoclimate datasets by the biogeographical community, it is essential to assess their applicability in mountains from the perspective of biogeographical reconstructions.

In this study, we assess paleoclimate data variability in mountains by comparing estimated Mean Annual Temperature (MAT) at the LGM derived from five different downscaled paleoclimate datasets — Beyer, CHELSA-TraCE21k, EcoClimate, PALEO-PGEM-series, WorldClim — (Table 1) against a validation compilation of paleoenvironmental proxy records. Our objective is to investigate how these temperature estimations differ worldwide and test if substantial difference exists between paleoclimate datasets and within and outside mountain ranges. We further compare outcomes with a so-called ‘global grid cooling’ method where we lower the world’s temperature by a global LGM estimate. To showcase the impact of inter-dataset variability and its potential implications for biogeographical research, we reconstruct treelines at the LGM using the MAT estimations of the paleoclimate datasets. Treelines, which typically occur at high elevation in mountain ranges, but also include the polar treeline, delineate different biomes and serve as important ecological boundaries and transition zones (Jobbágy and Jackson 2000; Körner 2007; Cairns 2013; Greenwood and Jump 2014). This makes treelines an excellent subject for comparing paleoclimate datasets because they are climate-delimited with a particular sensitivity to temperature (Paulsen and Körner 2014; Körner 2021), thus relying heavily on accurate temperature estimations. By focusing on the treeline in mountains globally, we aim to provide insights into the variations and consistency of paleoclimate datasets in mountainous areas, informing future biogeographic and ecological research.

Our analysis reveals several key findings about LGM temperature estimations from downscaled paleoclimate datasets and their impact on treeline reconstructions. Firstly, paleoclimate datasets exhibit substantial variability in spatial resolution, indicating that the precision of these datasets can vary considerably, particularly in mountain ranges. Secondly, we observe that paleoclimate datasets tend to generally estimate higher temperature values for the LGM than proxy records, thus an underestimation of LGM cooling by these datasets. Although not significant, this difference appears to be larger within mountains than outside mountains, meaning that LGM temperatures within mountain ranges are too warm compared to proxies. Thirdly, the bias in simulated LGM temperature estimates varies substantially across the globe. Paleoclimate datasets generally overestimate LGM cooling within mountains towards the Northern Hemisphere, compared to a relatively stable relationship towards the Southern Hemisphere, indicating a hemispheric bias in temperature estimations. Still, the direction of dataset-proxy disagreement often concurs amongst datasets, highlighting a consistent pattern in dataset behaviour despite regional differences. Fourthly, treeline reconstructions exhibit large elevational ranges within datasets and stark differences between datasets. Finally, coarser resolution datasets fail to capture smaller mountain ranges and islands in treeline reconstructions, underscoring the limitations of these datasets in representing finer geographical features.

Paleoclimate datasets vary vastly in spatial resolution. It is evident from our spatial comparison that not all datasets that cover the LGM have a suitable resolution for studies in mountain ranges (Fig. 3A). For example, coarser resolution datasets (Beyer, EcoClimate and PALEO-PGEM) fail to capture the gradual change of temperature with elevation, leading to counter-intuitive results (Fig. 3B). Limitations of paleoclimate models in mountains have been noted before, particularly their challenges in accurately simulating “regions with high gradients” (Brown et al. 2018), “regions with complex topography” (Lima-Ribeiro et al. 2015) or “topographic peculiarities” (Gamisch 2019). The more complex the topography, the more temperature information is lost with decreasing spatial resolution of the climate simulation. As mountain regions are of particular interest for the biogeographical community (Perrigo et al. 2020), it is important to consider the effect of the underlying climate simulation’s spatial resolution on downscaled temperature estimations. Unfortunately, the varying temporal range and resolution of the datasets (see Table 1) could mean that the highest-resolution dataset available for the period of interest is insufficient to perform a sensible analysis.

Higher resolution, however, does not guarantee a more accurate representation of reality. For instance, downscaling methods allow researchers to computationally increase the spatial resolution, resulting in pseudo-replicated information (Taylor et al. 2012; Lima-Ribeiro et al. 2015). This is observed in Fig. 4, where CHELSA-TraCE21k showed the largest deviations from proxy estimates, despite having the highest spatial resolution. Additionally, the elevation-temperature relationship shown by the CHELSA-TraCE21k dataset has a surprisingly sharp cut-off (Fig. 3B). The clustering of datapoints on the left side of the trend line in CHELSA-TraCE21k’s scatterplot and the absence of data points on the right could reflect specific downscaling parameters such as thresholding, filtering or interpolation techniques. Such specificities are handled in all datasets and influence the distribution and range of temperature values across elevations (Varela et al. 2015). Furthermore, statistical downscaling (Box 1) relies on several assumptions that differ from reality, such as a stationary relationship between macro-and microclimates and a constant lapse rate over time, which can lead to inaccuracies in downscaled simulations (Salvi et al. 2016). This is also the case for CHELSA-TraCE21k, which is a statistical downscaling of the ERA-Interim reanalysis, with the temperature downscaling based on mean lapse rates and elevation (Karger et al. 2017). While these assumptions and techniques can improve spatial resolution, they also risk introducing biases that compromise the accuracy of temperature-elevation relationships, highlighting the need for caution when interpreting downscaled datasets.

We observe substantial spatial variability in LGM temperature estimates amongst different paleoclimate datasets, with a pronounced discrepancy between datasets and proxy records, especially in mountainous regions. These mismatches are not surprising, given the inherent difficulty of modelling and downscaling with a multitude of specific mechanisms and spatial variations at play (Harrison and Bartlein 2012). Overall, LGM cooling is underestimated by paleoclimate datasets (Fig. 4), which is in line with earlier findings of Kageyama et al. (2006), Masson-Delmotte et al. (2006), Bartlein et al. (2011) and Braconnot et al. (2012). Although not significant, this effect of underestimation is observed to be stronger within than outside mountain ranges for paleoclimate datasets. Such new insights are intuitive given that proxy records have indicated larger temperature changes at higher elevations (Harrison 2005; Loomis et al. 2017) and topographic-rich regions have been notoriously challenging to model (Lima-Ribeiro et al. 2015; Brown et al. 2018; Gamisch 2019). Furthermore, the dataset’s temperature differences from proxy records are not randomly distributed in geographical space (Fig. 5). In South America, Central Europe and South Africa, LGM temperatures estimates from datasets are generally higher than proxy records (i.e. underestimation of cooling). In Alaska, Norway and Russia, LGM temperature estimates from datasets are generally colder than proxy records (i.e. overestimation of cooling). Varela et al. (2015) observed similar patterns for temperate regions when comparing GCM predictions for the LGM. Our study highlights that, also for downscaled paleoclimate datasets and especially for mountain regions, regional biases exist and should be considered when interpreting research outcomes, based on such datasets.

We also observe some outliers in the proxy data, located in the northernmost part of the Northern Hemisphere, with exceptionally large negative ∆T_difference values. These indicate that, for these locations, proxy records estimated higher LGM temperatures than the datasets, namely 8 and 13 °C. Many of these outliers occur at or near the former position of the Northern Hemisphere ice sheets, suggesting that a potential cause for these outliers is a modelled misspecification of the extent of these ice sheets, likely due to limited spatial resolution. This effect was possibly captured by the proxy records, but not by the simulations that underlie the datasets. Since these outlier records were present in every dataset we included, we believe it did not affect the observed inter-dataset variability.

Our results reveal that the global grid cooling method, despite its simplicity, has smaller incongruences within mountain ranges than outside in comparison to the proxy records (Figs 4, 6). A possible explanation for this is that the global grid cooling method applies a uniform temperature reduction across the globe without the introduction of dataset-specific biases, leading to more consistent and predictable results in regions with complex topography, such as mountain ranges. Still, it is surprising to see that this method has a closer similarity to the proxy-inferred LGM temperatures than CHELSA-TraCE21k’s LGM temperature estimations. The difficulties in reconstructing temperature in mountains stem from the challenges of modelling complex topography, the scarcity of proxy data, variability in temperature lapse rates and microclimatic effects, amongst others (Gavin et al. 2014). It would be of interest to further investigate the potential of this method for mountain ranges, for example, by using regional temperature reconstructions derived from paleoenvironmental records (e.g. Chala et al. (2017); Rangel et al. (2018); Flantua et al. (2019)) — albeit as a heuristic approach until more robust paleoclimate datasets for mountains become available. Using regional temperature reconstructions has already shown great potential for spatial mapping of ecosystem shifts, for example, in island (Rijsdijk et al. 2014; Norder et al. 2018) and mountain biogeography (Flantua et al. 2019). Our results underscore that simplicity does not necessarily have to come at the expense of quality, offering many opportunities for environmental and biogeographical research depending on paleoclimatic reconstructions.

Regional differences in LGM cooling estimates, especially around the Tropics, are much less pronounced for the global grid cooling method compared to the paleoclimate datasets (Fig. 5). This might be related to the fact that there is often a lack of high-quality, long-term climatic and proxy data in tropical regions compared to temperate regions (Harrison 2005; Tovar et al. 2022). This scarcity of data can hinder the calibration and validation of simulations, making it more difficult to produce accurate LGM simulations. Additionally, tropical mountain regions, such as the Northern Andes, are characterised by complex climate dynamics with diverse vegetation, making it difficult to accurately capture these interactions (Tovar et al. 2022). In addition, temperature changes between the present and LGM were generally larger in temperate regions compared to the Tropics, due to the extensive ice sheets covering large parts of North America and Eurasia and polar amplification (Cao et al. 2019). Addressing these challenges by improving the calibration of paleoclimate simulations and incorporating more diverse proxy records from tropical mountains could significantly enhance our understanding of LGM climate dynamics in these complex regions.

There are substantial disagreements amongst paleoclimate datasets and across mountain ranges in treeline reconstructions, as shown by our results (Fig. 7B). The treeline reconstructions follow the well-known general pattern of lower treelines towards the poles and higher treelines around the Equator. Yet, the differences between datasets for a single mountain range are substantial (Fig. 7B). Depending on the dataset, the final reconstructed treeline can vary thousands of metres. Although some of the broad ranges between the minimum and maximum treeline reconstructions can be attributed to mountain ranges spanning over many latitudinal degrees (e.g. Southern Andes), others cannot (e.g. Selenga Highlands).

Underestimating LGM cooling places treelines too high, reducing the estimated alpine biome surface area, while overestimating LGM cooling places treelines too low, thus inflating the estimated alpine biome surface area. The magnitude of these over- and underestimations vary from dataset to dataset and per region within a dataset. Considering how dataset choice influences the outcomes of the treeline reconstructions, such variations are equally expected to occur with any other boundary or biome reconstruction in mountain ranges. This has far-reaching implications for research concerning, amongst others, SDMs and ENMs, where species and biome range reconstructions are an important outcome to derive conclusions on, for example, refugia, climatic stability and velocity (e.g. Seo et al. (2008); Schorr et al. (2012); Gavin et al. (2014); Domic and Capriles (2021)). Inter-dataset variability (i.e. different datasets estimating different treeline elevations) and within-dataset divergence (i.e. broad estimated elevation range within a dataset) stresses the critical role of dataset choice on ecological reconstructions.

Modelling limitations are ubiquitous and perhaps even inevitable (Morrison and Lawrence 2023). The aim of this research is not to determine which dataset is better, but rather to highlight the magnitude and direction of intra- and inter-dataset variability and illustrate the implications for biogeographical research. It is difficult to determine whether the observed incongruences are caused by the underlying simulations of GCMs or RCMs, observations, statistical downscaling methods, or something else entirely. It is possible, for example, that the downscaling methods fail to capture elevation-dependent cooling in mountains ranges, which could explain the underestimated LGM cooling. Since the original simulations from ESMs, GCMs and RCMs are relatively coarse, the downscaled products could give false precision by trying to simulate elevational variations for high-resolution topographies that are absent in the original grids, leading to artefacts and biases in the downscaled data. It is also possible, however, that the biases were not introduced through downscaling, but already present in the original simulations and simply passed on to the downscaled dataset. Harrison et al. (2013), for example, found that paleoclimate models consistently underestimate land cooling and overestimate ocean surface cooling during the glacial period. Harrison et al. (2013) also found that models generally capture the direction of climate change, but struggled to reproduce spatial patterns, which is in line with the findings of this research. The datasets that were compared in this research, however, are based on various climate models and downscaling approaches (see Table 2), making it harder to pinpoint from where exactly the discrepancies come. Additional research is needed to investigate the underlying causes of the observed discrepancies.

Regarding paleoenvironmental proxy records, this study included proxies from a broad geographic distribution of records. Nonetheless, not all mountain ranges were represented by proxy records and, generally, only a single record was available per mountain range. While this limitation does not affect the observed differences between proxy records and datasets for a specific location, the addition of further records, especially for high elevations, could potentially refine the trends’ strength and orientation currently observed. Furthermore, there are other potential proxy records (e.g. chironomids) and multi-proxy approaches that could improve the temperature simulations of paleoclimate models. When it comes to reconstructing past environments in mountains, regional paleoenvironmental proxy records can provide crucial information, but we should also consider the following: 1) despite global availability of such records, researchers cannot always rely on having sufficient records available in or near their study areas for the required time period. This is one of the reasons why downscaled paleoclimate datasets are so important and popular (Haywood et al. 2019); 2) proxy records also come with uncertainties and methodological limitations and cannot be treated as absolute truths (Harrison et al. 2016; Tierney et al. 2020a). For example, lower pCO₂ during the LGM could have had a negative effect on the water-use efficiency of plants and, thus, affect the climate-vegetation relationship (Street-Perrott et al. 1997). This could potentially make pollen assemblages appear colder (Cowling and Sykes 1999). However, while there seems to be an agreement on the impact of lower pCO₂ on pollen-based reconstructions of moisture availability (the reconstructions may be too dry relative to reality; Chevalier et al. (2020); Scheff et al. (2017)), there is no consensus about whether lower pCO₂ affects pollen-based reconstructions of temperature (Williams et al. 2000).

Combining a suite of proxy data with (multiple) paleoclimate simulations, a method called data assimilation, has been suggested as the best way to obtain LGM paleotemperature estimations for specific regions of interest (Schneider von Deimling et al. 2006; Worth et al. 2014). Paleoclimate data assimilation allows the refinement of model states, theoretically leading to more reliable climate reconstructions, as shown by, for example, Tierney et al. (2020a) and Osman et al. (2021). Although the integration of multiple simulations with proxy information helps reduce uncertainties in the climate estimates, this method remains dependent on the quality and spatial resolution of the underlying simulations. Furthermore, Williams et al. (2022) warn of circularity issues, where proxy networks, derived from certain biological proxies included in assimilation, can preclude the use of the resultant temperature estimates for studying, for example, ecoclimate sensitivity of the same proxies. Additionally, the inherent challenges of reconstructing climate at high resolution in topographically complex regions (i.e. mountains) or for regions or periods with sparse proxy records, also apply for data assimilation. Methods to tackle these issues are already in development, (see Tierney et al. (2020a) and Osman et al. (2021)), but still require further analysis.

We have highlighted the discrepancies for the LGM temperatures and illustrated the impacts for treeline reconstructions, but it is likely that the observed issues occur in other regions, moments in time and for other biome boundaries as well. Until paleoclimate datasets provide more reliable LGM temperature estimates for mountainous regions, the biogeographical community should remain cautious of these discrepancies. Researchers should consider these dataset limitations when interpreting climate reconstructions in mountain regions and use multiple datasets and methods to cross-validate their findings. We further suggest avoiding datasets with strong variation or bias in regions of interest. Continued efforts to refine paleoclimate datasets and incorporate newly available proxy compilations from mountains and tropical regions will be essential for enhancing the reliability of these critical reconstructions. Such improvements are key to advancing our understanding of past biogeographical patterns and ecological dynamics that have shaped mountain ecosystems through time.

Downscaled paleoclimate datasets are extensively used in environmental and biogeographical research, with the Last Glacial Maximum (LGM) being a critical period due to its unique climatic conditions and significant impact on global ecosystems. However, the LGM temperature estimations from these datasets often differ from paleoenvironmental proxy records, with larger discrepancies observed within mountain ranges. A global grid cooling method, which lowers local temperatures by a global LGM estimate, provides closer-to-proxy estimations within mountains compared to downscaled paleoclimate datasets. Furthermore, coarse resolution datasets fail to capture the temperature gradient with elevation in mountainous regions.

We also highlight the disagreement within and amongst datasets and demonstrate that inter-dataset variability can substantially affect the outcomes of temperature-delimited boundary reconstructions of biomes. Not all paleoclimate datasets are suitable for studies on mountain ranges and significant incongruences and inconsistencies must be considered. Therefore, we call for greater awareness amongst the interdisciplinary community regarding the limitations of paleoclimate datasets, especially in mountainous areas. We advocate for more transparency and guidance in the methodologies of downscaled paleoclimate datasets for a broad audience and suggest that further research should focus on developing methods specifically tailored to capture the complex heterogeneity of mountain ranges, providing precise and reliable temperature estimations across time and space.

ESR, LS, ATW, JS, JAG and SGAF acknowledge financial support from Trond Mohn Research Foundation (TMF) and the University of Bergen for the start-up grant ‘TMS2022STG03’ to S.G.A. Flantua. JS additionally acknowledges financial support from Research Foundation – Flanders (FWO) Odysseus Type II project GlaciersMD. OM is funded by the Czech Science Foundation PIF grant (GN23-06386I), by the Charles University Research Centre programme (UNCE/24/SCI/006) and by the Institutional Support for Science and Research of the Ministry of Education, Youth and Sports of the Czech Republic. We also thank two anonymous referees for their detailed feedback, which led to significant improvements of the publication.

The authors have declared that no competing interests exist.

ESR: Conceptualisation, methodology, software, validation, formal analysis, investigation, data curation, writing-original, visualisation. OM: Investigation, writing-review and editing. CP: Investigation, writing-review and editing. LS: Investigation, writing-review and editing. JS: Investigation, writing-review and editing. ATW: Methodology, software, data curation, writing-review and editing. JAG: Investigation, writing-review and editing, supervision. SGAF: Conceptualisation, methodology, resources, investigation, writing-review and editing, supervision, project administration, funding acquisition.

Code and data required to perform the ∆T comparison, treeline reconstruction and all visualisations are available as a GitHub release on Zenodo: https://doi.org/10.5281/zenodo.14228272.