Late Quaternary fluctuation in upper range limit of trees shapes endemic flora diversity on the Tibetan Plateau

Domain: earth_science
Published: 2025-02-20
Editorial outcome: major_revision

Abstract
The influence of paleoclimate in shaping current biodiversity pattern is widely acknowledged. However, it remains unclear how the upper paleo-range limit of trees, which dictated the habitat of endemic alpine species, affects the variability in endemic alpine species composition across space over the Tibetan Plateau. We integrated satellite-derived upper range limit of trees, dendrochronological data, and fossil pollen records with a paleoclimate dataset in a climate-driven predictive model to reconstruct the spatio-temporal upper range limit of trees at 100-year intervals since the Last Glacial Maximum. Our results show that trees distributed at the lowest elevations during the Last Glacial Maximum (~3426 m), and ascended to the highest elevations during the Holocene Climatic Optimum (~4187 m), a level ~180 m higher than the present-day (~4009 m). The temporal fluctuations in paleo-range limits of trees play a more important role than paleoclimate in shaping the current spatial pattern of beta-diversity of endemic flora, with regions witnessing higher fluctuations having lower beta-diversity. We therefore suggest that anthropogenic-caused climate change on decadal-to-centennial timescales could lead to higher fluctuations in range limits than orbitally-forced climate variability on centennial-to-millennium timescales, which consequently could cause spatial homogenization of endemic alpine species composition, threatening Tibetan endemic species pool.

Introduction
The high mountain ranges of the Tibetan Plateau are global biodiversity hotspots with a rich endemic diversity in alpine flora1–3, which constitute a high proportion of regional diversity4,5. It is known that local endemic alpine species richness (alpha diversity), which typically have small effective population sizes and limited dispersal capacities6,7, is threatened by climate change8,9 and the upward movement of trees10–13. The upward migration of tree species threatens alpine ecosystems through habitat loss and fragmentation, as well as alterations in microclimate (e.g., increased sheltering effect). To understand how best to conserve endemic alpine species diversity across spatial scales, it is imperative to understand the driving mechanisms of beta-diversity, as it is the component of regional gamma-diversity that accumulates as a result of compositional differences between alpha diversity14–18. The pattern of contemporary beta-diversity, which quantified the variation in species composition within a region, is generally interpreted from past climate regimes and their oscillations19,20, but the role of past changes in the upper range limit of trees remains unknown. Specifically, in areas with stable climate conditions, low fluctuation in the position of the upper range limit of trees would result in reduced habitat loss and fragmentation, thereby creating stable environmental conditions that are more conducive to the survival of endemic alpine species, compared to areas where the position of upper range limit of trees fluctuates more significantly. We should therefore expand our analysis of the drivers behind the current spatial patterns of beta-diversity, shifting focus from climate change to include a more immediate factor—range limit dynamics. This broader understanding could offer valuable insights for more accurately predicting future shifts in alpine endemism diversity under a rapidly warming world.

Here the upper range limit of trees, which is defined as a connection between the uppermost group of upright trees of a certain height12,21, could theoretically reach their physiological limits set by temperature (referred to as the treeline position)21,22 but it is often confined to lower elevations due to environmental stresses and stochastic effects such as human disturbance12,21–24. Our previous understandings of paleo changes in the upper range limit of trees are mostly derived from analysis of fossil plant data including pollen, macrofossils, and charcoal. For example, classical paleoecological records and ancient DNA in lake sediments indicate that, the upper range limit of trees at some locations may have reached 3500 m elevation during the LGM25, and 400–600 m higher than today in the Holocene Climate Maximum26,27. While, relying solely on fossil data at the site level to reconstruct the spatio-temporal evolution of the upper range limit of trees remains highly uncertain, especially on the Tibetan plateau where a handful of sites were available for a few time periods since the Last Glacial Maximum (LGM, ~22 kyr BP). Evidence from paleoecological records show that modern trees survived the extremely harsh climatic conditions during the most recent glacial period25. This knowledge opens up the possibility of reconstructing the upper paleo-range limit of trees from paleoclimate data based on the modern climatic responses of the upper-range limit of trees. Notably, the existence of cryptic forest refugia, which provide long-term shelters for trees, could facilitate the application of a climate-driven modeling approach for paleo reconstruction28. These refugia, found in heterogeneous topographies within high mountain areas and intersecting valleys of larger rivers on the Tibetan Plateau25,29, could result in high dispersal rates and enable tree species to respond within decades to deglacial climate change30.

Here we derived past dynamics in upper range limits of trees using paleoclimate data at a spatial resolution of ~1 km at 100-year intervals since LGM, using the modern relationship between the satellite-derived upper range limit of trees and climate. The period since the LGM was used because its significant paleoclimatic changes led to notable vegetation and floristic shifts, as evidenced by fossil records preserved in sediments19,31–33, potentially driving legacies in modern diversity patterns. In addition, our paleo-reconstruction of the upper range limit of trees was based on the assumption that the relationship between the upper range limit of trees and climate, extracted from modern observations, has remained constant through time since LGM21,22,34. This assumption, which is well known as the uniformitarian principle, has been widely used in previous pollen-based paleoclimate reconstructions35,36. We further estimated how paleo fluctuation in the upper range limit of trees affects the current spatial pattern of beta-diversity in alpine endemic flora over the Tibetan Plateau.

Results
Spatial pattern of current upper range limit of trees and its drivers
We delineated the spatial distribution of the current upper range limit of trees using satellite-derived forest extent data at a high spatial resolution of 30 m12 (see Methods). To validate the robustness of this satellite-derived upper range limit of trees, we generated 602,879 manually interpreted samples using Google Earth high-resolution (<5 m) images as an independent validation dataset, with 5.1% (N=30,746), 84.6% (N=510,035), and 10.3% (N=62,098) located in the northern, southern and inner Tibetan Plateau, respectively12,37. Across all samples, the satellite-derived upper range limit of trees showed a highly accurate alignment with manual interpretations of Google Earth high-resolution images (R2=0.98, Slope=0.98, ME=17 m, N=602,879) (Fig. 1b).

Our analysis of the satellite-derived upper range limit of trees reveals that almost 94.0% (N=11,284,982 pixels at a 30-m resolution) of pixels depicting the upper range limit of trees are distributed in the northeastern and southeastern margins of the Tibetan Plateau, including the southern Qilian (3294±314 m, mean elevation±SD), eastern Himalaya (4206±342 m) and southern Hengduan (3995±355 m) mountain ranges. The remaining 6.0% (N=721,452) of pixels depicting the upper range limit of trees are scattered throughout the headwater region of the three rivers (Yangtze, Yellow, and Lancang River) (3971±353 m) and Tanggula mountain range (4211±222 m) (Fig. 1a). The upper range limit of tree elevation displays a latitudinal gradient, with elevation in the southern region (4039±359 m) higher than that in the northern (3295±315 m) and inner Tibetan Plateau (3859±387 m). The mean elevation of the upper range limit of trees over the Tibetan Plateau has reached 4009±383 m, which is much higher than that of the Himalayas12 (3633±475 m), and 17 locations in four counties (e.g., Bayi, Gongbu Jiangda, Milin and Lang) have higher elevations than highest records according to field observations38 (4900 m).

Furthermore, we explored the drivers of spatial changes in the current upper range limit elevation over the Tibetan Plateau. Trees do not always reach their thermal treeline positions approximated by an isotherm of mean growing-season temperature18,21,22,24 (see Methods), and they often remain confined to lower elevations due to human activities12,21,39 (such as grazing, logging and/or burning) or other disturbances, or critical water shortage. Our results show that nearly 85.2% of pixels have agreement between the upper range limit of the tree and the thermal treeline21,37 positions over the Tibetan Plateau, suggesting a predominant role of temperature in determining upper range limits of trees. In the southern and inner Tibetan Plateau, nearly 87.5% and 51.8% of pixels have their upper range limit of trees reaching thermal treeline positions, respectively, while only 27.9% of pixels in the northern Tibetan Plateau reached their thermal positions. We then explored which factors could explain the remaining 14.8% of pixels that fall below the thermal treeline, especially in the northern Tibetan Plateau (Supplementary Fig. 1 and Supplementary Table 1). To do so, we calculated the elevational difference between the satellite-derived upper range limit of trees and thermal treeline12,21,22, denoted as D_Treeline, and used a random-forest algorithm40 to rank the importance of several variables (N=19) belonging to categories of environmental stresses and disturbances on D_Treeline (see Methods). Our results showed that D_Treeline has a low value in the southern (54±162 m) and inner Tibetan Plateau (186±223 m), but a high one in the northern Tibetan Plateau (271±217 m). Both cumulative climatic water deficit (CWD) during spring (March to May) and vapor pressure deficit (VPD) were the most important variables driving the high D_Treeline values in the northern Tibetan Plateau, suggesting a regional drought-induced absence of trees from the thermal treeline position12,41 (Supplementary Fig. 1). Furthermore, we assessed how spatial changes of the leading factors (CWD and VPD) modulated the spatial distribution of D_Treeline. For example, the abrupt D_Treeline transition occurring at the CWD threshold of ~100 mm yr−1 (Supplementary Fig. 1b) separated the southern Tibetan Plateau, with low CWD and D_Treeline values, from northern Tibetan Plateau with high ones at ~35°N (Fig. 1a). The transition of D_Treeline along the VPD gradient further sharpens the observed south–north contrast, with D_Treeline exhibiting a pronounced increase when the VPD exceeded 0.4 hPa (Supplementary Fig. 1c).

Developing a predictive model of the upper range limit of trees
We developed a predictive model to reconstruct and project changes in the upper range limit of trees. As there are time lags in seed dispersal and tree establishment, the shift in the upper range limit of trees may not occur in synchrony with climate change21,34,42–44. We included both climate change and migration lag in our predictive model. We used growing season temperature along with climatic stressors (especially spring CWD and VPD) to parameterize the effect of climate change on temporal shifts in the upper range limit of trees, and developed a climate-based predictive model using a training dataset (N=116,337) (see Methods). Our climate-driven predictive model could explain nearly 84% of spatial variation of the upper range limit of trees using an independent validation dataset (N=11,634) (Supplementary Fig. 2).

We then incorporated the impact of migration lag on shifts in the upper range limit of trees into the climate-driven predictive model. To determine the optimal lag, we tested time lags from 10 years to 100 years in our climate-driven predictive model, by minimizing the difference between simulated and in situ shifts in the upper range limit of trees over the past century across 20 sites of the Tibetan Plateau. These range limit elevation dynamics were reconstructed from tree ring data using the dendrochronological method39,45–52 (Supplementary Fig. 3 and Supplementary Table 2). Among all-time lags tested, incorporating a 50-year time lag into the climate-driven predictive model maximizes its explanatory power in representing in situ shifts in the upper range limit of trees over the last century. Specifically, the predictive model explains nearly 70% of the spatial variation in range limit shifts across 20 sites on the Tibetan Plateau (R2=0.68, P<0.01) (Supplementary Fig. 4). Furthermore, after including the time lag, the trees across all sites are predicted to shift upward at a rate of 20 ± 21 m decade−1, which becomes closer to in situ shifts (18 ± 23 m decades−1). This confirms the robustness of our constructed model in the prediction of the upper range limit of trees at the century timescale (Supplementary Fig. 3b).

Prediction of upper range limit of trees changes from the LGM to the end of this century
We used this predictive model to reconstruct spatial-temporal changes in the upper range limit of trees at a temporal resolution of 100 years from the LGM to the present, based on a paleoclimate dataset53. We then projected changes to the end of this century (2080–2099) under the three emissions scenarios (see Methods). This model mapped upper range limit of trees at a spatial resolution of ~1 km over time.

According to the predictive model, the reconstruction of upper range limit elevation spanning the past 22 kyr BP to the end of the century varied by more than 850 m, which is around one-fifth of the present-day range limit elevation (4009±383 m) (Fig. 2). Across the past 22 kyr BP, the elevation of the upper range limit of trees reached the lowest position during the LGM (22–19 kyr BP), gradually peaking in the warmest period of the Holocene (HCO, roughly 10–6.5 kyr BP), and then gradually decreasing towards the present day (4009±383 m). By the end of this century, trees are predicted to migrate to values comparable to those during HCO (4185±362 m) under the SSP2–4.5 scenario, while under the SSP5–8.5 scenario, they are expected to reach elevations ~240 m higher than the present day.

Our model predicted that elevational changes in the upper range limit of trees reveal a distinct four-phase pattern, generally corresponding to changes in growing-season temperature. The first phase covers the LGM period (22–19 kyr BP), when the mean growing-season temperature was 2.7 °C lower than the present (Supplementary Fig. 5). The upper range limit of trees then reached the lowest elevation (3426±347 m), which is ~583 m lower than the present day (4009 m). During the second phase from the Last Deglaciation period (19–12 kyr BP) to the middle Holocene (9–5 kyr BP), the climate transitioned from a glacial cold to a Holocene warm period, and the upper trees gradually increased towards the highest predicted elevations (4187±377 m), which is 178 m higher than the present (4009 m). Note that this gradual increase in the upper range limit of trees was interrupted by the two millennial-scale events, the Bølling-Allerød warm event (14.7–12.9 kyr BP) and Younger Dryas cold event (12.9–11.7 kyr BP), which caused a marked increase and decline in elevations, respectively.

Throughout the third phase from Late Holocene (5–0 kyr BP) until the preindustrial era, the climate underwent a period of steady cooling and fluctuating warm and cold phases, resulting in a gradual decrease of the upper range limit of trees that eventually stabilized at a value ~61 m lower than the present (4009 m). The fourth phase, referred to the anthropogenic warming period from the pre-industrial era to the end of this century, has and will continue to experience a steady upslope of the upper range limit of trees. Relative to the present, the mean elevational increase in the range limit over the Tibetan Plateau are 102 m, 176 m, 243 m for SSP1–2.6, SSP2–4.5, and SSP5–8.5, respectively (Fig. 2 and Supplementary Table 3). We also calculated the alpine land area above the mean reconstruction of upper range limit of trees (see Methods). Our results show that the alpine land area reached a minimum extent of 503,035 km2 during HCO and a maximum extent of 849,574 km2 during LGM, 75% and 150% of the present-day area of alpine habitat (597,788 km2), respectively (Supplementary Fig. 6).

To verify the reliability of our model in reconstructing historical changes in the upper range limit of trees, we used fossil pollen abundance records from lake sediments to determine the presence and absence of trees at sampling sites throughout the Holocene. By integrating data from six lakes situated below, at and above the present upper range limit of trees26,35,54–57 (Supplementary Fig. 7 and Supplementary Table 4), this elevational transect of fossil pollen records provides observational evidence of changes in the upper range limit of trees over the Holocene. Specifically, for a given paleo time period, we obtained the elevational range of the upper range limit of trees by identifying the presence and absence of trees from fossil pollen records in lake sediments at different elevations. For example, during the early to middle Holocene (~10–4 kyr BP), trees were detected in Wuxu Lake (3710 m), Muge Co Lake (3780 m), and Yidun Lake (4470 m), but were absent in Koucha Lake (4540 m) and Saiyong Co Lake (4545 m) (Supplementary Fig. 7b). We then obtained that the elevation of paleo-upper range limit was situated between 4470 m (Yidun Lake) and 4540 m (Koucha Lake) during ~10–4 kyr BP. By further comparing elevational changes between different paleo time periods, we could obtain paleo changes in upper range limit of trees using fossil pollen records. To compare such pollen-derived data with our model reconstruction, the same grid cells covered by the pollen-derived data are selected from our high-spatial-resolution reconstructed data (Supplementary Fig. 7). The pollen-derived data are consistent with our model reconstruction showing that the elevation of upper range limit of trees increased in the early Holocene, held a consistent elevation between ~10 kry BP and ~4 kry BP during the middle Holocene and then decreased in the late Holocene (Supplementary Fig. 7b). The pollen-derived data further reveals a descent of upper range limit of trees by nearly 270 m between the early Holocene and the present over the Tibetan Plateau, which is also consistent with our model results (239 m). Interestingly, this range shift is also consistent with that of the European Alps during the Holocene (200–250 m) using macrofossil data24,58. However, determining precise elevational shifts in the upper range limit of trees remains a challenge due to the relatively coarse elevational resolution of fossil pollen records along our constructed elevation transect. Thus, it becomes necessary to either increase the elevational and temporal resolution of fossil pollen records or include analysis of terrestrial plant macrofossil and macroscopic charcoal data at and above the upper range limit of trees.

Impact of paleo-range limit fluctuation on current spatial pattern of beta-diversity
Next, we assessed the legacies of paleo changes in the upper range limit of trees on the spatial pattern of beta diversity over the Tibetan Plateau. For this purpose, we collated elevational and spatial distributions for 2111 endemic alpine species in around 100 counties over the Tibetan Plateau (Supplementary Data 1), with the upper limit of species distributions extending past the current upper range limit of trees (see “Methods”). We used Sørensen-based multiple-county dissimilarity59 to calculate the beta-diversity (β-diversity) of alpine endemic flora within a region that is defined to consist of a target county and its surrounding counties (see “Methods”) (Supplementary Fig. 8). To understand mechanisms underlying the spatial pattern of β-diversity, we obtained the regional-mean value of potential drivers (e.g., paleo-range limit fluctuation and paleoclimatic fluctuation), by averaging them across the target county and its surrounding ones. The paleo-range limit fluctuation, represented by the standard deviation of the high-frequency component of time series of the upper range limit of trees (see Methods), was negatively correlated with β-diversity (P<0.01), with a lower β-diversity found in regions with greater temporal variability in paleo range limits. Paleo-range limit fluctuation (R=−0.53, P<0.01) is more important than both paleo-growing-season temperature (R=−0.33, P<0.01) and precipitation (R=−0.09, P>0.1) fluctuation in explaining the current spatial pattern of β-diversity (Fig. 3a–c). By contrast, β-diversity did not vary significantly across gradients in mean climate and seasonality (Supplementary Fig. 9). Our results suggest the important role of both paleo-range limit and paleo-temperature fluctuations in determining the current spatial pattern of β-diversity of alpine endemic flora across regions. Specifically, both increased paleo-range limit and paleo-temperature fluctuations lead to spatial homogenization of endemic species composition.

To understand why β-diversity reduced as paleo fluctuations intensified, we partitioned β-diversity into two components: dissimilarity due to species replacement (species present at one county are absent at another, but are replaced by other species absent from the first) and dissimilarity due to nestedness (species present at one country are absent at another, but are not replaced by additional species) (Supplementary Fig. 8) (see Methods). The species replacement component predominated in regions with low paleo-range limit fluctuations, while the nestedness component predominated in those with higher fluctuations (Supplementary Fig. 8 and Supplementary Fig. 10). Consistent with β-diversity patterns, the species replacement component decreased with increasing fluctuation, while the nestedness component showed the opposite trend (Fig. 3d, g).

This pattern is expected because gradients in paleo-range limit fluctuation could not only drive species replacement but also species loss and gain. Specifically, decreasing paleo-range limit fluctuations suggest increasingly stable alpine environmental conditions. This increased stability could enhance the capacity of alpine environments to promote niche differentiation and allow diverse assemblages with more specialized habitat preferences19,60–64, thereby increasing β-diversity through species replacement. In contrast, increasing paleo-range limit fluctuation could lead to increased fragmentation (e.g., reduced fragment area and increased isolation) and the contraction of suitable habitats for the survival of endemic alpine species. Increased fragmentation leaves narrowly-ranged endemic alpine species, which typically have small effective population sizes and limited dispersal capacities1,4–7, particularly vulnerable to extirpation and recolonizations, thereby increasing β-diversity through nestedness65–67. Therefore, as paleo fluctuation intensified, the decline in β-diversity resulted from a larger decrease in the species replacement component offsetting an increase in the nestedness component.

In addition, we observed a stronger response of β-diversity to paleo-range limit fluctuation than paleo-temperature fluctuation across regions (Fig. 3a, b). This result is expected because shifts in the upper range limit of trees could induce fragmentation and expansion of suitable habitats for alpine endemism to a larger extent than temperature changes alone. Tree expansion into alpine areas could directly lead to the loss of suitable habitats, and also cause habitat displacement by altering microclimate conditions through increased sheltering effects12,13. Although warming could potentially reduce the availability of suitable habitats for endemic alpine species, topographically diverse environments in alpine areas may provide suitable microhabitats68. These microhabitats allow endemic alpine species to track suitable habitats over short distances69. The presence of these microhabitats could then potentially buffer large-scale atmospheric warming impacts on the survival of alpine endemism. In the coming few decades, projected rates of the upper range limit of trees are several orders of magnitude larger than those occurring over the last millennia. In addition to drastic warming, accelerated tree expansion70 would exacerbate alpine habitat fragmentation and spatial isolation of remaining microhabitats12. Consequently, alpine endemics are not equipped to respond to the dual force from both accelerated tree expansion and drastic warming. This likely would accelerate spatial homogenization of endemic alpine species diversity, threatening the total endemic species pool on the Tibetan Plateau.

Methods
Application of multi-source satellite data in mapping upper range limit of trees
To map the spatial distribution of the upper range limit of trees over the Tibetan Plateau, we first developed a fusion method to map forested extent for the year 2020 (30 m spatial resolution) by integrating multiple-data sources including optical, microwave, and LiDAR data. Forested regions were defined as regions with tree canopy cover larger than 10% because the extracted upper range limit of trees using this value had the highest correlation with in situ observations12.

Generating the spatial distribution of forest extent over the Tibetan Plateau. Existing publicly-accessible datasets on forest extent systematically underestimate the presence of sparse tree cover76, leading to a direct underestimation of the upper range limit of trees. As such, we have constructed a tripartite upscaling method that integrates in situ identification of landcover types, visual interpretations from Google Earth’s high-resolution (<5 m) images, and data from Landsat multi-spectral images, Sentinel-1A microwave data, and GEDI-derived canopy height. These data sources were integrated using a Convolutional Neural Networks (CNN)77 classification algorithm to produce a high-precision map of forest extent for the year 2020 at a 30 m resolution.

We used space-borne LiDAR measurements of canopy height over 39,675,123 valid laser shots covering a resolution of 25 m on the Tibetan Plateau, collected by the GEDI instrument aboard the International Space Station78. We first randomly selected 47,658 sites and visually interpreted their corresponding vegetation classes using high-resolution (<5 m) Google Earth images produced between May and September 2020 and finally identified a total of 26,541 forested sites. In addition, we also surveyed 346 forested sites in the field during the autumn of 2020 and measured their tree heights.

We then built a CNN classification model linking the landcover types (presence of forest) to structural (GEDI-derived tree height), texture79 (Sentinel-1A microwave-derived co-polarization and cross-polarization data), and spectral attributes80 (Normalized Difference Vegetation Index, and seven spectral reflectance bands from Landsat images) during the period from June to September across GEDI sites. We retained 10% of the data (4765 sites) for independent validation of our constructed CNN model (R2=0.89, P<0.01). We then use the remote sensing data (except GEDI data) as predictors in the CNN classification model to generate a high-resolution (30 m) map of forests for the year 2020.

Generating spatial distribution of upper range limit of trees over the Tibetan Plateau. We then used the developed forest extent map as the input for forest edge detection based on the Canny algorithm81, which extracted the forest boundary by finding abrupt changes in the forest image. We removed the lower forest boundary and only maintained the upper forest edges by defining a local adaptive elevation threshold, below which pixels were removed following Wang et al.12. Specifically, the process begins with constructing a histogram of forest edge elevation within a specified window size centered on a geographic point. The choice of window size is based on empirical analysis (ranges from 0.1 to 1°), balancing the inclusivity of upper and lower edges. We then assessed the unimodality of the distribution within the window size by smoothing the histogram using a Savitzky-Golay filter and determining the number of Gaussian functions in the data. If the number of Gaussian functions exceeds one, we found a threshold to separate the upper and lower ones, based on the position that contains a 99% signal in a Gaussian distribution. Subsequently, we eliminated the lower forest edge and iterated through these steps until a unimodal distribution was achieved (i.e., the number of Gaussian functions equals one).

Validating satellite-derived upper range limit of trees. We validate the satellite-derived upper range limit elevation of trees using manually interpreted samples from high-resolution (HR) (<5 m) satellite imagery provided by the Google Earth platform. To create an equal probability sample, the number of sample points is directly related to the number of pixels that have their upper range limit tree observed by the satellite. We generated a geographical 1 km grid over the Tibetan Plateau, and then randomly selected ten validation samples at a 30 m resolution for each 1 km grid to ensure even coverage across the region. Six well-trained workers visually interpreted and cross-validated all samples using Google Earth and the ArcGIS platform. Furthermore, a single quality controller checked all the results. We therefore obtained 602,879 validation samples evenly distributed across the Tibetan Plateau, with 30,746 (5.1%), 510,035 (84.6%), and 62,098 (10.3%) distributed in the northern, southern, and inner regions, respectively. Validation shows that the satellite-derived range limit has a high degree of consistency with manual interpretations from Google Earth images (R2=0.98, Slope=0.98, ME=17 m; Fig. 1b).

Identification of drivers of current spatial patterns in the upper range limit of trees
Not all upper range limits of trees reached the thermal treeline position determined by temperature12,21,22. In reality, the upper range limit of trees is often lower than the thermal treeline position due to environmental stresses12,41 (e.g., critical water shortage) and disturbances12 (e.g., human activities). We first compared the satellite-derived upper range limit of trees to the thermal treeline position, and then evaluated which factors were responsible for the deviation of satellite-derived range limit from thermal treeline positions.

To do so, we first calculated the thermal treeline elevation using growing-season temperature at a spatial resolution of 1 km. The thermal treeline is the cold edge of the fundamental niche of trees21, and was calculated as the isotherm of the mean growing-season air temperature at 2 m above the ground (T_air) of 6.4±0.7 °C37 using the WorldClim dataset and the mean growing-season land surface temperature (T_surface; the radiative temperature of the land derived from infrared radiation emitted from the surface) of 7.6±1.0 °C12 using the Terra Moderate Resolution Imaging Spectroradiometer (MODIS). A global meta-analysis of more than 30 thermal treeline sites from different climate zones showed that thermal treeline positions generally occur at a mean growing-season ground temperature of 6.7 °C at 10 cm depth, with an narrow amplitude of 2.2 °C23,37. Since there was a strong correlation12 between growing-season mean T_surface and ground temperature at a 10 cm depth across these global treeline sites21, the treeline was then estimated to occur at the growing-season mean T_surface of 7.6±1.0 °C. In addition, a meta-analysis of thermal treeline sites at the global scale also showed that the start and end of growing season was defined as date when the ground temperature at 10 cm depth exceeds 3.2 °C23. Such ground-temperature threshold corresponds to MODIS-derived land surface temperature threshold of 0.7 °C. Therefore, the start and end of the growing season were determined as the date when daily T_air rises above 0.9 °C, the date at which it falls below 0.9 °C, and the daily T_surface passes 0.7 °C, respectively. Furthermore, the growing season length must be longer than 94 days37.

We then calculated the deviation of satellite-derived upper range limit of tree elevation from thermal treeline elevation as D_Treeline (the difference between the edges of the realized and the fundamental niche)12,22. Before calculating D_Treeline, we aggregated satellite-derived upper range limit of trees elevation from the original 30 m to 1 km resolution by calculating the mean value in each 1 km grid cell. To understand the drivers of D_Treeline, we compiled several variables (N=19) belonging to the four different categories (climate limitation, disturbance, soil, and topography factors) (Supplementary Table 1), and then used a random forest model to rank the importance of these variables and calculate the partial contribution of each variable on D_Treeline. The entire dataset, including cases where D_Treeline ≥ 0 (N=116,337), was used in our analysis. We generated the partial-contribution plot using the “forestFloor” package in R statistical software (http://cran.r-project.org/), to visualize how each independent variable impacts the prediction while controlling changes in all other variables (Supplementary Fig. 1b to t).

Constructing a climate-based predictive model of upper range limit of trees
In constructing the model, we first used a growing-season temperature threshold to represent the thermal treeline position, and then included several variables including environmental stresses and disturbances, to account for the observed deviation of the satellite-derived upper range limit from the treeline position (D_Treeline). Thus, the modeled upper range limit of trees is the sum of the thermal treeline position and D_Treeline due to environmental stresses and stochastic effects. We separated temperature from other climatic variables when modeling D_Treeline, as there was significant covariation between temperature and other climatic variables, such as indicators of critical water shortages (e.g., vapor pressure deficit and cumulative water deficit), both of which are influenced by temperature.

In modeling D_Treeline, we employed a recursive feature elimination method to identify the most relevant variables12. The error decreased when using the six most important variables in determining the spatial variation of D_Treeline (Supplementary Fig. 11) but decreased by less than 1% when additional variables were included. We then constructed a parsimonious climate-D_Treeline model by selecting the following six variables: spring cumulative climatic water deficit (CWD), spring vapor pressure deficit (VPD), cloud cover, minimum temperature of the coldest month, summer dehydration, and surface curvature (Supplementary Fig. 1 and Supplementary Table 1). The mechanisms for the potential impact of these variables on tree growth and establishment at thermal treeline positions were formulated as follows. For example, the water availability (represented by CWD and VPD) could constrain tree growth and become decisive for seedling establishment, therefore affecting the tree establishment at the thermal treeline position. The CWD, which is calculated as potential evapotranspiration minus precipitation, indicates the quantity of water needed to meet evaporative demand and measures the drought stress on plants, and VPD, as a function of relative humidity and temperature, represents an independent moisture stress on plant photosynthesis. Furthermore, summer dehydration, which is calculated as a product of wind speed, temperature, and relative humidity82, generally occurs during daylight hours under conditions characterized by high temperatures, low relative humidity and strong winds, thereby leading to trees transpiring moisture at a faster rate than their roots take it up from the soil. Such summer dehydration could further aggravate the moisture stress (CWD and VPD) on tree growth and tree establishment at thermal treeline positions. The relationship between cloud cover and tree growth at thermal treeline position is complex. The high cloud cover could reduce the solar radiation, and potentially undermine the plant photosynthesis and tree growth at thermal treeline positions. While the cloudy days are generally associated with the relatively high-water availability and reduced transpiration during the growing season; and this could reduce the radiative cooling and trap heat during the night, potentially allowing trees to survive extreme cold temperatures at the relatively high elevation24. The minimum temperature of the coldest month is used as the freezing stressor that could directly control the absence of trees from the thermal treeline position21. In addition to these macroclimatic drivers, we also included local topographical factors (e.g., surface curvature) as a surrogate of microclimate conditions, because of the potential deviation of microclimate from macroclimate at the high elevation. After including these variables, the parsimonious model achieved a high prediction power (R2=0.84, P<0.01; Supplementary Fig. 2), and is less susceptible to errors in model extrapolation than full models considering all potential drivers. We then applied our parsimonious climate-D_Treeline model to predict changes in D_Treeline at a spatial resolution of ~1 km at 100-year intervals from LGM to the end of this century (mean over the period 2080–2099). All of these predictors from LGM to the end of this century were derived from CHELSA-TraCE21k, WorldClim, and CMIP6 model projections (see “Assembling climate variables”).

Reconstruction and prediction of changes in the upper range limit of trees
We reconstructed spatio-temporal dynamics of the upper range limit of trees spanning from 22 kyr BP to the end of this century at a spatial resolution of 1 km for every 100 years. Specifically, we determined the thermal treeline elevation using the T_air of 6.4 ± 0.7 °C37 and the T_surface of 7.6 ± 1.0 °C12, with the growing season length longer than 94 days37, and then predicted changes in D_Treeline due to regional environmental stresses using the parsimonious D_Treeline model.

Assembling climate variables. Climate data from 22 kyr BP to 2100 CE was compiled from three different data sources: (1) monthly climate data from CHELSA-TraCE21k (Climatologies at High resolution for the Earth’s Land Surface Areas) at a temporal and spatial resolution of 100 years and 1 km, respectively, for the last 22,000 years53; (2) current climatology (2020 CE) from the WorldClim dataset83; and (3) future projections from 9 Earth system models at the end of this century (mean over 2080–2099) participating in CMIP6 under low (SSP1–2.6), intermediate (SSP2–4.5) and high emission scenarios (Supplementary Table 3). Here, SSPx–y represents an integrated scenario of future climate and societal change, with a peak in radiative forcing at y W m−2 before 2100.

We used monthly CHELSA-TraCE21k at a spatial resolution of ~1 km and a temporal resolution of 100 years. This CHELSA-TraCE21k dataset was downscaled from the transient climate simulations (CCSM3 TraCE-21k) from the LGM to present day using the CHELSA V1.2 algorithm53. CHELSA is a high-resolution climate dataset for the time period 1979–201384, which uses statistical downscaling of the ERA-interim reanalysis with a bias correction of monthly precipitation from Global Precipitation Climatology Centre. The downscaled paleoclimate data have included the impact of paleo-orography on temperature and precipitation fields given that topography could have changed drastically due to the ice retreat along poles and high mountain areas through time85. The CHELSA TraCE21k was shown to produce a reasonable representation of high-spatial-resolution temperature and precipitation through time, and simulations using this dataset could be capable of detecting effective LGM species refugia53. Additionally, we analyzed changes in orography since LGM and found that a decrease in elevation is mainly located in the river valleys of the southeast Tibetan Plateau and over Hengduan mountain ranges (Supplementary Fig. 12). By contrast, orography changes within the 1 km buffer around the current range limit of trees are almost negligible.

Note that we used WorldClim data to map the present-day upper range limit of trees, and estimated the paleo-range limit by combining the current climatology from WorldClim with paleo-climatic changes derived from CHELSA. Furthermore, we generated future climate data at a spatial resolution of ~1 km by downscaling and bias-correcting CMIP6 model simulations using WorldClim data. Therefore, discrepancies among different data sources to estimate dynamics in the upper range limit of trees could be well resolved through bias correction.

Validating changes in the upper range limit of trees using fossil pollen data analysis at the millennium scale. We then used fossil pollen assemblages preserved in lake sediments over the Tibetan Plateau to determine if the trees were present in the vicinity of the sampling site throughout the Holocene. By incorporating data from an elevational transect, this approach can be used to delineate the upper range limit of trees and their elevational changes over time. We compiled a list of taxonomically harmonized and temporally standardized fossil pollen data from lake sediments over the Tibetan Plateau86 (Supplementary Fig. 7). We only included fossil pollen abundance data that falls within a 500 km radius of the present forest communities, to minimize the risk of misinterpreting the occurrence of trees due to long-distance transport of exogenous arboreal pollen into lake sediments87. This filtering resulted in a total of six lakes distributed across a range of elevation from 3700 to 4500 m26,35,54–57 that had data for at least 500 years (Supplementary Table 4).

We concentrated on the four most common high-elevational arboreal tree taxa—Betula, Abies, Picea, and Pinus. These tree species have different production and dispersal strategies, which could greatly affect the representation and abundance of tree taxa in fossil pollen records88–90. For example, Betula and Pinus have high pollen production rates, good wind transportability and good preservation characteristics88,91–93. Therefore, even a relatively high percentage at a given elevation might report a false presence. Since there were no available estimates of pollen production for the main tree taxa on the eastern Tibetan Plateau, we then used the reliable threshold of pollen abundance to determine if one of the tree taxa was present in the vicinity of the sampling site. Specifically, the pollen abundance thresholds with different probabilities were determined by performing a logistic regression using pollen abundance (0–100%) as the predictor and taxon classification (0 and 1 represent negative and positive classes, respectively) as the response variable using 2434 modern pollen samples of the 14 key arboreal taxa from China and satellite-based observation and the Vegetation Atlas of China of modern trees taxa distribution89. We used the pollen abundance threshold at a probability of 0.9 to determine the presence and absence of trees. The thresholds for Betula, Abies, Picea, and Pinus are 16%, 5%, 9%, and 37%, respectively89. Upon identifying any one of the four tree taxa, we assumed that the trees were present in the vicinity of the sampling site. Not all pollen grains will be equally well preserved throughout the time due to their differing sensitivities to environmental changes94,95, and thus, some tree taxa with very low pollen production rates and/or poor preservation could be missed. However, the probability of such underrepresentation is unlikely in this study. For example, common tree species including Corylus, Nitraria, Tamarix, Cyperaceae, Poaceae, and Fabaceae, which were widely reported to be greatly under-represented in China based on 898 modern pollen sampling sites and 2115 pollen data96, were not found at the current realized range limit of trees over the Tibetan Plateau.

Quantifying beta-diversity of endemic alpine species and its drivers
We aimed to estimate the role of paleo fluctuations in the upper range limit of trees on the spatial pattern of beta-diversity. We initially obtained geographical coordinates and elevational ranges for endemic alpine species at the county level, using the data from the Flora of China97 and the Flora of Tibet (Supplementary Data 1). The selection of endemic species is based on the Chinese Endemic Species list, and our focus was primarily on alpine species whose uppermost distribution surpassed the current upper range limit of trees. In this way, we have identified 2111 endemic alpine species in over 100 counties of the Tibetan Plateau.

First, we calculated the β-diversity of endemic alpine species using indices from the Sørensen-based multiple-county dissimilarity measures59 (Supplementary Fig. 8a). The Sørensen-based multiple-county dissimilarity indices were calculated using the following Eq. (1),

β_diversity = (Σ a + b) / (2 × [Σ i S_i − S_T] + a + b)

a = Σ i<j min(b_ij, b_ji) and b = Σ i<j max(b_ij, b_ji)

For a given focal county i with N neighbors, S_i is defined as the total number of endemic alpine species in i (alpha diversity from endemic species) and S_T is the total number of endemic alpine species in all neighbors (gamma diversity). b_ij, b_ji are the number of endemic alpine species only present in the county i and j, respectively. We calculated β-diversity for all 100 counties selected.

We further partitioned β-diversity into two separate and antithetic components: dissimilarity due to species replacement and dissimilarity due to nestedness using Eqs. 2 and 359. The species replacement component refers to the replacement of some species with others among counties98, and there is no difference in species richness among counties. On the other hand, the nestedness component reflects differences in species richness, when species-poor assemblages are nested in species-rich assemblages. Different counties would then have different species compositions in the absence of species replacement.

β_species replacement = a / (Σ i S_i − S_T + a)

β_nestedness = (b − a) / (2 × [Σ i S_i − S_T] + a + b) × (Σ i S_i − S_T + a)

The accuracy of the Sørensen index to estimate beta-diversity was recently suggested to be heavily dependent on species prevalence, with high species prevalence generally yielding misleading inferences about species associations99. We then quantified the extent to which species prevalence affected the Sørensen index in estimating β-diversity on the Tibetan Plateau, using a statistical mathematical metric generally used for evaluating pairwise similarities instead of assessing similarities among multiple (>2) species/sites. The prevalence of pairwise species is estimated by calculating the log odds ratio of the conditional occurrence probability (the occurrence rate of pair-wise species in total counties)99. The endemic alpine species on the Tibetan Plateau have a generally low prevalence of 0.04±0.05 (Supplementary Fig. 13a), and thus, it is unlikely that pairwise species would invalidate the interpretations of the traditional Sørensen-based metrics of species co-occurrence. This is confirmed by our analysis showing that there is a good correspondence (positive correlation) between the affinity metric of co-occurrence (Alpha MLE) that is not affected by species prevalence99 and the most popular co-occurrence metrics such as those due to Jaccard (R=0.75, P<0.01), Sørensen-Dice (R=0.75, P<0.01), and Simpson (R=0.83, P<0.01) (Supplementary Fig. 13b–d). We thus confirmed the robustness of using Sørensen-based multiple-county dissimilarity measures to estimate β-diversity.

Second, in terms of paleo-temporal fluctuations in the upper range limit of trees and climate drivers, we only retained the high-frequency component of the paleo time series. Specifically, we decomposed the original time series into a finite set of oscillatory components called Intrinsic Mode Functions (IMFs) and a residual trend component, using Empirical Mode Decomposition method100. IMF1 stands for the first IMF obtained during the decomposition process, which is the dominant oscillatory component extracted from the original time series of each variable. The standard deviation of the IMF1 for each county was used to represent paleo temporal fluctuations in the upper range limit of trees and climatic drivers (growing-season temperature and precipitation). In addition, to assess whether warming-induced upslope tree expansion would potentially lead to habitat loss, we also calculated the alpine land area by summing the surface area above the mean elevation of the reconstructed upper range limit of trees in each county, using a 30-m spatial resolution digital elevation model101. For each mountain defined in the GMBA mountain inventory v2.0 dataset102, we created a 100 km radius buffer around a target upper range limit pixel, and then summed up the surface area of all 30 m resolution pixels with elevations higher than that of the target pixel.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This paper uses multiple data sources to create a highly advanced map that shows changes in three-line elevation over time. The authors improved on Hansen et al.'s (2013) Science approach to mapping forest distribution and focusing on the altitudinal gradient, specifically the treeline, instead of the global distribution pattern, which provides a novel take on the idea of forest distributions. The paper also considers a Late-Pleistocene perspective by reconstructing the treeline location over time and evaluates these estimates using paleoecological and dendrochronological records. This is a commendable feature of the paper.\nAlthough the methods used to map tree lines and evaluate their accuracy are clear, I could not determine the exact model employed to reconstruct the altitudinal range limits of trees during the Late Pleistocene. As a result, I could not assess how they converted the precise maps of existing treeline distribution into a model that can be applied to forecast treeline variations based on factors such as water availability, cloud cover, and temperature values.\nI understand that the paper provides a valuable mapping exercise, but I find the presentation of the ecological side of the study to be unclear. Specifically, I couldn't find an explanation or hypothesis as to why changes in tree line location or their variability would determine the replacement of endemic species over space. This lack of explanation makes me wonder about the point of using turnover in endemic species as the ecological response variable. I believe there could be a better link between endemic diversity and habitat stability.\nOther points:\nL45-47: This formulation read weird: abiotic determines are paces in both limits of the fundamental and realised niche.\nL48-51: A missing point is why log term change in treeline location would affect species composition. Also, why is looking to the LGM relevant in this region?\nL72: Replace “obtained” with “derived” and “from” with “using”.\nL93: Saying “0.6 million” reads like an overly inflating number of observations, in my opinion. Just say “600,000.”\nL106: Please clarify what “the remaining one” is. It is unclear to me.\nL116: Insert “the” before “drivers”.\nL123-125: Sorry, but I don't understand what these percentages refer to. Are these the percentages of cells that agree between the remote sensing and the climate-derived limits in a given region?\nL134: Please add here the explanation of what Dtreeline is measuring (the same one you have in L429-432).\nL145-153: Sorry, but this phrasing does not provide enough details to understand what type of predictive model was built to estimate the upper range limits of trees from climatic drivers. Not even after you look into this point in the methods (L439-448)\nAnother missing detail is how you incorporated the possible impacts of migration, specifically, how you included “time-lags from 10 years to 100 years in changes of range limits (L156-157)”.\nL238-239: While I agree with this point, it is not a sufficient explanation as to why we should expect that variability in the three line locations affects the turnover of endemic species. In this regard, it is unclear why you focused on endemics (Shouldn’t this affect the complete set of species?).\nAlso, why did you choose Turnover as your variable of interest? Also, why did you expect treeline stability to be the two facets of beta diversity?\nThese key points are not fully developed in the text (neither in the introduction nor the results), so the purpose of these analyses from a mechanism/hypothesis-based perspective is unclear.\nFigure 2: Add a colour ramp to indicate the values of the red to blue colours in the main panel.\nFigure 3: What are the y-axis units in panel a (%?)?\nL421-424: Can you please state how did you define these temperature limits? Are they coming from Korner (2021)?\nL439-448: This section lacks many details about what type of model was developed and how it was developed. So, I can't assess whether what was done is correct.\nL470-481: I don't think you did this. CHELSA-TRaCE21 is already at this resolution.\nL540-550: An unclear point is how you summarise the “variability” in treeline location and other environmental variables per county.", "reviewer_2": "In this study, Xu Jinfeng and colleagues examined the temporal fluctuation of the upper range limit of trees during the late Quaternary and used a suite of data and models to predict the drivers of spatial turnover of endemic flora. This is a comprehensive study, aggregating various types of data, methods, and models. The work is commendable. Having reviewed several manuscripts from the same lab and institute, I am consistently impressed by the rigor and scientific insight this lab brings to their research. The authors of this manuscript have done a great job in synthesizing an interesting storyline. Beyond providing valuable facts about treeline and tree limit dynamics over geological timescales, this paper substantially contributes to the understanding of treeline dynamics.\nMAJOR COMMENTS:\n1. Drivers of Discrepancy Between Tree Limit and Thermal Treeline: The authors mention in lines 140-141: “and explained nearly 85.3% of spatial variability in DTreeline across the Tibetan Plateau”. Additional details about this RF model would be valuable, including a comprehensive list of variables used and the number of instances/data points utilized.\nThe manuscript lacks sufficient information on the explained variance. It is important to note that variable importance in Random Forest models does not equate directly to variance explained as it does in linear regression models. While variable importance measures provide insight into the relative significance of each predictor in making accurate predictions, these measures do not sum to 100%. They offer a ranking or score based on the impact of removing or altering a variable on model accuracy, but they do not represent a proportion of variance in the outcome variable explained by the predictors in the model.\nIn comparing the position of the tree limit with the thermal treeline, it was found that 85.2% of the upper limits match the thermal treeline position. Consequently, the RF model was applied to the remaining 14.8% of locations where discrepancies exist between the tree limit and thermal treeline. This analysis identified climatic water deficit and vapor pressure deficit as the primary drivers of discrepancies in the position of the tree limit versus the thermal treeline. Was the RF algorithm also applied to the entire dataset for a more generalizable model? Given that CWD and VPD are posited as the underlying drivers of discrepancies across the study area, they should predict the lack of discrepancy at some values not corresponding with high-discrepancy sites. While RF models accommodate highly non-linear responses, an analysis of the entire dataset could provide a more robust ecological inference beyond what has already been done.\n2. Language and clarity: The manuscript’s presentation and language require enhancements for better clarity. It frequently contains lengthy and convoluted sentences, alongside instances of awkward phrasing. The authors would benefit from the assistance of a native English speaker with a scientific background to refine the language. For instance, the sentence in lines 54-59 is overly long and could be simplified for clarity: “There is evidence from classical paleoecological records and ancient DNA in lake sediments showing that the upper range limit of trees, which were nowadays reported in a range of 3389–4631 m based on limited field reports17-19, could potentially reach up to just 3500 m elevation during the Last Glacial Maximum (LGM) and 400–600 m higher than today in Holocene Climate Maximum over some locations of the Tibetan Plateau17,19,20.”\n3. Beta diversity measurement: In the final section of Results, you estimated beta diversity and then partitioned it to turnover and nestedness. Interestingly, an unusual situation arises when the null distribution of an event is unknown, making it impossible for us to compute the P-value or any reliable quantity characterizing some property of the system (called a statistic or metric). One glaring example of this situation is beta diversity, which has been traditionally measured with popular indices such as Jaccard, Sørensen-Dice, and Simpson. An inherent and a systematic problem with using these indices (see many papers by Ulrich, Gotelli, etc) has been repeatedly reported by ecologists and biogeographers who attempted to fix the problem with various techniques of normalization and standardization.\nThis longstanding issue has been addressed with the introduction of a new metric founded on a solid mathematical basis and a null model developed with well-established assumptions, as detailed in Mainali et al. 2022 and Mainali and Slud 2024. Figure 1E from the Science Advances article conclusively illustrates the inaccuracies of the Sorenson index in estimating beta diversity. The ecological interpretations derived could significantly differ, highlighting the discrepancies between the traditional use of the Sorensen index and a statistically robust metric.\nMainali, K. P., E. Slud, M. C. Singer, and W. F. Fagan. 2022. A better index for analysis of co-occurrence and similarity. Science Advances 8:1–9.\nMainali, K.P. and Slud, E., 2024. Null model analyses are not adequate to summarize strong associations: Rebuttal to Ulrich et al.(2022). Journal of Biogeography, 51(2), pp.263-268.\n4. Beta diversity pattern: Starting line 274, the authors fitted a Generalized Linear Model (GLM) to explain beta diversity patterns. It was found that paleo fluctuations in upper range limit predicts the beta diversity pattern more efficiently than paleo fluctuations in growing season temperature. This finding presents some challenges in interpretation, as the primary mechanism through which growing season temperature could influence beta diversity patterns is by affecting species distribution. In contrast, fluctuations in the upper range limit is a direct predictor of beta diversity patterns. Please elaborate on the rationale and reasoning behind this analysis, particularly in understanding the significance of these predictors.\nMINOR COMMENTS:\n• Line 45-46: The statement that \"Shifts in the upper range limit of trees, whose fundamental niche is determined by seasonal mean temperature\" requires substantiation. If references 12-15 support this claim, it would be appropriate to cite these immediately following the statement for stronger evidence.\n• Line 49: You mean \"spatio-temporal map\"?\n• Line 53: Please specify the timeframe of the Last Glacial Maximum (LGM).\n• Line 72 and related paragraph: The term “upper range limit of tree dynamics” comes across as unclear. Can you rephrase it? Additionally, your statement “using the relationship between satellite-derived upper range limit of trees and climate in modern observations, under the common assumption that this relationship has remained constant through time,” needs to be further elaborated and justified.\n• Line 119-120: The reference to Paulsen and Korner in defense of your methodology prompts a question about the spatial grain size of your isotherm. If your unit of analysis (e.g, square grid cell) is too large, then a better match between expected treeline and the observed one is likely due to methodological artifact.\n• Line 127-128: The assertion that \"The upper range limits of trees below this thermal treeline position are considered to be a result of land use\" requires further explanation and citation.\n• Line 133-134: The phrase “the deviation of the upper range limit of trees from treeline position” seems to be missing a reference to the \"thermal treeline position.\" Additionally, clarifying whether this deviation is measured as a surface distance or an elevational gradient is important for interpretation.\n• Line 145-150: Can you add here or in the Methods the total number of variables included from the four categories? What metric was used to evaluate model performance, especially given the claim of over 80% spatial variation explained?\n• Line 154-160: Is \"improve the explanatory power\" meant to imply \"maximize\"? In the rest of the paragraph, the authors presented the results of space for time analysis. I have trouble understanding how a 50-year lag was evidenced from the space-for-time results.\n• Line 170: The first definition of LGM should appear at its initial mention in both the abstract and the main text.\n• Line 183 and related paragraph: I was delighted to read the discussion on the four-phase pattern.\n• Line 206-234: The utilization of fossil pollen abundance records from six lakes to support your argument is commendable for its methodological diversity and the range of approaches and data types employed to construct a narrative on tree limit dynamics in the region. However, it would be beneficial to discuss potential biases in your methodology, particularly given that different tree species present in the area at various times may not have pollen that is equally well-preserved. A discussion on how your approach mitigates the risk of underrepresentation of certain species in pollen fossils due to differential preservation would strengthen the reliability of your findings.\n• Section \"Reconstruction and prediction of changes in the upper range limit of trees\" (Line 450): The methodology for reconstructing spatiotemporal dynamics using climate data from CHELSA-TraCE21k and CMIP6 future projections is outlined. How did you ensure the robustness of these projections, particularly given the potential for significant changes in topography and climate conditions over the considered time period?", "reviewer_3": "This paper presents important research regarding drivers of the upper range limit of trees, as well as beta-diversity within a global biodiversity hotspot. It produces important findings both for our understanding of paleo dynamics at treeline but also for understanding future treeline shifts and implications for biodiversity in alpine mountain environments. It is most immediately important to understanding treeline dynamics of the Mongolian Plateau, but findings may also be generalizable to alpine mountain regions globally. It produces important findings for ecologists, earth systems scientists, biogeographers, and wildland managers. The conclusions and claims are well-supported, and I see no flaws in analysis nor interpretation. It is methodologically sound. The methods provide sufficient detail and code for the work to be reproduced. With minor revisions I approve its publication.\nKey results – this paper examines spatiotemporal drivers of the upper range limit of trees, and models the upper range limit over the last several thousand years, for the Tibetan Plateau. It also finds that treeline fluctuation is an important predictor of beta-diversity – areas with more fluctuation have lower levels of beta-diversity, and areas with less fluctuation have higher levels.\nValidity – the authors employ robust methods in their examination of treeline shift and its effects on beta-diversity over the Holocene. Their methods include remote sensing, modeling, machine-learning, in situ validation, paleo-records, and so on.\nSignificance – this paper has great significance not only to mountain regions of the Mongolian Plateau, but to mountain regions globally. It provides important insight into the effects of treeline stability on beta-diversity in alpine environments.\nData and methodology – Their methods are multidisciplinary and rigorous.\nAnalytical approach – the analytical approach is novel, valid, and comprehensive. I have only a few questions regarding their methodology (see major/minor comments below).\nSuggested improvements – I have grouped my suggestions into “major comments” and “minor comments” below. They pertain mainly to organization and language.\nClarity and context – overall the manuscript is clear, but I do make several suggestions regarding improving the flow of the paper, fixing grammatical issues, maintaining consistent tense (mainly past tense), and so on.\nReferences – yes, I believe the manuscript appropriately references previous literature.\nYour expertise – I teach mountain geography, GIS, and environmental modeling courses at my institution. I am a biogeographer by training. My research is not in treeline dynamics, but I do use many modeling techniques described in the manuscript, and I do have fair understanding of treeline and of paleo records used throughout this analysis. I know remote sensing concepts. I understand most/all of the modeling decisions and techniques made and their rationale.\nMajor comments\nIn the introduction, please more formally define which kind of “treeline” you are reconstructing. The terms timberline, treeline, forestline, and scrubline are all terms used for describing different aspects of the forest to alpine tundra transition. Please include a definition somewhere of which you are reconstructing. I am guessing you are reconstructing treeline, the upper limit of upright tree growth. Scrubline would describe the upper limit of prostrate trees.\nIn the results, please be more explicit regarding which timeframe your results are applicable to. For example, in line 100-115, you are referring to current upper range limit of trees, correct? In such locations in the paper, please add words like “present” or “current” upper range limit of trees to remind the reader.\nFor your model of upper range of trees, why not make temperature a variable and directly model upper range, instead of modeling the deviation from the temperature-based estimated limits? I think there could be a better explanation of the decision-making here – perhaps firmer statements regarding how temperature is a major driver of treeline.\nOverall, I think the authors should read over the results section and find ways to move some of the material into the methods section. I think the results should be more results-driven. In addition, some of the results stray into discussion, so re-read the results and find places to move material into the discussion. Overall the results section needs to be more focused on plain presentation of results.\nIf you have multiple models in this analysis, then consider naming and abbreviating them in the results, rather than working to describe them every time they are mentioned. E.g. “the current forest distribution model” or the “presence/absence of trees model” etc.\nIn the minor comments below, I often suggest adding a few clue-words throughout the paper to describe when results are derived from models rather than some other source. Read the results over and please sprinkle these hints to remind the reader that models are producing many of these quantitative results (e.g. treeline elevation, extent of alpine tundra)\nMinor comments\n25 “position” – perhaps simply say “elevation”\n29-30 “fluctuation in range limit of trees” – are you writing about fluctuation over time or over space? I think you should add a few words in the latter portion of the abstract to be clear as to whether you are talking about fluctuation from time-to-time or from location-to-location\n39 should say “patterns” (plural)\n42 could cut “While, in addition to past climate change,”\n42-43 “endemic species with a restricted range” sound like species found in alpine habitats. If this is mostly true, perhaps you should add “alpine” to this part.\n49-51 reword. “we required a spatial-temporally map” is confusing. Change to “it is essential to map the upper range limit of trees over time” or something similar.\n59 use “however,” instead of “but”\n60-61 “to create a spatial-temporally map” is awkward – say “to map reconstructions over time of the upper range limit…”\n62 cut “while”\n66 “provide” not plural “provides”\n67 “tree existence” – simply say “trees”\n72 use “developed” or “estimated” instead of writing “obtained” and also write “limits” (plural)\n100 – “Our spatial distribution” … change to “Our analysis”101-104 this part is somewhat confusing to me. What does it mean that “upper trees are currently concentrated in the northeastern and southeastern margins of the Tibetan Plateau”? Please reword here. Are you saying this is where treeline is highest within your study area?\n106 – “The remaining one” – what are you referring to? Add a few words of clarification\n106 “scatteredly distributed” – awkward language – rewrite here – perhaps simply write “scattered throughout the headwater region”\n108 “tree” not “trees” I think you mean\n111 – “mean elevation” “of treeline” you mean to say?\n116-142 – consider placing some of this material in the methods section. This paragraph should be more results-oriented.\n145 “we take the knowledge from the”… rewrite as “we applied the”\n140 versus 148 – there are different measures of deviance explained. Be clearer as to the differences between them. Is one calculated on training data, and the other is calculated on the independent validation dataset?\n173 should you add “according to our model” to the beginning of this sentence? Where are these results coming from exactly – paleorecords or a model? Add a few words to clarify.\n183 again, in this paragraph, remind the reader that these are model-driven results, not based on fossils, pollen, etc., by adding a few words\n218 should “grids” be instead written “grid cells”? You are extracting pollen-derived data at lake sites, which correspond to specific grid cells.\n239-243 remind the reader that these are model-driven results? E.g. “modeled minimum extent”, etc.\n243 cut “While”\n238-252 some of this material sounds like introduction, some sounds like methods, and some sounds like results. Do consider how you might edit/revise this paragraph to re-distribute and better organize the material\n300-307 – I like the ideas here. Yes, perhaps treeline position and dynamics is more crucial to understanding biodiversity patterns in alpine tundras than climate directly. Has anyone else made this point before? I would suggest reviewing the literature some more to see if your statements here support previous research.\n306 cut “most promising”\n313 put “such as sedimentary ancient DNA” in parentheses, I think\n337 “four distinct four-phases” …is this a typo? Rewrite as “four distinct phases”\n363 add hyphen – “publicly-accessible”\n407 “ones” – replace with a more descriptive word\n411-412 who exactly were these interpreters? Students? Labor from Amazon Mechanical Turk?\n412-415 why so much validation? Could you spend some time explaining your rationale?\n421-427 can you provide more justification for this method? Point to previous literature? I’m curious as to why you can use temperature in this way, rather than using temperature in a model of upper range limit of trees.\n442-445 here and elsewhere, explain the sources of the data. Or, if you don’t explain the sources here, mention that tables in the supplementary material explain where the datasets are from, their original spatial resolution, etc. Just ensure that somewhere in the paper you direct readers to sources of the variables and the characteristics of the data.\n484-485 “the life-form tree was” – just rewrite as “trees were”\n495 add “m” after “4500” for meters\n496 rewrite as “We concentrated on”\n499 rewrite as “were determined” (past tense)\n505-506 the amount for Pinus seems high, and Betula seems somewhat high also, as if those trees should still be “present” at these pollen levels. Is there are reason why pollen levels could be this high yet the tree species could be “absent”? It’d be good to cite some literature explaining perhaps that pollen of some taxa is transported longer distances, and/or that these taxa produce more pollen than others – if that’s the case.\n507 again simplify “the life-form tree was” to “trees were”\n511 “We aimed to estimate” (past tense)\n516 “was” not “is”\n511-518 and elsewhere – make sure you are writing mainly in past tense\n551-558 please be specific about where the alpine land area data are coming from. Is this initially based on your model of forest cover (from earlier in the methods section)?\nIs it at all worth mentioning somewhere in the paper that treeline position changes albedo? I imagine Arctic treeline latitude is more crucial to altering Earth’s albedo, but perhaps you could somewhere expand upon the significance of your work by stating that shifting mountain treelines could affect albedo (at least locally) – it’s up to you.\nSupplementary Fig. 4 – edit this caption because there are some awkward sentences and misspellings.\nSupplementary Table 1 – please add the original spatial resolution of each variable. I think it is important to be clear that some variables had finer or coarser resolution.\nSupplementary Table 5 – please explain the source(s) of alpine flora and their elevational ranges", "reviewer_author_response": "[Comment 1] This paper uses multiple data sources to create a highly advanced map\nthat shows changes in three-line elevation over time. The authors improved on\nHansen et al.'s (2013) Science approach to mapping forest distribution and focusing\non the altitudinal gradient, specifically the treeline, instead of the global distribution\npattern, which provides a novel take on the idea of forest distributions. The paper also\nconsiders a Late-Pleistocene perspective by reconstructing the treeline location over\ntime and evaluates these estimates using paleoecological and dendrochronological\nrecords. This is a commendable feature of the paper.\nAlthough the methods used to map tree lines and evaluate their accuracy are clear, I\ncould not determine the exact model employed to reconstruct the altitudinal range\nlimits of trees during the Late Pleistocene. As a result, I could not assess how they\nconverted the precise maps of existing treeline distribution into a model that can be\napplied to forecast treeline variations based on factors such as water availability,\ncloud cover, and temperature values.\nI understand that the paper provides a valuable mapping exercise, but I find the\npresentation of the ecological side of the study to be unclear. Specifically, I couldn't\nfind an explanation or hypothesis as to why changes in tree line location or their\nvariability would determine the replacement of endemic species over space. This lack\nof explanation makes me wonder about the point of using turnover in endemic species\nas the ecological response variable. I believe there could be a better link between\n1\nendemic diversity and habitat stability.\n[Response] Thank you for your positive and constructive comments. In order to well\nresolve the reviewer’s comments, we have made substantial revisions to the main\nmanuscript. First, we have rewritten the “Results and discussions” and “Methods”\nsections to clarify how we incorporated our mechanistic understandings of the current\nrange limit of trees into a climate-driven model reconstructing the upper range limit of\ntrees since the Last Glacial Maximum (see detailed responses to Comment 8 by\nReviewer #1). Second, we have included detailed ecological explanations on how the\nstability in upper range limit of trees would determine the β-diversity and its two\ncomponents (see detailed responses to Comment 9 by Reviewer #1). To further\nstrengthen the ecological aspect of this study, we have generated a null expectation for\nβ-diversity and assessed the extent to which the observed β-diversity pattern along\ngradients of paleo-range limit fluctuation is driven by ecological processes versus\nneutral assembly processes. Third, we have rewritten the first paragraph of the\nintroduction to emphasize the significance of our focus on endemic alpine species and\ntheir β-diversity, and have included a discussion on how our findings can inform\nconservation planning for endemic alpine species on the Tibetan Plateau. We do hope\nthat our responses could well address the reviewer’s comments.\nSpecific Comments\n[Comment 2] L45-47: This formulation read weird: abiotic determines are paces in\nboth limits of the fundamental and realised niche.\nspecies composition. Also, why is looking to the LGM relevant in this region?\n[Response] Following your constructive suggestion, we have rephased the definition\nof upper range limit of trees as “The upper range limit of trees, which is defined as a\nconnection between the uppermost group of upright trees of a certain height12,19, could\ntheoretically reach their physiological limits set by temperature (referred to as the\ntreeline position)19,20 but it is often confined to lower elevations due to environmental\nstresses and stochastic effects such as human disturbance12,19-22.” (Lines 54-58 on\nPages 3-4)\nWe have added the following text to explain why changes in upper range limit of\ntrees would affect alpine species composition: “It is known that local diversity of\nendemic alpine species, which typically have small effective population sizes and\nlimited dispersal capacities6,7, is threatened by climate change8,9 and the upward\nmovement of trees10-13. The upward migration of tree species threatens alpine\necosystems through habitat loss and fragmentation, as well as alterations in\nmicroclimate (e.g., increased sheltering effect).” (Lines 39-43 on Page 3)\nLast, we have added the following text to explain why we considered the period\nsince LGM in the analysis: “The period since the LGM was used because its\nsignificant paleoclimatic changes led to notable vegetation and floristic shifts, as\nevidenced by fossil records preserved in sediments29-32, potentially driving legacies in\nmodern diversity patterns.” (Lines 80-82 on Page 5)\n3\n[Comment 3] L72: Replace “obtained” with “derived” and “from” with “using”.\nL93: Saying “0.6 million” reads like an overly inflating number of observations, in\nmy opinion. Just say “600,000.”\n[Response] Done as suggested.\n[Comment 4] L106: Please clarify what “the remaining one” is. It is unclear to me.\n[Response] We have rephrased it as “The remaining 6.0% (N = 721,452) of pixels\ndepicting the upper range limit of trees” (Lines 109-110 on Page 6)\n[Comment 5] L116: Insert “the” before “drivers”.\n[Response] Done as suggested.\n[Comment 6] L123-125: Sorry, but I don't understand what these percentages refer\nto. Are these the percentages of cells that agree between the remote sensing and the\nclimate-derived limits in a given region?\n[Response] Here the percentage refers to the proportion of pixels with satellite-\nderived upper range limit of trees reaching climate-driven (or thermal) treeline\nposition over the Tibetan Plateau. In the revised manuscript, we have rephrased the\nsentence as “Our results show that nearly 85.2% of pixels have agreement between\nthe upper range limit of tree and the thermal treeline19,36 positions over the Tibetan\nPlateau, suggesting a predominant role of temperature in determining upper range\nlimits of trees. In the southern and inner Tibetan Plateau, nearly 87.5% and 51.8% of\npixels have their upper range limit of trees reaching thermal treeline positions,\nrespectively, while only 27.9% of pixels in the northern Tibetan Plateau reached their\n4\nthermal positions.” (Lines 124–130 on Page 7)\n[Comment 7] L134: Please add here the explanation of what Dtreeline is measuring\n(the same one you have in L429-432).\n[Response] Following your kind suggestion, we have added the following sentence\ninto the revised manuscript: “Trees do not always reach their thermal treeline\npositions approximated by an isotherm of mean growing-season temperature18-21 (see\nMethods), and they often remain confined to lower elevations due to human\nactivities12,19,38 (such as grazing, logging and/or burning) or other disturbances, or\ncritical water shortage.” (Lines 120–124 on Page 7)\nFurthermore, we have mentioned how D was calculated: “To do so, we\nTreeline\ncalculated the elevational difference between the satellite-derived upper range limit of\ntrees and thermal treeline12,19,20, denoted as D\nTreeline\n” (Lines 132–134 on Page 7)\n[Comment 8] L145-153: Sorry, but this phrasing does not provide enough details to\nunderstand what type of predictive model was built to estimate the upper range limits\nof trees from climatic drivers. Not even after you look into this point in the methods\n(L439-448) Another missing detail is how you incorporated the possible impacts of\nmigration, specifically, how you included “time-lags from 10 years to 100 years in\nchanges of range limits (L156-157)”.\n[Response] Sorry for our unclear expression. We have rewritten the previous two\nparagraphs as the following text: “We developed a predictive model to reconstruct and\nproject changes in the upper range limit of trees. As there are time lags in seed\n5\ndispersal and tree establishment, the shift in the upper range limit of trees may not\noccur in synchrony with climate change19,33,41,42. We included both climate change\nand migration lag in our predictive model. We used growing season temperature along\nwith climatic stressors (especially spring CWD and VPD) to parameterize the effect\nof climate change on temporal shifts in the upper range limit of trees, and developed a\nclimate-based predictive model using a training dataset (N = 116,337) (see Methods).\nOur climate-driven predictive model could explain nearly 84% of spatial variation of\nthe upper range limit of trees using an independent validation dataset (N = 11,634)\n(Supplementary Fig. 2).\nWe then incorporated the impact of migration lag on shifts in the upper range\nlimit of trees into the climate-driven predictive model. To determine the optimal lag,\nwe tested time lags from 10 years to 100 years in our climate-driven predictive model,\nby minimizing the difference between simulated and in situ shifts in the upper range\nlimit of trees over the past century across 20 sites of the Tibetan Plateau. These range\nlimit elevation dynamics were reconstructed from tree ring data using the\ndendrochronological method38,43-50 (Supplementary Fig. 3 and Supplementary Table\n2). Among all time lags tested, incorporating a 50-year time lag into the climate-\ndriven predictive model maximizes its explanatory power in representing in situ shifts\nin the upper range limit of trees over the last century. Specifically, the predictive\nmodel explains nearly 70% of the spatial variation in range limit shifts across 20 sites\non the Tibetan Plateau (R2 = 0.68 P < 0.01) (Supplementary Fig. 4). Furthermore,\nafter including the time lag, the trees across all sites are predicted to shift upward at a\n6\nrate of 20 ± 21 m decade-1, which becomes closer to in situ shifts (18 ± 23 m decades-\n1). This confirms the robustness of our constructed model in the prediction of the\nupper range limit of trees at the century timescale (Supplementary Fig. 3b).” (Lines\n145-170 on Pages 8-9)\nIn addition, we have rephrased the previous method as follows. “In constructing\nthe model, we first used a growing-season temperature threshold to represent the\nthermal treeline position, and then included several variables including environmental\nstresses and disturbances, to account for the observed deviation of the satellite-derived\nupper range limit from the treeline position (D ). Thus, the modeled upper range\nTreeline\nlimit of trees is the sum of the thermal treeline position and D due to\nTreeline\nenvironmental stresses and stochastic effects. We separated temperature from other\nclimatic variables when modeling D , as there was significant covariation\nTreeline\nbetween temperature and other climatic variables, such as indicators of critical water\nshortages (e.g., vapor pressure deficit and cumulative water deficit), both of which are\ninfluenced by temperature.\nIn modeling D , we employed a recursive feature elimination method to\nTreeline\nidentify the most relevant variables12. The error decreases when using the six most\nimportant variables in determining the spatial variation of D (Supplementary\nTreeline\nFig. 12), but decreased less than 1% when additional variables were included. We\nthen constructed a parsimonious climate-D model by selecting the following six\nTreeline\nvariables: spring cumulative climatic water deficit (CWD), spring vapor pressure\ndeficit (VPD), cloud cover, minimum temperature of the coldest month, summer\n7\ndehydration, and surface curvature (Supplementary Fig. 1 and Supplementary Table\n1). This parsimonious model achieved a high prediction power (R2 = 0.84, P < 0.01;\nSupplementary Fig. 2), and is less susceptible to errors in model extrapolation than\nfull models considering all potential drivers.” (Lines 504-524 on Pages 28-29)\nFig. R1 (also shown as Supplementary Fig. 2 in revised manuscript) Validation of\nthe parsimonious climate-D model over the Tibetan Plateau. 90% of the\nTreeline\nD at a 1 km spatial resolution were employed in parsimonious climate-D ,\nTreeline Treeline\nwhile the remaining 10% were reserved for model validation. The dotted black is 1:1\nline, and the continuous black line is regression slope between observed and predicted\nD . R2 and ME indicate the coefficient of determination and the mean error,\nTreeline\nrespectively.\ndetermination between in situ shifts in realized range limit of trees and model\nsimulations using different time lags ranging from 10 years to 100 years across\nall 20 sites on the Tibetan Plateau.\nFig. R3 (also shown as Supplementary Fig. 12 in revised manuscript) Changes in\nmodel error with the number of variables used to construct the climate-D\nTreeline\nmodel. Model error fluctuated and decreased less than 1% when more than 6\nvariables were included. Only the 6 most important variables were used to construct a\nparsimonious climate-D model.\nTreeline\n[Comment 9] L238-239: While I agree with this point, it is not a sufficient\nexplanation as to why we should expect that variability in the three line locations\n9\naffects the turnover of endemic species. In this regard, it is unclear why you focused\non endemics (Shouldn’t this affect the complete set of species?). Also, why did you\nchoose Turnover as your variable of interest? Also, why did you expect treeline\nstability to be the two facets of beta diversity? These key points are not fully\ndeveloped in the text (neither in the introduction nor the results), so the purpose of\nthese analyses from a mechanism/hypothesis-based perspective is unclear.\n[Response] In order to well resolve the reviewer’s constructive comments, we\nhave rewritten the main paragraphs related to the β-diversity. First, we have rewritten\nthe first paragraph to explain our focus on endemic alpine species and their β-\ndiversity, as well as to propose potential mechanisms linking shifts in the upper range\nlimit of trees to current diversity patterns: “The high mountain ranges of the Tibetan\nPlateau are global biodiversity hotspots with a rich endemic diversity in alpine flora1-\n3, which constitute a high proportion of regional diversity4,5. It is known that local\ndiversity of endemic alpine species, which typically have small effective population\nsizes and limited dispersal capacities6,7, is threatened by climate change8,9 and the\nupward movement of trees10-13. The upward migration of tree species threatens alpine\necosystems through habitat loss and fragmentation, as well as alterations in\nmicroclimate (e.g., increased sheltering effect). To understand how best to conserve\nendemic alpine species diversity across spatial scales, it is imperative to understand\nthe driving mechanisms of beta-diversity, as it is the component of regional gamma-\ndiversity that accumulates as a result of compositional differences between local\nendemic alpine species diversity14-18. The pattern of contemporary diversity is\n10\ngenerally interpreted from past climate regimes and their oscillations, since\npaleoclimatic changes govern processes of speciation, extinction and biogeographic\ndispersal. However, it remains unknown how paleo changes in the upper range limit\nof trees affect the current beta-diversity pattern. This knowledge could provide a\nframework for understanding the impacts of rapid tree expansion with the goal of\nconserving alpine endemism in a rapidly warming world.” (Lines 37-53 on Page 3)\nSecond, we have generated a null expectation for β-diversity to show that the\nobserved β-diversity pattern is mainly driven by ecological processes, and confirmed\nthat paleo fluctuations in both temperature and upper range limit of trees served as\nstrong environmental filters structuring the current β-diversity pattern. We then added\nthe detailed reasoning to resolve the question of why stability in the upper range limit\nof trees is expected to drive the two components of β-diversity. Specifically, we have\nadded the following paragraphs into the main text: “To assess the extent to which β-\ndiversity was driven by ecological processes resulting from these paleo fluctuations,\nwe removed the amount of β-diversity attributed solely to neutral assembly processes.\nSpecifically, we generated a null expectation for β-diversity under a random assembly\nprocess that controls for alpha diversity and sizes of regional species pools, and then\nsubtracted observed β-diversity from this null expectation (hereafter as β-deviation)\n(see Methods). We show that β-deviation was significantly correlated with both paleo-\nrange limit (R = -0.58, P < 0.01) and paleo-temperature fluctuations (R = -0.30, P <\n0.01) (Supplementary Fig. 10), confirming that paleo fluctuations serve as strong\nenvironmental filters structuring the current β-diversity pattern.” (Lines 262-271 on\n“To understand why β-diversity reduced as paleo fluctuations intensified, we\npartitioned β-diversity into two components: dissimilarity due to spatial turnover\n(species present at one county are absent at another, but are replaced by other species\nabsent from the first) and dissimilarity due to nestedness (species present at one\ncountry are absent at another, but are not replaced by additional species)\n(Supplementary Fig. 8) (see Methods). The spatial turnover component predominated\nin counties with low paleo-range limit fluctuations, while the nestedness component\npredominated in those with higher fluctuations (Supplementary Fig. 8 and\nSupplementary Fig. 11). Consistent with β-diversity patterns, the spatial turnover\ncomponent decreased with increasing fluctuation, while the nestedness component\nshowed the opposite trend (Fig 3d and g).\nThis pattern is expected because gradients in paleo-range limit fluctuation could\nnot only drive species replacement but also species loss and gain. Specifically,\ndecreasing paleo-range limit fluctuations suggest increasingly stable alpine\nenvironmental conditions. This increased stability could enhance the capacity of\nalpine environments to promote niche differentiation and allow diverse assemblages\nwith more specialized habitat preferences29,59-63, thereby increasing β-diversity\nthrough spatial turnover. In contrast, increasing paleo-range limit fluctuation could\nlead to increased fragmentation (e.g., reduced fragment area and increased isolation)\nand the contraction of suitable habitats for the survival of endemic alpine species.\nIncreased fragmentation leaves narrowly-ranged endemic alpine species, which\n12\ntypically have small effective population sizes and limited dispersal capacities1,4-7,\nparticularly vulnerable to local extirpation and recolonizations, thereby increasing β-\ndiversity through nestedness64-66. Therefore, as paleo fluctuation intensified, the\ndecline in β-diversity resulted from a larger decrease in the spatial turnover\ncomponent offsetting an increase in the nestedness component.” (Lines 272-297 on\nPages 14-15)\nThird, focusing research on β-diversity and its generation mechanisms could also\nhelp provide guidelines for designing conservation networks for endemic alpine\nspecies. We have therefore added a discussion on how our findings could assist\nconservation planning for endemic alpine species on the Tibetan Plateau: “The\nenhanced paleo-range limit fluctuation, coupled with paleo-temperature fluctuation,\ncould shift the pattern of dissimilarity derived from spatial replacement to a pattern of\ndissimilarity derived from nestedness, which has important implications for\nconservation planning. Specifically, prioritizing conservation policies should target\nmultiple counties in regions with low fluctuation, and the richest counties in high-\nfluctuated region, respectively. These fluctuations are expected to intensify in a\nwarmer world, potentially leading to the homogenization of endemic alpine species\ndiversity across counties, thereby threatening regional endemic species pool. Thus, in\nregions with high paleo fluctuations, we should also design corridors with optimal\nlocations and orientations to enhance alpine habitat connectivity and facilitate the\nmixing of endemic species among counties. These corridors may initially reduce β-\ndiversity in these regions; however, they would decrease the likelihood of endemic\n13\nspecies extinctions and help maintain the regional endemic species pool in the long\nterm.” (Lines 325-339 on Pages 16-17)\nLast but not the least, here we focused on alpine species as those with their main\ndistribution above the upper range limit of trees, and for which a reduction in\nabundance is noticeable just below that limit. The endemic alpine species typically\nhave small effective population sizes and limited dispersal capacities. Here we\nfocused on endemic alpine species primarily because they are urgently threatened by\ndrastic warming and uphill movement of trees. Furthermore, the high mountain ranges\nof the Tibetan Plateau have a rich endemic diversity in alpine flora that constitutes an\nimportant proportion of regional biodiversity. Therefore, there is an urgent need to\ndecipher main processes underlying current endemic alpine species diversity patterns,\nwhich could also yield more accurate projections of their changes in a rapidly\nwarming world. In addition, the important role of paleo-range limit fluctuation in\nregulating current β-diversity pattern is likely to apply to these non-endemic alpine\nspecies as well. Unfortunately, this hypothesis could not be tested here due to a lack\nof data on alpine non-endemic species.\n[Comment 10] Figure 2: Add a colour ramp to indicate the values of the red to blue\ncolours in the main panel.\n[Response] Done as suggested. We added a color bar to represent the four distinct\nphases on the top of main panel. In addition, we improved the coherence of the color\ngradient.\nPaleoreconstruction of upper range limit of trees since 22 kyr to 1950 CE using\nCHELSA TraCE-21k (left) and projection to the end of this century (right) using\nCMIP6 multi-model ensemble mean under the SSP1-2.6 (yellow), SSP2–4.5 (blue)\nand SSP5-8.5 (red) scenario, respectively. The solid line represents the mean elevation\nof the upper range limit of trees on the left, with the dotted line and shaded regions\nrepresent the 5th, 25th, 75th, and 95th percentiles. The gradients of red and blue\nrepresent the four distinct phases. The gray box on the right represents satellite-based\nobservation. The whiskers and boxes represent the 5th, 25th, 75th, and 95th\npercentiles and the white lines show mean values.\n[Comment 13] Figure 3: What are the y-axis units in panel a (%?)?\n[Response] Many thanks for pointing out this error. The y-axis unit in panel is percent\n(%), that represents relative importance of predictor variables (% of R2) in generalized\nlinear mixed model. In the revised manuscript, we have removed this graph.\n15\n[Comment 14] L421-424: Can you please state how did you define these temperature\nlimits? Are they coming from Korner (2021)?\n[Response] Here we used the two temperature thresholds to define the thermal\ntreeline: the mean growing-season air temperature at 2 m above the ground (T ) of\nair\n6.4 ± 0.7 °C (Paulsen & Körner 2014) and the mean growing-season land surface\ntemperature (T ; the radiative temperature of the land derived from infrared\nsurface\nradiation emitted from the surface) of 7.6 ± 1.0 °C (Wang et al., 2022). The first\nthreshold is directly taken from Paulsen and Körner (2014), while the second one is\nfrom our recent study (Wang et al., 2022). Specifically, Paulsen and Körner (2014)\nidentified 376 global treeline sites from satellite images, and then developed a\nclimate-based modeling approach to show that the thermal treeline would occur at a\nT of 6.4 °C by minimizing the model-data misfit in treeline elevation. A global\nair\nmeta-analysis of more than 30 sites from different climate zones showed that the\nthermal treeline positions generally occurred at a mean growing-season ground\ntemperature of 6.7 °C at 10 cm depth, with a narrow amplitude of 2.2 °C (Körner &\nPaulsen 2004). Since there was a strong correlation between growing-season mean\nT and ground temperature at a 10 cm depth across these global treeline sites\nsurface\n(Wang et al., 2022), the treeline was then estimated to occur at the growing-season\nmean T of 7.6 ± 1.0 °C. All of this information has been added in the revised\nsurface\nmanuscript.\nReferences:\nKörner, C. & Paulsen, J. A world-wide study of high altitude treeline temperatures. J.\nPaulsen, J. & Körner, C. A climate-based model to predict potential treeline position\naround the globe. Alp. Bot. 124, 1-12 (2014).\nWang, X. Y. et al. Enhanced habitat loss of the Himalayan endemic flora driven by\nwarming-forced upslope tree expansion. Nat. Ecol. Evol. 6, 890–899 (2022).\n[Comment 15] L439-448: This section lacks many details about what type of model\nwas developed and how it was developed. So, I can't assess whether what was done is\ncorrect.\n[Response] Sorry for the unclear expression in climate-D predictive model. We\nTreeline\nhave rewritten this section (Please see detailed responses to Comment 8 by Reviewer\n#1).\n[Comment 16] L470-481: I don't think you did this. CHELSA-TRaCE21 is already at\nthis resolution.\n[Response] In order not to mislead the readers, we have rephrased these sentences as\nfollows: “We used monthly CHELSA-TraCE21k at a spatial resolution at ~1 km and a\ntemporal resolution of 100 years. This CHELSA-TraCE21k dataset was downscaled\nfrom the transient climate simulations (CCSM3 TraCE-21k) from the LGM to present\nday using the CHELSA V1.2 algorithm51.” (Lines 544-547 on Pages 29-30)\n[Comment 17] L540-550: An unclear point is how you summarise the “variability”\n17\nin treeline location and other environmental variables per county.\n[Response] Sorry for the unclear description. We added a detailed description in the\nmethod section: “Third, in terms of paleo temporal fluctuations in the upper range\nlimit of trees and climate drivers, we only retained the high-frequency component of\nthe paleo time series. Specifically, we decomposed the original time series into a finite\nset of oscillatory components called Intrinsic Mode Functions (IMFs) and a residual\ntrend component, using Empirical Mode Decomposition (EMD) method99. IMF1\nstands for the first IMF obtained during the decomposition process, which is the\ndominant oscillatory component extracted from the original time series of each\nvariable. The standard deviation of the IMF1 for each county was used to represent\npaleo temporal fluctuations in the upper range limit of trees and climatic drivers\n(growing-season temperature and precipitation).” (Line 669-678 on Pages 35-36)\nGeneral Comments\n[Comment 1] In this study, Xu Jinfeng and colleagues examined the temporal\nfluctuation of the upper range limit of trees during the late Quaternary and used a\nsuite of data and models to predict the drivers of spatial turnover of endemic flora.\nThis is a comprehensive study, aggregating various types of data, methods, and\nmodels. The work is commendable. Having reviewed several manuscripts from the\nsame lab and institute, I am consistently impressed by the rigor and scientific insight\nthis lab brings to their research. The authors of this manuscript have done a great job\nin synthesizing an interesting storyline. Beyond providing valuable facts about\ntreeline and tree limit dynamics over geological timescales, this paper substantially\ncontributes to the understanding of treeline dynamics.\n[Response] Many thanks for your encouragement, constructive comments and\nsuggestions. The point-by-point responses are listed following each\ncomment/suggestion.\nSpecific Comments\n[Comment 2] 1. Drivers of Discrepancy Between Tree Limit and Thermal Treeline:\nThe authors mention in lines 140-141: “and explained nearly 85.3% of spatial\nvariability in DTreeline across the Tibetan Plateau”. Additional details about this RF\nmodel would be valuable, including a comprehensive list of variables used and the\nnumber of instances/data points utilized.\n19\n[Response] Following your kind suggestion, we have added a table in the\nsupplementary material to give a comprehensive list of variables used in random\nforest modeling. In addition, we have included the number of data points “(N =\n116,337)” (Line 152 on Page 8; Line 497 on Page 27) used to rank the importance of\nvariables in explaining the spatial distribution of the deviation of the upper range limit\nof trees from the treeline position.\n[Comment 3] The manuscript lacks sufficient information on the explained variance.\nIt is important to note that variable importance in Random Forest models does not\nequate directly to variance explained as it does in linear regression models. While\nvariable importance measures provide insight into the relative significance of each\npredictor in making accurate predictions, these measures do not sum to 100%. They\noffer a ranking or score based on the impact of removing or altering a variable on\nmodel accuracy, but they do not represent a proportion of variance in the outcome\nvariable explained by the predictors in the model.\n[Response] Thanks a lot for pinpointing this mistake. In the previous manuscript, we\nnormalized the sum of importance values of all variables to 100%. As the reviewer\ncorrectly pinpointed, these values only offered the ranking based on how much the\naccuracy decreases when the variable is excluded. We primarily used this ranking\napproach to identify the main drivers of spatial variation in the upper range limit of\ntrees, and then used this information to develop a climate-based predictive model.\nIn order not to mislead the readers, we therefore removed the quantitative description of\n20\nthe mean decrease in accuracy. The proportion of variance explained by predictors\nwas only retained in the evaluation of constructed climate-based predictive model for\nthe upper range limit of trees.\n[Comment 4] In comparing the position of the tree limit with the thermal treeline,\nit\nwas found that 85.2% of the upper limits match the thermal treeline position.\nConsequently, the RF model was applied to the remaining 14.8% of locations where\ndiscrepancies exist between the tree limit and thermal treeline. This analysis identified\nclimatic water deficit and vapor pressure deficit as the primary drivers of\ndiscrepancies in the position of the tree limit versus the thermal treeline. Was the RF\nalgorithm also applied to the entire dataset for a more generalizable model? Given\nthat CWD and VPD are posited as the underlying drivers of discrepancies across the\nstudy area, they should predict the lack of discrepancy at some values not\ncorresponding with high-discrepancy sites. While RF models accommodate highly\nnon-linear responses, an analysis of the entire dataset could provide a more robust\necological inference beyond what has already been done.\n[Response] Thanks for pointing out this in our random forest modeling. Here, we do\nuse the entire data to construct a random forest model to understand the drivers of\nspatial variation in D . As the reviewer correctly mentioned, if excluding points\nTreeline\nwith D being equal to be zero, we would omit scenarios when drought is not a\nTreeline\nlimiting factor, potentially leading to an overestimation of past and future changes in\nD . Our current method, which used the entire data, consider varying degree of\nTreeline\n21\ndrought (D > 0) as well as no drought conditions (D = 0). This method\nTreeline Treeline\nallows for a relative robust reconstruction of historical changes in D , providing a\nTreeline\nbalanced perspective on changes in upper range limit of trees across diverse\nenvironmental conditions. To address the reviewer's concerns, we have clarified the\ndatasets used in the random forest models was the entire dataset: “The entire dataset,\nincluding cases where D ≥ 0 (N = 116,337), was used in our analysis.” (Lines\nTreeline\n496-497 Page 27).\n[Comment 5] 2. Language and clarity: The manuscript’s presentation and language\nrequire enhancements for better clarity. It frequently contains lengthy and convoluted\nsentences, alongside instances of awkward phrasing. The authors would benefit from\nthe assistance of a native English speaker with a scientific background to refine the\nlanguage. For instance, the sentence in lines 54-59 is overly long and could be\nsimplified for clarity: “There is evidence from classical paleoecological records and\nancient DNA in lake sediments showing that the upper range limit of trees, which\nwere nowadays reported in a range of 3389–4631 m based on limited field reports17-\n19, could potentially reach up to just 3500 m elevation during the Last Glacial\nMaximum (LGM) and 400–600 m higher than today in Holocene Climate Maximum\nover some locations of the Tibetan Plateau17,19,20.”\n[Response] We asked a native English speaker (Sydney.Hedberg@colostate.edu) to\nhelp us rewrite some sentences. For example, we revised the previous lengthy\nsentence to “For example, limited classical paleoecological records and ancient DNA\n22\nin lake sediments indicate that, the upper range limit of trees at some locations may\nhave reached 3,500 m elevation during the LGM23, and 400–600 m higher than today\nin the Holocene Climate Maximum24,25.” (Lines 60-64 on Page 4)\n[Comment 6] 3. Beta diversity measurement: In the final section of Results, you\nestimated beta diversity and then partitioned it to turnover and nestedness.\nInterestingly, an unusual situation arises when the null distribution of an event is\nunknown, making it impossible for us to compute the P-value or any reliable quantity\ncharacterizing some property of the system (called a statistic or metric). One glaring\nexample of this situation is beta diversity, which has been traditionally measured with\npopular indices such as Jaccard, Sørensen-Dice, and Simpson. An inherent and a\nsystematic problem with using these indices (see many papers by Ulrich, Gotelli, etc)\nhas been repeatedly reported by ecologists and biogeographers who attempted to fix\nthe problem with various techniques of normalization and standardization.\nThis longstanding issue has been addressed with the introduction of a new metric\nfounded on a solid mathematical basis and a null model developed with well-\nestablished assumptions, as detailed in Mainali et al. 2022 and Mainali and Slud\n2024. Figure 1E from the Science Advances article conclusively illustrates the\ninaccuracies of the Sorenson index in estimating beta diversity. The ecological\ninterpretations derived could significantly differ, highlighting the discrepancies\nbetween the traditional use of the Sorensen index and a statistically robust metric.\nanalysis of co-occurrence and similarity. Science Advances 8:1–9.\nMainali, K.P. and Slud, E., 2024. Null model analyses are not adequate to summarize\nstrong associations: Rebuttal to Ulrich et al. (2022). Journal of Biogeography, 51(2),\npp.263-268.\n[Response] Many thanks for insightful comments. The accuracy of the traditional\nSørensen index in estimating beta diversity is dependent on the species prevalence,\nwith the high species prevalence generally yielding misleading inferences about\nspecies associations (Mainali et al., 2022). We then quantified the extent to which\nspecies prevalence affected the Sørensen index in estimating pair-wise beta diversity\non the Tibetan Plateau. The prevalence of pair-wise species is estimated by calculating\nthe log odds ratio of the conditional occurrence probability (the occurrence rate of\npair-wise species in total counties) (Mainali et al., 2022). As show in Fig. R4a,\nendemic alpine species is expected to have a very low prevalence of 0.04 ± 0.05 over\nthe Tibetan Plateau, which is four times lower than that used as an example in\nMainali’s study (0.19 ± 0.19) (Mainali et al., 2022). The low species prevalence\nwould then have the least possibility to invalidate the interpretations of the common\nmetrics of species co-occurrence, such as the Sørensen index. This is confirmed by\nour analysis showing that there is a good correspondence (positive correlation)\nbetween the newly-developed metric of co-occurrence that is insensitive to the species\nprevalence (Mainali et al., 2022) and the most popular co-occurrence metrics such as\nthose due to Jaccard (R = 0.75, P < 0.01), Sørensen-Dice (R = 0.75, P < 0.01), and\nWhile, this new metric built upon a solid mathematical basis is generally used for\nevaluating pairwise similarities instead of assessing similarities among multiple (more\nthan two) species/sites. Although the similarities of multiple (more than two)\nspecies/sites could principally be reduced to multiple pairwise similarities, this simple\nreduction could not well capture the information shared by multiple sites (more than\ntwo) (Chao et al., 2008), potentially resulting in the loss of the degree of nestedness\n(components of beta-diversity) of species among sites (Baselga, 2010).\nIn the revised manuscript, we retained the use of the traditional Sørensen multi-\nsite dissimilarity index to calculate the beta-diversity but added the comparison\nbetween Sørensen-based method and alpha method. “The accuracy of the Sørensen\nindex to estimate beta-diversity was recently suggested to be heavily dependent on\nspecies prevalence, with high species prevalence generally yielding misleading\ninferences about species associations96. We then quantified the extent to which species\nprevalence affected the Sørensen index in estimating β-diversity on the Tibetan\nPlateau, using a new mathematical metric generally used for evaluating pairwise\nsimilarities instead of assessing similarities among multiple (> 2) species/sites. The\nprevalence of pair-wise species is estimated by calculating the log odds ratio of the\nconditional occurrence probability (the occurrence rate of pair-wise species in total\ncounties)96. The endemic alpine species on the Tibetan Plateau have a generally low\nprevalence of 0.04 ± 0.05 (Supplementary Fig. 14a), and thus, it is unlikely that\npairwise species would invalidate the interpretations of the traditional Sørensen-based\n25\nmetrics of species co-occurrence. This is confirmed by our analysis showing that there\nis a good correspondence (positive correlation) between the affinity metric of co-\noccurrence (Alpha MLE) that is not affected by species prevalence96 and the most\npopular co-occurrence metrics such as those due to Jaccard (R = 0.75, P < 0.01),\nSørensen-Dice (R = 0.75, P < 0.01), and Simpson (R = 0.83, P < 0.01)\n(Supplementary Fig. 14b to d).” (Lines 636-652 on Pages 34-35)\nReferences:\nChao, A. et al. A Two-Stage Probabilistic Approach to Multiple-Community\nSimilarity Indices. Biometrics 64, 1178-1186 (2008).\nBaselga, A. Partitioning the turnover and nestedness components of beta diversity.\nGlob. Ecol. Biogeogr. 19, 134-143 (2010).\nFig. R4 (also shown as Supplementary Fig. 14 in revised manuscript) The\ncomparison of beta-diversity calculated using different metrics over the Tibetan\ndeveloped metric and other traditional metrics of species beta-diversity (b) Jaccard,\n(c) Sørensen-Dice, (d) Simpson. The black lines were fitted with linear regressions.\nThe asterisks indicate that the correlations are statistically significant (P < 0.01 with\nthe F-test).\n[Comment 7] 4. Beta diversity pattern: Starting line 274, the authors fitted a\nGeneralized Linear Model (GLM) to explain beta diversity patterns. It was found that\npaleo fluctuations in upper range limit predicts the beta diversity pattern more\nefficiently than paleo fluctuations in growing season temperature. This finding\npresents some challenges in interpretation, as the primary mechanism through which\ngrowing season temperature could influence beta diversity patterns is by affecting\nspecies distribution. In contrast, fluctuations in the upper range limit is a direct\npredictor of beta diversity patterns. Please elaborate on the rationale and reasoning\nbehind this analysis, particularly in understanding the significance of these\npredictors.\n[Response] Many thanks for your insightful comment. The tree advance into alpine\ntreeless areas could be a direct threat for alpine endemism. Specifically, this negative\nimpact primarily occurred through the loss and fragmentation of alpine habitats, as\nwell as alterations to microclimate conditions (e.g., increased sheltering effect). These\neffects could reduce the persistence of endemic alpine species while promote upward\nexpansion of more widespread species from lower elevations (Jump et al., 2012;\nsuitable habitats for endemic alpine species. While topographically diverse\nenvironments in alpine areas may also provide fine-scale suitable microhabitats,\nallowing endemic alpine species to locally shift their physical habitat to track these\nsuitable microhabitats over short-distances (Scherrer and Körner, 2011). Therefore,\nthe presence of these fine-scale habitats could potentially buffer large-scale\natmospheric warming impacts on the survival of alpine endemism. We indeed found\nthat paleo-range limit fluctuation (R = -0.53, P < 0.01) is more important than both\npaleo growing-season temperature (R = -0.33, P < 0.01) and precipitation (R = -0.09,\nP > 0.1) fluctuation in explaining current β-diversity pattern. Therefore, we suggested\nthat shifts in upper range limits of trees could induce fragmentation and expansion of\nsuitable habitats for alpine endemism to a larger extent than temperature changes\nalone. While, we strongly agreed with the reviewer that the role of paleo-temperature\nfluctuation in reducing β-diversity should not be overlooked in projecting changes in\nβ-diversity in a warmer world.\nIn order to well resolve the reviewer’s comment, we have added the\naforementioned reasoning and rewritten implications for future projections in the\nrevised manuscript as follows: “In addition, we observed a stronger response of β-\ndiversity to paleo-range limit fluctuation than paleo-temperature fluctuation (Fig. 3a\nand b). This result is expected because shifts in the upper range limits of trees could\ninduce fragmentation and expansion of suitable habitats for alpine endemism to a\nlarger extent than temperature changes alone. Tree expansion into alpine areas could\n28\ndirectly lead to the loss of suitable habitats, and also cause habitat displacement by\naltering microclimate conditions through increased sheltering effects12,13. Although\nwarming could potentially reduce the availability of suitable habitats for endemic\nalpine species, topographically diverse environments in alpine areas may provide\nsuitable microhabitats67. These microhabitats allow endemic alpine species to track\nsuitable habitats over short-distances68. The presence of these microhabitats could\nthen potentially buffer large-scale atmospheric warming impacts on the survival of\nalpine endemism. In the coming few decades, projected rates of the upper range limit\nof trees are several orders of magnitude larger than those occurring over the last\nmillennia. In addition to drastic warming, accelerated tree expansion69 would\nexacerbate alpine habitat fragmentation and spatial isolation of remaining\nmicrohabitats12. Consequently, alpine endemics are not equipped to respond to the\ndual force from both accelerated tree expansion and drastic warming. This likely\nwould accelerate spatial homogenization of local endemic alpine species diversity,\nthreatening the regional endemic species pool.” (Lines 298-317 on Pages 15-16)\nReferences:\nJump, A. S. et al. Rapid altitudinal migration of mountain plants in Taiwan and its\nimplications for high altitude biodiversity. Ecography 35, 204-210 (2012).\nGreenwood, S. & Jump, A. S. Consequences of treeline shifts for the diversity and\nfunction of high altitude ecosystems. Arct. Antarct. Alp. Res. 46, 829-840 (2014).\nScherrer, D. & Körner, C. Topographically controlled thermal-habitat differentiation\nbuffers alpine plant diversity against climate warming. J. Biogeogr. 38, 406-416\n29\n(2011).\n[Comment 8] Line 45-46: The statement that \"Shifts in the upper range limit of trees,\nwhose fundamental niche is determined by seasonal mean temperature\" requires\nsubstantiation. If references 12-15 support this claim, it would be appropriate to cite\nthese immediately following the statement for stronger evidence.\n[Response] In the revised manuscript, we have rewritten this paragraph as “The upper\nrange limit of trees, which is defined as a connection between the uppermost group of\nupright trees of a certain height12,19, could theoretically reach their physiological\nlimits set by temperature (referred to as the treeline position)19,20 but it is often\nconfined to lower elevations due to environmental stresses and stochastic effects such\nas human disturbance12,19-22.” (Lines 54-58 on Pages 3-4)\n[Comment 9] Line 49: You mean \"spatio-temporal map\"?\n[Response] We have corrected it as “spatio-temporal map” (Line 23 and 321 on Page\n2 and 16).\n[Comment 10] Line 53: Please specify the timeframe of the Last Glacial Maximum\n(LGM).\n[Response] We have added the timeframe of the Last Glacial Maximum (LGM, ~22\nkyr BP) in the revised text. (Lines 24-25 and 60 on Page 2 and 4)\n[Comment 11] Line 72 and related paragraph: The term “upper range limit of tree\n30\ndynamics” comes across as unclear. Can you rephrase it? Additionally, your\nstatement “using the relationship between satellite-derived upper range limit of trees\nand climate in modern observations, under the common assumption that this\nrelationship has remained constant through time,” needs to be further elaborated and\njustified.\n[Response] Following your kind suggestions, we have changed the “past upper range\nlimit of tree dynamics” into “past dynamics in upper range limits of trees”. (Line 77\non Page 4) In addition, we have rephrased this sentence to make it clearer as follows:\n“Here we derived past dynamics in upper range limits of trees using paleoclimate data\nat a spatial resolution of ~1 km at 100-year intervals since LGM, using the modern\nrelationship between the satellite-derived upper range limit of trees and climate. The\nperiod since the LGM was used because its significant paleoclimatic changes led to\nnotable vegetation and floristic shifts, as evidenced by fossil records preserved in\nsediments29-32, potentially driving legacies in modern diversity patterns. In addition,\nour paleo-reconstruction of the upper range limit of trees was based on the assumption\nthat the relationship between the upper range limit of trees and climate, extracted from\nmodern observations, has remained constant through time since LGM19,20,33. This\nassumption, which is well known as the uniformitarian principle, has been widely\nused in previous pollen-based paleoclimate reconstructions33-35.” (Lines 77-88 on\nPages 4-5)\n[Comment 12] Line 119-120: The reference to Paulsen and Korner in defense of your\n31\nmethodology prompts a question about the spatial grain size of your isotherm. If your\nunit of analysis (e.g, square grid cell) is too large, then a better match between\nexpected treeline and the observed one is likely due to methodological artifact.\n[Response] The most-robust high-resolution climate data currently available is\nlimited to a spatial resolution of 1 km. Therefore, for the thermal treeline, we used\ntwo temperature thresholds (air temperature and land surface temperature) to delineate\nits position at the best available 1 km resolution. To ensure the comparability with the\nthermal treeline, we aggregated satellite-derived upper range limit of trees from the\noriginal 30 m to 1 km resolution by calculating the mean value in each 1 km grid cell.\nIn fact, we used the same approach but with less advanced satellite data (which did\nnot include vegetation structure information such as vegetation height and foliage\ndistribution) to map upper range limit of trees across the Himalayan region (Wang et\nal., 2022), and found that the trees in western and central Himalaya are absent from\nthe treeline position due to stochastic effects such as human activities. Therefore, a\nbetter match between the expected treeline and satellite-derived upper range limit of\ntrees on the Tibetan Plateau is unlikely due to methodological artifacts, but rather to\nfewer stochastic effects such as human activities.\nReference:\nWang, X. Y. et al. Enhanced habitat loss of the Himalayan endemic flora driven by\nwarming-forced upslope tree expansion. Nat. Ecol. Evol. 6, 890–899 (2022).\n[Comment 13] Line 127-128: The assertion that \"The upper range limits of trees\n32\nbelow this thermal treeline position are considered to be a result of land use\"\nrequires further explanation and citation.\n[Response] Thanks for your suggestion. Humans could remove trees from the thermal\ntreeline by grazing, logging and/or burning (Körner 2021; Wang et al., 2019).\nTherefore, to be more accurate, we have replaced land use with human activities (such\nas grazing, logging and/or burning).\nReferences:\nKörner, C. The cold range limit of trees. Trends Ecol. Evol. 36, 979-989 (2021).\nWang, Y. et al. Fire facilitates warming-induced upward shifts of alpine treelines by\naltering interspecific interactions. Trees 33, 1051-1061 (2019).\n[Comment 14] Line 133-134: The phrase “the deviation of the upper range limit of\ntrees from treeline position” seems to be missing a reference to the \"thermal treeline\nposition.\" Additionally, clarifying whether this deviation is measured as a surface\ndistance or an elevational gradient is important for interpretation.\n[Response] We have changed “the deviation of the upper range limit of trees from\ntreeline position” to “To do so, we calculated the elevational difference between the\nsatellite-derived upper range limit of trees and thermal treeline12,19,20, denoted as\nD ”. (Lines 132-134 on Page 7)\nTreeline\n[Comment 15] Line 145-150: Can you add here or in the Methods the total number\nof variables included from the four categories? What metric was used to evaluate\n33\nmodel performance, especially given the claim of over 80% spatial variation\nexplained?\n[Response] Following your kind suggestion, we modified the following sentence to:\n“used a random-forest algorithm39 to rank the importance of several variables (N =\n19) belonging to categories of environmental stresses and disturbances on D ”\nTreeline\n(Lines 134-136 on Page 7), and the Method section “To understand the drivers of\nD , we compiled several variables (N = 19) belonging to the four different\nTreeline\ncategories” (Lines 492-493 on Page 27). In addition, we used coefficient of\ndetermination between the observed D and modeled D to evaluate the\nTreeline Treeline\nmodel performance. This information has been added in the Method section.\n[Comment 16] Line 154-160: Is \"improve the explanatory power\" meant to imply\n\"maximize\"? In the rest of the paragraph, the authors presented the results of space\nfor time analysis. I have trouble understanding how a 50-year lag was evidenced from\nthe space-for-time results.\n[Response] Following your kind suggestion, we have replaced “improved” with\n“maximizes”, (Line 163 on Page 9) and added Fig. R2 into Supplementary material.\nIn addition, we have rewritten the section “We then incorporated the impact of\nmigration lag on shifts in the upper range limit of trees into the climate-driven\npredictive model. To determine the optimal lag, we tested time lags from 10 years to\n100 years in our climate-driven predictive model, by minimizing the difference\nbetween simulated and in situ shifts in the upper range limit of trees over the past\n34\ncentury across 20 sites of the Tibetan Plateau. These range limit elevation dynamics\nwere reconstructed from tree ring data using the dendrochronological method38,43-50\n(Supplementary Fig. 3 and Supplementary Table 2). Among all time lags tested,\nincorporating a 50-year time lag into the climate-driven predictive model maximizes\nits explanatory power in representing in situ shifts in the upper range limit of trees\nover the last century. Specifically, the predictive model explains nearly 70% of the\nspatial variation in range limit shifts across 20 sites on the Tibetan Plateau (R2 = 0.68\nP < 0.01) (Supplementary Fig. 4). Furthermore, after including the time lag, the trees\nacross all sites are predicted to shift upward at a rate of 20 ± 21 m decade-1, which\nbecomes closer to in situ shifts (18 ± 23 m decades-1). This confirms the robustness of\nour constructed model in the prediction of the upper range limit of trees at the century\ntimescale (Supplementary Fig. 3b).” (Lines 155-170 on Pages 8-9)\n[Comment 17] Line 170: The first definition of LGM should appear at its initial\nmention in both the abstract and the main text.\n[Response] Done as suggested.\n[Comment 18] Line 183 and related paragraph: I was delighted to read the\ndiscussion on the four-phase pattern.\n[Response] Thanks for your encouragement.\n[Comment 19] Line 206-234: The utilization of fossil pollen abundance records from\n35\nsix lakes to support your argument is commendable for its methodological diversity\nand the range of approaches and data types employed to construct a narrative on tree\nlimit dynamics in the region. However, it would be beneficial to discuss potential\nbiases in your methodology, particularly given that different tree species present in\nthe area at various times may not have pollen that is equally well-preserved.\nA discussion on how your approach mitigates the risk of underrepresentation of\ncertain species in pollen fossils due to differential preservation would strengthen the\nreliability of your findings.\n[Response] Following your constructive suggestion, we have added the following\ndiscussion into the method section: “We concentrated on the four most common high-\nelevational arboreal tree taxa - Betula, Abies, Picea, and Pinus. These tree species\nhave different production and dispersal strategies, which could greatly affect the\nrepresentation and abundance of tree taxa in fossil pollen records34,86,87. For example,\nBetula and Pinus have high pollen production rates, good wind transportability and\ngood preservation characteristics 34,88-90. Therefore, even a relatively high percentage\nat a given elevation might report a false presence. We then used a threshold of pollen\nabundance to determine if one of the tree taxa was present in the vicinity of the\nsampling site. Specifically, the pollen abundance thresholds with different\nprobabilities were determined by performing a logistic regression using pollen\nabundance (0 to 100%) as the predictor and taxon classification (0 and 1 represent\nnegative and positive classes, respectively) as the response variable using 2,434\nmodern pollen samples of the 14 key arboreal taxa from China and satellite-based\n36\nobservation and the Vegetation Atlas of China of modern trees taxa distribution86. We\nused the pollen abundance threshold at a probability of 0.9 to determine the presence\nand absence of trees. The thresholds for Betula, Abies, Picea, and Pinus are 16%,\n5%,\n9%, and 37%, respectively86. Upon identifying any one of the four tree taxa, we\nassumed that the trees were present in the vicinity of the sampling site. Not all pollen\ngrains will be equally well preserved throughout the time due to their differing\nsensitivities to environmental changes91,92, and thus, some tree taxa with very low\npollen production rates and/or poor preservation could be missed. However, the\nprobability of such underrepresentation is unlikely in this study. For example,\ncommon tree species including Corylus, Nitraria, Tamarix, Cyperaceae, Poaceae and\nFabaceae, which were widely reported to be greatly under-represented in China based\non 898 modern pollen sampling sites and 2,115 pollen data93, were not found at the\ncurrent realized range limit of trees over the Tibetan Plateau.” (Lines 574-599 on\nPages 31-32)\n[Comment 20] Section \"Reconstruction and prediction of changes in the upper range\nlimit of trees\" (Line 450): The methodology for reconstructing spatiotemporal\ndynamics using climate data from CHELSA-TraCE21k and CMIP6 future projections\nis outlined. How did you ensure the robustness of these projections, particularly given\nthe potential for significant changes in topography and climate conditions over the\nconsidered time period?\n[Response] We strongly agreed with the reviewer that the paleoclimate data from\nsimulations without constraints from paleoecological data. The CHELSA-TraCE21k\nincorporates the impact of orographic changes on paleoclimate by generating a paleo\norography at 1 km resolution for each timestep. We further analysed changes in\norography since LGM and found that a decrease in elevation is mainly located in the\nriver valleys of southeast Tibetan Plateau and over Hengduan mountain ranges\n(Figure R7). By contrast, orography changes within the 1 km buffer around the\ncurrent range limit of trees are almost negligible. Additionally, for future projections,\nwe used a multi-model ensemble to account for uncertainties in climate projections\nand applied the delta approach to correct biases in future climate projections on a\nmonthly timescale, using historical observations from the period 1970–2000. This\ninformation has been added in the revised manuscript (Lines 569–578 on Page 31).\nFurthermore, we have also included the following discussion into the revised\nmanuscript: “Furthermore, the robustness of our paleoreconstruction model relies\nupon the accuracy of the paleoclimatic forcing inputs. For example, annual\ntemperature changes during the Holocene are currently the subject of a broad debate,\nwith proxy-based reconstructions indicating a cooling trend70-72 while climate models\nsuggest warming73. Thus, reconstructing particular climate variables from\npaleoecological data is challenging over the Tibetan Plateau, due to difficulties in\ndistinguishing moisture from thermal signals35. The tracking ability of upper trees\nwith rising temperatures has given us the confidence to use paleoecological data\nlocated near their current upper range limits for reconstructing warm-season\n38\ntemperatures. After developing a robust seasonal temperature reconstruction, we\ncould then separate the moisture signal from the temperature one in paleoecological\ndata. This approach allows us to overcome the limitations of traditional\npaleoecological proxies and climate transfer functions that are notoriously\ncomplicated or distorted by the potential high collinearity between moisture and\ntemperature74.” (Lines 350-363 on Pages 17-18)\nFig. R5 (also shown as Supplementary Fig. 13 in revised manuscript) Elevation\nchanges over the Tibetan Plateau between LGM and present day. Spatial pattern\nof the elevational difference between LGM and 1950 CE. The blue lines represent\nrivers originated on the Tibetan Plateau. The bar represents elevational changes from\nLGM to 1950 CE within a 1 km buffer around the upper range limit of trees region\nbased on CHELSA TraCE-21k Digital Elevation Model (DEM) data.\nGeneral Comments\n[Comment 1] This paper presents important research regarding drivers of the upper\nrange limit of trees, as well as beta-diversity within a global biodiversity hotspot. It\nproduces important findings both for our understanding of paleo dynamics at treeline\nbut also for understanding future treeline shifts and implications for biodiversity in\nalpine mountain environments. It is most immediately important to understanding\ntreeline dynamics of the Mongolian Plateau, but findings may also be generalizable to\nalpine mountain regions globally. It produces important findings for ecologists, earth\nsystems scientists, biogeographers, and wildland managers. The conclusions and\nclaims are well-supported, and I see no flaws in analysis nor interpretation. It is\nmethodologically sound. The methods provide sufficient detail and code for the work\nto be reproduced. With minor revisions I approve its publication.\nKey results – this paper examines spatiotemporal drivers of the upper range limit of\ntrees, and models the upper range limit over the last several thousand years, for the\nTibetan Plateau. It also finds that treeline fluctuation is an important predictor of\nbeta-diversity – areas with more fluctuation have lower levels of beta-diversity, and\nareas with less fluctuation have higher levels.\nValidity – the authors employ robust methods in their examination of treeline shift and\nits effects on beta-diversity over the Holocene. Their methods include remote sensing,\nmodeling, machine-learning, in situ validation, paleo-records, and so on.\nSignificance – this paper has great significance not only to mountain regions of the\ninto the effects of treeline stability on beta-diversity in alpine environments.\nData and methodology – Their methods are multidisciplinary and rigorous.\nAnalytical approach – the analytical approach is novel, valid, and comprehensive. I\nhave only a few questions regarding their methodology (see major/minor comments\nbelow).\nSuggested improvements – I have grouped my suggestions into “major comments”\nand “minor comments” below. They pertain mainly to organization and language.\nClarity and context – overall the manuscript is clear, but I do make several\nsuggestions regarding improving the flow of the paper, fixing grammatical issues,\nmaintaining consistent tense (mainly past tense), and so on.\nReferences – yes, I believe the manuscript appropriately references previous\nliterature.\nYour expertise – I teach mountain geography, GIS, and environmental modeling\ncourses at my institution. I am a biogeographer by training. My research is not in\ntreeline dynamics, but I do use many modeling techniques described in the\nmanuscript, and I do have fair understanding of treeline and of paleo records used\nthroughout this analysis. I know remote sensing concepts. I understand most/all of the\nmodeling decisions and techniques made and their rationale.\n[Response] Many thanks for your careful review and constructive suggestions. We\ngreatly appreciate your insightful comments, which have significantly improved the\nmanuscript. We hope our revisions and explanations have satisfactorily addressed\n41\nyour questions and comments.\nSpecific Comments\n[Comment 2] In the introduction, please more formally define which kind of\n“treeline” you are reconstructing. The terms timberline, treeline, forestline, and\nscrubline are all terms used for describing different aspects of the forest to alpine\ntundra transition. Please include a definition somewhere of which you are\nreconstructing. I am guessing you are reconstructing treeline, the upper limit of\nupright tree growth. Scrubline would describe the upper limit of prostrate trees.\n[Response] Here we have added the following definition in the revised manuscript.\n“The upper range limit of trees, which is defined as a connection between the\nuppermost group of upright trees of a certain height12,19, could theoretically reach\ntheir physiological limits set by temperature (referred to as the treeline position)19,20\nbut it is often confined to lower elevations due to environmental stresses and\nstochastic effects such as human disturbance12,19-22.” (Lines 54-58 on Pages 3-4)\n[Comment 3] In the results, please be more explicit regarding which timeframe your\nresults are applicable to. For example, in line 100-115, you are referring to current\nupper range limit of trees, correct? In such locations in the paper, please add words\nlike “present” or “current” upper range limit of trees to remind the reader.\n[Response] We have added “current” in front of upper range limit of tress to remind\nthe reader. (Line 16, 93, 94, 119 etc.)\n42\n[Comment 4] For your model of upper range of trees, why not make temperature a\nvariable and directly model upper range, instead of modeling the deviation from the\ntemperature-based estimated limits? I think there could be a better explanation of the\ndecision-making here – perhaps firmer statements regarding how temperature is a\nmajor driver of treeline.\n[Response] Not all upper range limits of trees reach the treeline position determined\nby low temperature (Körner 2021; Wang et al., 2022). In the real world, much of the\nupper range limit of trees is often lower than the treeline position due to\nenvironmental stresses (e.g., critical water shortage) and disturbances (e.g., human\nactivities). Here we will firstly use growing-season temperature to delineate the\ntreeline position, and then include environmental stresses and disturbances to model\nthe observed deviation of satellite-derived range limit from the treeline position.\nThe modeled upper range limit of trees is the sum of the thermal treeline position and\nD caused by environmental stresses and stochastic effects. We separated\nTreeline\ntemperature from other climatic variables when modeling D , as there was\nTreeline\nsignificant covariation between temperature and indicators of critical water shortages,\nsuch as vapor pressure deficit and cumulative water deficit, both of which are\ninfluenced by temperature. This information has been added in the revised\nmanuscript. (Lines 468-474 on Page 26)\nReferences:\nKörner, C. The cold range limit of trees. Trends Ecol. Evol. 36, 979-989 (2021).\naltering interspecific interactions. Trees 33, 1051-1061 (2019).\n[Comment 5] Overall, I think the authors should read over the results section and\nfind ways to move some of the material into the methods section. I think the results\nshould be more results-driven. In addition, some of the results stray into discussion,\nso re-read the results and find places to move material into the discussion. Overall the\nresults section needs to be more focused on plain presentation of results.\n[Response] Done as suggested. We have moved some of the material into the\nmethods, and changed “Results” into “Results and discussions” section. For example,\nthe description of the satellite data, and the calculation of alpine land area and beta-\ndiversity were moved into the method.\n[Comment 6] If you have multiple models in this analysis, then consider naming and\nabbreviating them in the results, rather than working to describe them every time they\nare mentioned. E.g. “the current forest distribution model” or the “presence/absence\nof trees model” etc.\n[Response] Following your kind suggestion, the model we used to reconstruct upper\nrange limit of trees is named as the predictive model of range limit in the revised text.\n(Line 23, 144, 145, 148 etc.)\n[Comment 7] In the minor comments below, I often suggest adding a few clue-words\n44\nthroughout the paper to describe when results are derived from models rather than\nsome other source. Read the results over and please sprinkle these hints to remind the\nreader that models are producing many of these quantitative results (e.g. treeline\nelevation, extent of alpine tundra)\n[Response] Many thanks for this thoughtful comment. We have added a few clue-\nwords throughout the paper to distinguish model-based results from those derived\nfrom some other sources.\nMinor comments\n[Comment 8] 25 “position” – perhaps simply say “elevation”\n[Response] Done as suggested.\n[Comment 9] 29-30 “fluctuation in range limit of trees” – are you writing about\nfluctuation over time or over space? I think you should add a few words in the latter\nportion of the abstract to be clear as to whether you are talking about fluctuation\nfrom time-to-time or from location-to-location\n[Response] We have changed it into “paleo temporal fluctuations in the upper range\nlimit of trees”. (Line 669 on Page 35)\n[Comment 10] 39 should say “patterns” (plural)\n42 could cut “While, in addition to past climate change,”\n42-43 “endemic species with a restricted range” sound like species found in alpine\nhabitats. If this is mostly true, perhaps you should add “alpine” to this part.\n45\n[Response] Done as suggested.\n[Comment 11] 49-51 reword. “we required a spatial-temporally map” is confusing.\nChange to “it is essential to map the upper range limit of trees over time” or\nsomething similar.\n[Response] We have changed it into “Here we derived past dynamics in upper range\nlimits of trees using paleoclimate data at a spatial resolution of ~1 km at 100-year\nintervals since LGM, using the modern relationship between the satellite-derived\nupper range limit of trees and climate.” (Lines 77-80 on Pages 4-5)\n[Comment 12] 59 use “however,” instead of “but”\n60-61 “to create a spatial-temporally map” is awkward – say “to map\nreconstructions over time of the upper range limit…”\n62 cut “while”\n[Response] Done as suggested.\n[Comment 13] 62 I think a new paragraph could start here focusing on using refugia\nlocations vs. climate conditions as the basis for a modeling approach\n[Response] Done as suggested.\n[Comment 14] 66 “provide” not plural “provides”\n67 “tree existence” – simply say “trees”\n72 use “developed” or “estimated” instead of writing “obtained” and also write\n“limits” (plural)\n46\n[Response] Done as suggested.\n[Comment 15] 100 – “Our spatial distribution” … change to “Our analysis”\n101-104 this part is somewhat confusing to me. What does it mean that “upper trees\nare currently concentrated in the northeastern and southeastern margins of the\nTibetan Plateau”? Please reword here. Are you saying this is where treeline is highest\nwithin your study area?\n[Response] Following your kind suggestions, we have rephrased this sentence as\nfollows: “Our analysis of the satellite-derived upper range limit of trees reveals that\nalmost 94.0% (N = 11,284,982 pixels at a 30-m resolution) of pixels depicting the\nupper range limit of trees are distributed in the northeastern and southeastern margins\nof the Tibetan Plateau” (Lines 104-107 on Page 6). Here, with the percentage refers to\nthe proportion of pixels with satellite-derived upper range limit of trees.\n[Comment 16] 106 – “The remaining one” – what are you referring to? Add a few\nwords of clarification\n106 “scatteredly distributed” – awkward language – rewrite here – perhaps simply\nwrite “scattered throughout the headwater region”\n108 “tree” not “trees” I think you mean\n[Response] Following your kind suggestion, we have changed “The remaining one”\nto “The remaining 6.0% (N = 721,452) of pixels depicting the upper range limit of\ntrees are scattered throughout the headwater region of the three rivers” (Lines 109-110\non Page 6).\n47\n[Comment 17] 111 – “mean elevation” “of treeline” you mean to say?\n[Response] We have changed “The mean elevation over the Tibetan Plateau” into\n“The mean elevation of upper range limit of trees over the Tibetan Plateau”. (Lines\n117-118 on Page 6)\n[Comment 18] 116-142 – consider placing some of this material in the methods\nsection. This paragraph should be more results-oriented.\n[Response] Following your kind suggestion, we have moved some of the material\ninto the method section and kept this paragraph more focused on the results and\ndiscussions. (Lines 475-489 on Pages 26-27)\n[Comment 19] 145 “we take the knowledge from the”… rewrite as “we applied the”\n[Response] Done as suggested.\n[Comment 20] 140 versus 148 – there are different measures of deviance explained.\nBe clearer as to the differences between them. Is one calculated on training data,\nand the other is calculated on the independent validation dataset?\n[Response] We have clarified their differences in the revised manuscript.\nSpecifically, we have rewrote the paragraph as: “We used growing season temperature\nalong with climatic stressors (especially spring CWD and VPD) to parameterize the\neffect of climate change on temporal shifts in the upper range limit of trees, and\ndeveloped a climate-based predictive model using a training dataset (N = 116,337)\n(see Methods). Our climate-driven predictive model could explain nearly 84% of\n48\nspatial variation of the upper range limit of trees using an independent validation\ndataset (N = 11,634) (Supplementary Fig. 2).” (Lines 149-154 on Page 8)\n[Comment 21] 173 should you add “according to our model” to the beginning of this\nsentence? Where are these results coming from exactly – paleorecords or a model?\nAdd a few words to clarify.\n183 again, in this paragraph, remind the reader that these are model-driven results,\nnot based on fossils, pollen, etc., by adding a few words\n[Response] Done as you kindly suggested.\n[Comment 22] 218 should “grids” be instead written “grid cells”? You are\nextracting pollen-derived data at lake sites, which correspond to specific grid cells.\n[Response] Many thanks for pointing out this issue. We have modified it. (Line 231\non Page 12)\n[Comment 23] 239-243 remind the reader that these are model-driven results? E.g.\n“modeled minimum extent”, etc.\n243 cut “While”\n[Response] Done as you kindly suggested.\n[Comment 24] 238-252 some of this material sounds like introduction, some sounds\nlike methods, and some sounds like results. Do consider how you might edit/revise this\nparagraph to re-distribute and better organize the material\n[Response] Following your kind suggestion, we have rewritten the result section and\nkept this paragraph more focused on the results.\n[Comment 25] 294 I would write “unstable” rather than “fluctuated”\n[Response] Done as suggested.\n[Comment 26] 300-307 – I like the ideas here. Yes, perhaps treeline position and\ndynamics is more crucial to understanding biodiversity patterns in alpine tundras\nthan climate directly. Has anyone else made this point before? I would suggest\nreviewing the literature some more to see if your statements here support previous\nresearch.\n[Response] There is emerging evidence to show that the tree expansion could threaten\nlocal endemic alpine species diversity by reducing alpine habitat area (e.g.,\nGreenwood and Jump 2014; Dirnböck et al., 2011). To date, no research has\nquantified the impact of shifts in the upper range limit of trees on the current beta-\ndiversity of endemic alpine species. To inform the conservation of endemic alpine\nspecies diversity across spatial scales, we must consider beta-diversity which is the\ncomponent of regional diversity that accumulates from compositional differences\nbetween local endemic alpine species diversity, and understand processes that caused\nspecies composition to become more different (increased beta-diversity or biotic\n50\nheterogenization) or more similar (reduced beta-diversity or biotic homogenization)\nacross different localities. This information has been added in the introduction part.\nReferences:\nDirnböeck, T. et al. Disproportional risk for habitat loss of high-altitude endemic\nspecies under climate change. Glob. Change Biol. 17, 990-996 (2011).\nGreenwood, S. & Jump, A. S. Consequences of treeline shifts for the diversity and\nfunction of high altitude ecosystems. Arct. Antarct. Alp. Res. 46, 829-840 (2014).\n[Comment 27] 306 cut “most promising”\n313 put “such as sedimentary ancient DNA” in parentheses, I think\n[Response] Done as suggested. (Line 356 on Page 18)\n[Comment 28] 318 Fig 1 what is the difference between “Treeline” and “Upper\nrange limit of trees”? Explain earlier in the manuscript, and/or in the caption, And\nwhy are the lines different? Does it have something to do with land use at lower\nelevation / deforestation? Overall I am just confused here and I suspect others may\nalso be confused.\n[Response] The trees do not always reach their physiological limits set by\ntemperature (referred to as the treeline position), and they often remain confined to\nlower elevations due to human activities (such as grazing, logging and/or burning) or\nother disturbances and belong to the montane belt, or critical water shortage. In this\ncontext, the upper range limit of trees refers to the treeline or limits where trees are\n51\ndisplaced from the treeline position due to stochastic effects. This information has\nbeen added in the revised manuscript. (Lines 55-59 on Pages 3-4)\n[Comment 29] 337 “four distinct four-phases” …is this a typo? Rewrite as “four\ndistinct phases”\n363 add hyphen – “publicly-accessible”\n[Response] Done as suggested. (Line 419 on Page 24)\n[Comment 30] 407 “ones” – replace with a more descriptive word\n[Response] Following your kind suggestion, We changed “manually interpreted\nones” into “manually interpreted samples”. (Lines 451-452 on Page 25)\n[Comment 31] 411-412 who exactly were these interpreters? Students? Labor from\nAmazon Mechanical Turk?\n[Response] Sorry for the unclear expression. We have changed “Six well-trained\ninterpreters visually interpreted and cross-validated all samples.” into “Six well-\ntrained workers visually interpreted and cross-validated all samples using Google\nEarth and the ArcGIS platform.” (Lines 457-459 on Page 25)\n[Comment 32] 412-415 why so much validation? Could you spend some time\nexplaining your rationale?\n[Response] Thank you for your thoughtful question. Initially, we intended to map the\n52\nupper range limit of trees from Google Earth high-resolution images. However, this\nproved to be an overwhelming amount of work, so we opted to use Landsat images\ninstead. The upper range limit of trees interpreted from Google Earth images were\nthen used to validate those derived from Landsat images. This large number of\nsamples from Google Earth images demonstrates the robustness of the upper range\nlimit of trees derived from Landsat images.\n[Comment 33] 421-427 can you provide more justification for this method? Point to\nprevious literature? I’m curious as to why you can use temperature in this way, rather\nthan using temperature in a model of upper range limit of trees.\n[Response] Based on the spatial analysis, we identified main climatic drivers of the\nupper range limit of trees and used this information to develop a climate-based\npredictive model. In constructing the model, we first used a growing-season\ntemperature threshold to represent the thermal treeline position, and then included\nvariables such as environmental stresses and disturbances to account for the observed\ndeviation of the satellite-derived upper range limit from the treeline position (D ).\nTreeline\nThe modeled upper range limit of trees is the sum of the thermal treeline position and\nD caused by environmental stresses and stochastic effects. We separated\nTreeline\ntemperature from other climatic variables when modeling D , as there was\nTreeline\nsignificant covariation between temperature and indicators of critical water shortages,\nsuch as vapor pressure deficit and cumulative water deficit, both of which are\ninfluenced by temperature. This information has been added in the revised\nmanuscript. (Lines 468-489 on Pages 26-27; Lines 504-524 on Pages 28-29)\n53\n[Comment 34] 442-445 here and elsewhere, explain the sources of the data. Or, if\nyou don’t explain the sources here, mention that tables in the supplementary material\nexplain where the datasets are from, their original spatial resolution, etc. Just ensure\nthat somewhere in the paper you direct readers to sources of the variables and the\ncharacteristics of the data.\n[Response] Following your kind suggestion, we have added a table in the\nsupplementary material to give a comprehensive list of variables used in random\nforest model (Supplementary Table 1). This list of variables includes variable names,\ndata sources, original spatial resolution, and time period.\n[Comment 35] 484-485 “the life-form tree was” – just rewrite as “trees were”\n495 add “m” after “4500” for meters\n496 rewrite as “We concentrated on”\n499 rewrite as “were determined” (past tense)\n[Response] Done as suggested.\n[Comment 36] 505-506 the amount for Pinus seems high, and Betula seems\nsomewhat high also, as if those trees should still be “present” at these pollen levels.\nIs there are reason why pollen levels could be this high yet the tree species could be\n“absent”? It’d be good to cite some literature explaining perhaps that pollen of some\ntaxa is transported longer distances, and/or that these taxa produce more pollen than\nothers – if that’s the case.\n54\n[Response] Following your constructive suggestion, we have added the following\ndiscussion into the method: “We concentrated on the four most common high-\nelevational arboreal tree taxa - Betula, Abies, Picea, and Pinus. These tree species\nhave different production and dispersal strategies, which could greatly affect the\nrepresentation and abundance of tree taxa in fossil pollen records34,86,87. For example,\nBetula and Pinus have high pollen production rates, good wind transportability and\ngood preservation characteristics 34,88-90. Therefore, even a relatively high percentage\nat a given elevation might report a false presence. We then used a threshold of pollen\nabundance to determine if one of the tree taxa was present in the vicinity of the\nsampling site. Specifically, the pollen abundance thresholds with different\nprobabilities were determined by performing a logistic regression using pollen\nabundance (0 to 100%) as the predictor and taxon classification (0 and 1 represent\nnegative and positive classes, respectively) as the response variable using 2,434\nmodern pollen samples of the 14 key arboreal taxa from China and satellite-based\nobservation and the Vegetation Atlas of China of modern trees taxa distribution86. We\nused the pollen abundance threshold at a probability of 0.9 to determine the presence\nand absence of trees. The thresholds for Betula, Abies, Picea, and Pinus are 16%,\n5%,\n9%, and 37%, respectively86. Upon identifying any one of the four tree taxa, we\nassumed that the trees were present in the vicinity of the sampling site. Not all pollen\ngrains will be equally well preserved throughout the time due to their differing\nsensitivities to environmental changes91,92, and thus, some tree taxa with very low\npollen production rates and/or poor preservation could be missed. However, the\n55\nprobability of such underrepresentation is unlikely in this study. For example,\ncommon tree species including Corylus, Nitraria, Tamarix, Cyperaceae, Poaceae and\nFabaceae, which were widely reported to be greatly under-represented in China based\non 898 modern pollen sampling sites and 2,115 pollen data93, were not found at the\ncurrent realized range limit of trees over the Tibetan Plateau.” (Lines 574-599 on\nPages 31-32)\n[Comment 37] 507 again simplify “the life-form tree was” to “trees were”\n511 “We aimed to estimate” (past tense)\n516 “was” not “is”\n[Response] Done as suggested.\n[Comment 38] 511-518 and elsewhere – make sure you are writing mainly in past\ntense\n[Response] Thank you very much for the extensive and careful editing. We have\ncarefully revised the manuscript to ensure that we use past tense consistently.\n[Comment 39] 534 Cut “While” or replace with “on the other hand”\n[Response] Done as suggested. (Line 645 on Page 35)\n[Comment 40] 551-558 please be specific about where the alpine land area data are\ncoming from. Is this initially based on your model of forest cover (from earlier in the\n56\nmethods section)?\n[Response] Following your constructive suggestion, we have rewritten previous\nsentences as “In addition, to assess whether warming-induced upslope tree expansion\nwould potentially lead to habitat loss, we also calculated the alpine land area by\nsumming the surface area above the mean elevation of the reconstructed upper range\nlimit of trees in each county, using a 30-meter spatial resolution digital elevation\nmodel (DEM)100. For each mountain defined in the GMBA mountain inventory v2.0\ndataset101, we created a 100 km radius buffer around a target upper range limit pixel,\nand then summed up the surface area of all 30 m resolution pixels with elevations\nhigher than that of the target pixel.” (Lines 678-685 on Page 36)\n[Comment 41] Is it at all worth mentioning somewhere in the paper that treeline\nposition changes albedo? I imagine Arctic treeline latitude is more crucial to altering\nEarth’s albedo, but perhaps you could somewhere expand upon the significance of\nyour work by stating that shifting mountain treelines could affect albedo (at least\nlocally) – it’s up to you.\n[Response] Many thanks for this insightful suggestion. We agreed with the reviewer\nthat a shift in the upper range limit of trees generally decrease surface albedo, which\nwould enhance local warming and accelerate snow and glacial melt in alpine region.\nOn the other hand, tree expansion into alpine regions could enhance\nevapotranspiration that might cool the surface. It remains an open question whether or\nnot the tree expansion into alpine treeless region would amplify the warming.\ndynamics on endemic alpines species diversity pattern, we will not include this\nanalysis in the main text. Instead, we will plan to develop a separate paper on this\nexcellent idea using a coupled land-atmosphere model approach.\n[Comment 42] Supplementary Fig. 4 – edit this caption because there are some\nawkward sentences and misspellings.\n[Response] Following your kind suggestions, we have rephrased this caption to make\nit clearer as follows: “Supplementary Fig. 7 Validation of model-based\nreconstruction of upper range limit of trees using fossil pollen abundance records\nfrom lake sediments during the Holocene period. a, Spatial distribution of the\ncurrent upper range limit of trees at 1-km resolution, overlaid with six lakes\ncontaining Holocene fossil pollen records. b, The comparison between pollen- and\nmodel-based reconstruction of upper range limit of trees at each lake site during the\nHolocene. The gray shading represents pollen-derived arboreal pollen percentage,\nand\nthe red solid line represents paleo changes in growing-season temperature relative to\n1970–2000 mean. The dotted line with ticks represents the time span when trees\nwere present at the lake sites, while the black and blue dotted lines indicate pollen-\nbased and model-based reconstructions, respectively (see Methods).”\n[Comment 43] Supplementary Table 1 – please add the original spatial resolution of\neach variable. I think it is important to be clear that some variables had finer or\n58\ncoarser resolution.\n[Response] Done as suggested.\n[Comment 44] Supplementary Table 5 – please explain the source(s) of alpine flora\nand their elevational ranges\n[Response] Following your kind suggestions, we have rephrased this caption as\nfollows: “List of alpine endemic flora and their elevational range, derived from the\nFlora of China94 and the Flora of Tibet.”\n[Comment 1] I continue to be impressed by the authors' effort in developing this time\nseries of treeline dynamics on the Tibetan Plateau throughout the late Pleistocene.\nThe field effort, integration of current and paleo data, and sophisticated modelling\nare commendable. While I have only a few technical clarifications, my primary\nconcern remains with the ecological interpretation of the study, specifically the link\nbetween treeline fluctuations and endemic species turnover in the region. To begin, it\nremains unclear why treeline changes would impact endemic species turnover, a\nmetric that reflects compositional differences between sites. I am also uncertain why a\nvariable that reflects local processes (treeline variability) should be used to interpret\nchanges in a regional-scale biodiversity variable (e.g., turnover across counties). I\nunderstand the authors' rationale for this analysis—it serves as a case study for the\napplication of the impressive time series. However, I believe that the research question\nthey aim to address may be too ambitious given the spatial limitations. My suggestion\nis to focus more on the technical accomplishments of the treeline modelling,\nhighlighting its broader applicability beyond the Tibetan Plateau and its potential as\nan indicator of environmental instability in studies examining the impact of\nenvironmental change on biodiversity patterns.\n[Response] Thanks a lot for your constructive comments that have greatly improved\nthe clarity of manuscript. Instead of only focusing on the reconstruction of paleo-\nupper range limit of trees, we dared to retain the section on the role of paleo shifts in\nupper range limit of trees in present-day alpine endemic diversity. This is because\nshifts in the upper range limits of trees are directly connected to the expansion or\ncontraction of habitats for endemic alpine species that inhabit areas above these\nlimits. Unlike previous studies that have explored the influence of paleoclimatic\nlegacies on present-day biodiversity, no research to date has addressed the potential\nrole of paleo-range limit in alpine regions. Thanks to your constructive comments,\nwe\nrealized that our previous expression will definitely mislead readers that there is a\nspatial mismatch between shifts in upper range limit of trees and beta-diversity. In\nfact, although paleochanges in upper range limit of trees were limited to a relative\nnarrow range, these shifts represents significant changes in available habitat for alpine\nendemic species that live above these limits within a region. We aimed to understand\nthe potential contribution of spatial variation in fluctuations of available habitat area,\nwhich is connected to paleochanges in upper range limit of trees, to the spatial pattern\nof beta-diversity across regions. Here we used Sørensen-based multiple-county\ndissimilarity to calculate the beta-diversity within a region that is consisted of a target\ncounty and its surrounding counties. Importantly, the paleo-range limit fluctuation,\nalong with paleoclimatic fluctuation, were averaged within the region. Therefore, the\nbeta-diversity and its drivers are spatially matched. Last but not least, understanding\nhow paleo changes in the upper range limits of trees influenced the spatial patterns of\nbeta diversity could provide valuable insights for more accurately predicting future\nshifts in alpine endemism diversity in a rapidly warming world. All of this\ninformation has been clarified in the revised manuscript, and we hope that our\nresponses attached below could well resolve the reviewer’s concerns.\nMinor Comments\n[Comment 2] L32-33: This statement is somewhat argumentative. It would benefit\nfrom a distinction between natural climate variability and anthropogenic climate\nchange.\n[Response] We have rephrased this sentence as “Thus, our results suggest that\nanthropogenic-caused climate change on decadal-to-centennial timescales could lead\nto higher fluctuations in range limits than orbitally forced climate variability on\ncentennial-to-millennium timescales, which consequently could cause spatial\nhomogenization of endemic alpine species composition, threatening Tibetan endemic\nspecies pool.” (Lines 32-36 on Page 2)\n[Comment 3] L50-51: The statement that treeline changes affect community\ndissimilarity across sites lacks a clear explanation of the underlying mechanism.\n[Response] Following your constructive suggestion, we have added the mechanisms\nunderlying impacts of changes in paleoclimatic and paleo-upper range limit of trees\non the spatial pattern of beta-diversity into the revised manuscript: “The pattern of\ncontemporary beta-diversity, which quantified the variation in species composition\nwithin a region, is generally interpreted from past climate regimes and their\noscillations19,20, but the role of paleo changes in the upper range limit of trees remains\nunknown. The beta-diversity was generally contributed by two antithetic processes\n(species replacement and changes in species richness)21. Compared to regions that\nexperienced unstable climate and low fluctuation of tree range expansion/contraction\nover glacial-interglacial cycles, regions with a relatively stable environmental\ncondition would enhance niche differentiation. This would potentially allow more\ndiverse endemic alpine species to experience lower rates of species extinction and\nhigher rates of speciation, increasing beta-diversity through species replacement. By\ncontrast, regions with unstable past climates and high fluctuation of tree range\nexpansion/contraction would reduce suitable habitats and increase the fragmentation\nof these habitats for the survival of endemic alpine species, leading to higher\nextinctions and increasing beta-diversity through nestedness. Therefore,\nunderstanding how paleo changes in the upper range limits of trees influenced the\nspatial pattern of beta diversity could offer valuable insights for more accurately\npredicting future shifts in alpine endemism diversity under a rapidly warming world. ”\n(Lines 49-66 on Pages 3-4)\n[Comment 4] L64-66: This phrasing underplays the potential contributions of\npaleoecology. You leverage paleoecological data to validate treeline changes, and it\nis indeed challenging to create spatially continuous maps of treeline change using\npaleoecological information alone.\n[Response] We have rephrased it as “Our previous understandings of paleo changes\nin the upper range limit of trees is mostly derived from analysis of fossil plant data\nincluding pollen, macrofossils, and charcoal. For example, classical paleoecological\nrecords and ancient DNA in lake sediments indicate that, the upper range limit of\ntrees at some locations may have reached 3,500 m elevation during the LGM26, and\n400–600 m higher than today in the Holocene Climate Maximum27,28. While, relying\nsolely on fossil data at the site level to reconstruct the spatiotemporal evolution of the\nupper range limits of trees remains highly uncertain, especially on the Tibetan plateau\nwhere a handful of sites were available for a few time periods since the Last Glacial\nMaximum (LGM, ~22 kyr BP).” (Lines 71-80 on Pages 4-5)\n[Comment 5] L120-124: I have some concerns about potential underestimation of\npast treeline limits, particularly during periods with minimal human impact.\nThis is\nsomewhat contradicted by the high (84%) alignment of sites with growing season\ntemperature limits.\n[Response] We fully agreed with the reviewer’s concern. Before developing the\npredictive model for paleoreconstruction, we ranked the relative importance of all\npotential variables in modelling the spatial variation of D (the observed\nTreeline\ndeviation of the satellite-derived upper range limit from the treeline position) using all\nthe data at 1 km resolution. Our results showed that human intervention contributed\nlittle to explaining the spatial variation in D . The impact of human intervention\nTreeline\non D could only be detected when the intensity of human intervention\nTreeline\nrepresented by the human footprint index exceeded 5 (Supplementary Fig. 1n), while\nthe mean level of present-day human intervention is only 4.08 ± 2.65 (unitless)\naround the upper range limit of trees over the Tibetan Plateau.\nAlthough the role of human intervention in explaining the spatial variation of\nD is rather small, there are still some areas that have the relatively high level of\nTreeline\nhuman intervention, with the human footprint index exceeding 5. To overcome the\npotential underestimation of paleo-upper range limit of trees in these areas but with\nmuch less human intervention in the latest glacial-interglacial period, we therefore\nonly used the data from intact areas without significant human impacts to develop the\npredictive model for paleoreconstruction. This information has been clarified in the\nrevised manuscript.\n[Comment 6] L135-136 and Methods: Although a wide array of predictors is tested,\nit would be preferable to present hypotheses for each variable's potential impact.\nAdditionally, clarify if these environmental variables measure deviations between\nexpected and observed treeline locations.\n[Response] Following your constructive suggestion, we have added the following\nsentences into the Methods section: “The mechanisms for the potential impact of\nthese variables on tree growth and establishment at thermal treeline positions were\nformulated as follows. For example, the water availability (represented by CWD and\nVPD) could constrain tree growth and become decisive for seedling establishment,\ntherefore affecting the tree establishment at the thermal treeline position. The CWD,\nwhich is calculated as potential evapotranspiration minus precipitation, indicates the\nquantity of water needed to meet evaporative demand and measures the drought stress\non plants, and VPD, as a function of relative humidity and temperature, represents an\nindependent moisture stress on plant photosynthesis. Furthermore, summer\ndehydration, which is calculated as a product of wind speed, temperature and relative\nhumidity82, generally occurs during daylight hours under conditions characterized by\nhigh temperatures, low relative humidity and strong winds, thereby leading to trees\ntranspiring moisture at a faster rate than their roots take it up from the soil. Such\nsummer dehydration could further aggravate the moisture stress (CWD and VPD) on\ntree growth and tree establishment at thermal treeline positions. The relationship\nbetween cloud cover and tree growth at thermal treeline position is complex. The high\ncloud cover could reduce the solar radiation, and potentially undermine the plant\npotentiation and tree growth at thermal treeline positions. While the cloudy days\nare generally associated with the relatively high-water availability and reduced\ntranspiration during the growing season; and this could reduce the radiative cooling\nand trap heat during the night, potentially allowing trees to survive extreme cold\ntemperatures at the relatively high elevation25. The minimum temperature of the\ncoldest month is used as the freezing stressor that could directly control the absence of\ntrees from the thermal treeline position22. In addition to these macroclimatic drivers,\nwe also include local topographical factors (e.g., surface curvature) as a surrogate of\nmicroclimate conditions, because of the potential deviation of microclimate from\nmacroclimate at the high elevation.” (Lines 560-586 on Pages 30-31)\nFurthermore, we have added the leading environmental variables measure\ndeviations between expected and observed treeline locations: “Furthermore, we\nassessed how spatial changes of the leading factors (CWD and VPD) modulated the\nspatial distribution of D . For example, the abrupt D transition occurring at\nTreeline Treeline\nthe CWD threshold of ~100 mm yr−1 (Supplementary Fig. 1b) separated the southern\nTibetan Plateau, with low CWD and D values, from northern Tibetan Plateau\nTreeline\nwith high ones at ~35 °N (Fig. 1a). The transition of D along the VPD gradient\nTreeline\nfurther sharpens the observed south–north contrast, with D exhibiting a\nTreeline\npronounced increase when the VPD exceeded 0.4 hPa (Supplementary Fig. 1c).”\n(Lines 157-164 on Pages 8-9)\n[Comment 7] L175-187: This section addresses projected end-of-century treeline\nchanges rather than “Paleo reconstruction of upper range limit of trees since LGM.”\n[Response] We have changed it into “Prediction of upper range limit of trees changes\nfrom the LGM to the end of this century” (Lines 194-195 on Page 10)\n[Comment 8] L216-229: Please expand the explanation of how palynological data\nand treeline estimates were compared for congruence.\n[Response] Thanks for pointing out this omission. We have lacked a detailed\ndescription of how the pollen-inferred changes in paleo upper range limit of trees\nwere derived using six lakes situated below, at and above the present upper range\nlimit of trees. Such missing description led to an unclear communication regarding\nhow palynological data and upper range limit estimates were compared for\ncongruence. In order to well resolve the reviewer’s concern, we have therefore added\nthe following sentences into the revised manuscript. “Specifically, for a given paleo\ntime period, we obtained the elevational range of the upper range limit of trees by\nidentifying the presence and absence of trees from fossil pollen records in lake\nsediments at different elevations. For example, during the early to middle Holocene\n(~10–4 kyr BP), trees were detected in Wuxu Lake (3,710 m), Muge Co Lake (3,780\nm) and Yidun Lake (4,470 m), but was absent in Koucha Lake (4,540 m) and Saiyong\nCo Lake (4,545 m) (Supplementary Fig. 7b). We then obtained that the elevation of\npaleo-upper range limit was situated between 4,470 m (Yidun Lake) and 4,540 m\n(Koucha Lake) during ~10–4 kyr BP. By further comparing elevational changes\nbetween different paleo time periods, we could obtain paleochanges in upper range\nlimit of trees using fossil pollen records.” (Lines 245-255 on Pages 12-13)\n[Comment 9] L247-250: There seems to be a spatial mismatch between the\ncommunity dissimilarity metric at the county level and the local-scale treeline\nchanges. County-level compositional changes may be influenced by factors beyond\ntreeline dynamics.\n[Response] Sorry for our confused expression. Within each target county, we did not\nhave species data down to the sub-county or plot scale. We then calculated the\namount of differentiation between the target county and its surrounding counties to\nrepresent β-diversity for this region. Here paleochanges in upper range limit of trees\nrepresents changes in available habitat for alpine endemic species that live above\nthese limits within a region, and we aim to understand the potential contribution of\npaleo fluctuation in avaible habitat area to the spatial pattern of β-diversity for those\nalpine endemic species (across regions). To understand mechanisms underlying the\nspatial pattern of β-diversity, we obtained the regional-mean value of potential\ndrivers, e.g., paleo-range limit fluctuation and paleoclimatic fluctuation, by averaging\nthem across the target county and its surrounding ones. Therefore, the beta-diversity\nand its drivers are matched spatially. In order to well resolve the reviewer’s concern,\nwe have clarified that Sørensen-based multiple-county dissimilarity was used to\ncalculate the beta-diversity (β-diversity) of alpine endemic flora within a region that is\nconsisted of a target county and its surrounding ones. In addition, corresponding\nmodifications were also performed throughout the manuscript.\nWe totally agreed with the reviewer that, in addition to the paleo-range limit\nfluctuation, there are other factors that are responsible for the spatial pattern of β-\ndiversity. We have indeed explored the potential role of other drivers, including mean\nclimate seasonality, in contributing to the spatial pattern of β-diversity, while found\nthat β-diversity did not vary significantly across spatial gradients in mean climate and\nseasonality.\n[Comment 10] L264-268: This phrasing could be simplified to indicate that you\nestimated a standardized effect size of beta-diversity. However, given the\ndecomposition of beta-diversity into turnover and nestedness components, it’s unclear\nwhy this additional analysis is necessary.\n[Response] Thanks for your constructive comment and we fully agreed with the\nreviewer’s viewpoint. Calculating the deviation of observed β-diversity from a null-\nmodel generated stochastic expectation was mainly used to remove the component of\nβ-diversity resulting from differences in local alpha diversity and sizes of regional\nspecies. By decomposing β-diversity into the two components (nestedness and species\nreplacement), the nestedness component indeed reflects differences in local species\nrichness. In the revised manuscript, we have removed the analysis of β-deviation but\nretained the decomposition of β-diversity into nestedness and species replacement\ncomponents. Such decomposition should not only help us understand mechanisms\nunderlying the spatial organization of biodiversity, but also bring deeper insights into\nfuture projections of biodiversity changes if the drivers of spatial patterns of β-\ndiversity were assumed to drive their temporal changes.\n[Comment 11] L325-329: I understand the argument for local data reflecting\ncompositional changes over time, but what is the mechanism for observing this\npattern at larger spatial scales?\n[Response] Thank you for your constructive comment. Our results highlighted that\nthe spatial pattern of β-diversity and its two components (species replacement and\nnestedness) could be best described by both paleo-range limit and paleoclimatic\nfluctuation. Furthermore, these paleo fluctuations enhanced the nestedness component\nbut decreased the species replacement. These results provided mechanistic insights\ninto the spatial organization of biodiversity. If the drivers of spatial patterns of β-\ndiversity and its two components also drive their temporal changes, more fluctuated\nrange-limit and temperature in a warmer future would imply an increase in the\nnestedness component but a decline in the species replacement component of β-\ndiversity. Although such space-for-time assumption might be useful in predicting\nclimate-change impacts on biodiversity (Blois et al., 2013), it might only be valid\nwhen the spatial variation in environmental factors (e.g., climate and upper range\nlimit) could encompass its temporal gradient. Therefore, to be more accurate, we have\nrephrased the sentence as: “If the drivers of spatial patterns of β-diversity and its two\ncomponents also drive their temporal changes, the enhanced paleo-range limit\nfluctuation, coupled with paleo-temperature fluctuation, could shift the pattern of\ndissimilarity derived from spatial replacement to a pattern of dissimilarity derived\nfrom nestedness, which has important implications for conservation planning.” (Lines\n353-357 on Page 18).\nAdditionally, in the introduction section, we have highlighted that the spatial\nanalysis could potentially offer valuable insights into future predictions regarding the\nimpacts of climate change on alpine endemic species diversity. “Therefore,\nunderstanding how paleo changes in the upper range limits of trees influenced the\nspatial patterns of beta diversity could offer valuable insights for more accurately\npredicting future shifts in alpine endemism diversity under a rapidly warming world.”\n(Lines 63-66 on Page 4)\nReference\nBlois. J et al. Climate Change and the Past, Present, and Future of Biotic Interactions.\nScience 341, 499-504 (2013).\n[Comment 12] Additional Details and points related to the methods. Figure 1a:\nAdding the longitudinal axis for observed and predicted treeline changes would\nenhance clarity.\n[Response] Done as suggested. We added a longitudinal axis for the thermal treeline\n(calculated based on the mean growing-season temperature, see Methods) and the\nsatellite-derived upper range limit of trees.\nFig. 1 Elevation distribution of the multi-satellite-derived upper range limit of\ntrees across the Tibetan Plateau. (a) Spatial distribution of the upper range limit of\ntrees at 1 km2 resolution (left), and longitudinal/latitudinal distribution of treeline\nelevation and upper range limit of trees elevation (top/right). The solid lines represent\nthe median values for each elevational category, calculated at 1 km intervals\nlongitudinally. The curves have been smoothed using a Gaussian function, with the\nshaded regions indicating one standard deviation. (b) The relationships between the\nmulti-satellite-derived upper range limit of trees and visual interpretation from Google\nEarth images, with the 1:1 line (dashed) and the regression line (continuous). R2 and\nME indicate the coefficient of determination and the mean error, respectively.\n[Comment 13] L486-488: Further explanation of the thresholds defining climate\nbands as the growing season would strengthen this argument. Clarify why these limits\nwere chosen over those provided by indices like those from\n[CLIMDEX](https://www.climdex.org/learn/indices/), and indicate the dataset used.\n[Response] Following your kind suggestion, we have added the dataset, such as\nWorldClim and the Moderate Resolution Imaging Spectroradiometer (MODIS), that\nwere respectively used to calculate the mean growing-season air temperature at 2 m\nabove the ground and land surface temperature in the revised manuscript. (Lines 509-\n512 on Pages 27-28)\nOur definition of the temperature threshold for start and end of growing season is\nderived from continuous ground-temperature measurements from more than 30\nthermal treeline sites between 70 °N and 56 °S at the global scale. In contrast, in\nCLIMDEX data, air temperature threshold used to define the start and end of growing\nseason is 5°C. Since our focus is on the upper range limits of trees, we defined the\ngrowing season using the temperature threshold derived from a meta-analysis of\nthermal treeline sites, instead of relying on the CLIMDEX dataset.\nIn order to well resolve the reviewer’s concern, we have added the following\nsentence into the Method section: “In addition, a meta-analysis of thermal treeline\nsites at the global scale also showed that the start and end of growing season was\ndefined as date when the ground temperature at 10 cm depth exceeds 3.2°C24. Such\nground-temperature threshold corresponds to MODIS-derived land surface\ntemperature threshold of 0.7 °C. Therefore, the start and end of the growing season\nwere determined as the date when daily T rises above 0.9 °C, the date at which it\nair\nfalls below 0.9 °C, and the daily T passes 0.7 °C, respectively.” (Lines 518-524\nsurface\non Page 28).\n[Comment 14] L514-523: The modelling process for D-Tree is unclear. After using a\nrandom forest model to identify the relevant predictors, how were these predictors\napplied? Additionally, how did you handle resolution discrepancies between the\ntreeline model (30m) and environmental variables (typically at ~250m to 1km)?\n[Response] In the revised manuscript, we have added how these predictors were\ncalculated. For example: “The CWD, which is calculated as potential\nevapotranspiration minus precipitation, indicates the quantity of water needed to meet\nevaporative demand and measures the drought stress on plants, and VPD, as a\nfunction of relative humidity and temperature, represents an independent moisture\nstress on plant photosynthesis. Furthermore, summer dehydration, which is calculated\nas a product of wind speed, temperature and relative humidity82, generally occurs\nduring daylight hours under conditions characterized by high temperatures, low\nrelative humidity and strong winds, thereby leading to trees transpiring moisture at a\nfaster rate than their roots take it up from the soil. Such summer dehydration could\nfurther aggravate the moisture stress (CWD and VPD) on tree growth and tree\nestablishment at thermal treeline positions.” (Lines 564-574 on Page 30)\nWe further added the following sentence “We then applied our parsimonious\nclimate-D model to simulate changes in D at a spatial resolution of ~1 km\nTreeline Treeline\nat 100-year intervals since LGM to the end of this century (mean over the period\n2080–2099). All of these predicators from LGM to the end of this century were\nderived from CHELSA-TraCE21k, WorldClim and CMIP6 model projections,\nrespectively (see “Assembling climate variables”). ” (Lines 589-593 on Page 31)\nIn addition, in terms of issues related to resolution discrepancies, we have also\nadded the following sentence into the revised text: “Before calculating D , we\nTreeline\naggregated satellite-derived upper range limit of trees elevation from the original 30\nm to 1 km resolution by calculating the mean value in each 1 km grid cell.” (Lines\n528-530 on Page 28)\n[Comment 15] L533-540: Why were two different climate sources used? While the\nuse of CHALSA downscaling of the TRaCE experiment is understandable, why not use\nCHALSA climatology for current conditions instead of WORLDCLIM? Offsets\nbetween these sources could influence treeline predictions. Also, specify which CMIP6\nmodels were used, and confirm if they were downscaled by WORLDCLIM or\nCHALSA.\n[Response] Sorry for our confused expression. We used WorldClim data to map the\npresent-day upper range limit of trees, and estimated the paleo-range limit by\ncombining the current climatology from WorldClim with paleo-climatic changes\nderived from CHELSA. Furthermore, we generated future climate data at a spatial\nresolution of ~1 km by downscaling and bias-correcting CMIP6 model simulations\nusing WorldClim data. Therefore, discrepancies among different data sources in\nestimating dynamics in upper range limit of trees were resolved through bias\ncorrection. In addition, we have added a table in the Supplementary material to give a\ncomprehensive list of CMIP6 models (Supplementary Table 3). This information has\nbeen added in the revised manuscript.\nWe used WorldClim instead of CHELSA to map present-day upper range limit\nof trees, as WorldClim provides an observation-based dataset by integrating weather\nstation data with remote sensing observations. This makes it more suitable than\nCHELSA using topography-based downscaling of climate reanalysis data.\nSpecifically, the WorldClim aggregated global weather stations across 1970–2000,\nand used thin-plate spline interpolation method to generate climate data at a spatial\nresolution of ~1 km. The covariates used for interpolation includes elevation, distance\nto the coast and MODIS derived maximum and minimum land surface temperature as\nwell as cloud cover (Fick and Hijmans 2017). By contrast, the CHELSA at a spatial\nresolution of ~1 km is derived from the ERA-Interim (ECMWF Re-Analysis-Interim)\nreanalysis data using the topography-based downscaling algorithm (Karger et al.,\n2017). This topography-based downscaling algorithm incorporated the potential\nimpact of orographic predictors including wind effect (windward and leeward),\nvalley exposition (isolated valleys in high mountain), and lapse rates on air temperature at a\nhigh spatial resolution.\nReferences\nFick, S. E. & Hijmans, R. J. WorldClim 2: new 1-km spatial resolution climate\nsurfaces for global land areas. Int. J. Climatol. 37, 4302-4315 (2017).\nKarger, D. N. et al. Climatologies at high resolution for the earth's land surface areas.\nSci. Data 4, 170122 (2017).\n[Comment 16] L580-587: Why did you opt for this approach instead of correcting for\npollen production, as is standard in palynological analyses?\n[Response] We fully agreed with the reviewer’s comment. The pollen productivity\nwould be more accurate to assess the occurrence of a certain pollen taxon. While,\nthere is no available pollen production data for the main pollen taxa from the eastern\nof Tibetan Plateau, and such absence was mainly caused by the fragmented landform\nand landscape. In addition, the pollen production for a certain taxon cannot be simply\nused to estimate that from other regions. We then used a threshold of pollen\nabundance to determine if one of the tree taxa was present in the vicinity of the\nsampling site. This information has been clarified in the revised manuscript.\nTo Reviewer #2\nGreat job thoroughly revising the manuscript. I am satisfied with the revision and\nhave no additional comments.\n[Response] We sincerely thank the reviewer for their efforts on our manuscript and\ngreatly appreciate their contributions to its improvement.\nTo Reviewer #3\nI am satisfied with the revisions made by the authors, specifically those made in\nresponse to my major and minor comments.\n[Response] Thank you very much for your valuable and encouraging comments and\nsuggestions. Your suggestions have been incredibly helpful in improving the quality\nof our manuscript.\nTo Reviewer #1\nGeneral Comments\n[Comment 1] I appreciate the author’s patience regarding my ongoing need for\nfurther clarification—particularly concerning their beta-diversity analysis.\nMethodologically, I am now satisfied with the additions and clarifications. It is clear\nwhy certain variables were employed, how palynological data and treeline estimates\nwere compared, and the thresholds defining climate bands.\n[Response] We appreciate your efforts on constructive comments that have greatly\nimproved manuscript.\nSpecific Comments\n[Comment 2] Moreover, now that I understand their rationale for these analyses,\nI\nconcur that we must broaden our examination of the drivers of biodiversity change\nfrom the \"ultimate causes\" (e.g., climate change) to more immediate causes (range\ndynamics and environmental suitability). However, the text they propose to explain\nthis point (Lines 50-67) requires editing for clarity in the arguments", "reviewer_specifically": "nL53-54: The definition of beta diversity is output place as it’s very hard to link it\nto the points made below on the effects of changes in tree ranges.\nIn this section, L55-62 is unclear if you mean regions with relatively stable\nenvironmental conditions and low fluctuations in tree range changes. I recommend\ncondensing this section and focusing on stating that under unstable climate\nconditions, low fluctuation in tree range expansion or contraction would result in\nlower losses of endemic alpine species than in areas where tree ranges undergo\nsignificant fluctuation. This clarifies why we need to understand how paleo changes\nin the upper range limits of trees influenced diversity patterns.\n[Response] Following your constructive suggestion, we have condensed the\nparagraph into “Specifically, in areas with stable climate conditions, low fluctuation\nin the position of the upper range limit of trees would result in reduced habitat loss\nand fragmentation, thereby creating stable environmental conditions that are more\nconducive to the survival of endemic alpine species, compared to areas where the\nposition of upper range limit of trees fluctuates more significantly. We should\ntherefore expand our analysis of the drivers behind the current spatial patterns of beta-\ndiversity, shifting focus from climate change to include a more immediate factor –\nrange limit dynamics. This broader understanding could offer valuable insights for\nmore accurately predicting future shifts in alpine endemism diversity under a rapidly\nwarming world.” (Lines 50-59 on Pages 3-4)\n[Comment 3] While I acknowledge the argument that \"If the drivers of spatial\npatterns of β-diversity and its two components also drive their temporal changes, the\nenhanced paleo-range limit fluctuation, coupled with paleo-temperature fluctuation,\ncould shift the pattern of dissimilarity derived from spatial replacement to a pattern of\ndissimilarity derived from nestedness, which has important implications for\nconservation planning\". But this is a local effect as you mention that \"shifts in upper\nrange limit represent significant changes in available habitat for alpine endemic\nspecies,\" indicating that these are local rather than regional processes. Therefore, my\nquestion remains: Why would this approach apply to differences in county-level\ncomposition when examining processes that define species pools? Said differently,\nwhy will local tree line changes remove or add endemic species from the city-level\nspecies pool?\nThis issue is not a deal breaker for me, and I would be happy to leave it as a point of\nacademic debate as I see this as a small point in the larger framework of this study.\n[Response] We strongly agreed with the reviewer that this issue is not a deal breaker,\nand decided to delete this sentence in the revised manuscript."}, "author_response": ""}]

Decision letter