Divergent changes in the elevational gradient of vegetation activities over the last 30 years

The reported progressive change of vegetation activity along elevational gradients has important aesthetic and conservation values. With climate change, cooler locations are suggested to warm faster than warmer ones, raising concerns of a more homogenized landscape along the elevation. Here, we use global satellite data to investigate the spatio-temporal dynamics of the elevational gradient (EG) in vegetation greenness (NDVImax3), spring (SOS) and autumn phenology (EOS) during 1982–2015. Although we find clear geographical patterns of the EG in NDVImax3 and SOS, there are no prevalent trends of vegetation homogenization or phenology synchronization along elevational gradients. Possible mechanisms, including spatially heterogeneous temperature lapse rate changes, different vegetation sensitivities to climate change, and human disturbances, may play diverse roles across different regions. Our finding of mixed EG trends and no general rules controlling EG dynamics poses challenges for mitigating possible adverse impacts of climate change on mountainous biological diversity and ecosystem services.

V egetation activity often displays a clear elevational pattern 1,2 . Such an elevational pattern of vegetation activity can refer to vegetation phenology. For example, the dates of spring leaf unfolding and autumn leaf senescence have been observed to progressively change from low to high elevations for many regions 3,4 . It can also be related to vegetation growth and greenness 5 . At many places, vegetation growth and greenness have been found to change progressively with elevation as a result of differences in climate or in nutrient availability 6 (relative to an optimum). Spatially, progressively changing vegetation activity along elevational gradients have been reported for many temperate and boreal regions 7 , with some cases also observed for the tropics 8,9 . This progressive vegetation activity change creates fascinating natural beauties that have been subjects of admirations in poems and proses for many generations. For instance, the great Chinese poet Bai Juyi once famously wrote, "In the plains past April, flowers have all but gone; In the hills at the temple, 'tis the time for the peach blossoms to bloom". It also helps maintain a diverse landscape of ecosystem structures and functions, particularly in mountainous areas that harbor the highest terrestrial biodiversity and provide critical ecosystem services to society.
Elevation-dependent temperature differences are often regarded as the main cause for the observed differences in vegetation activity along elevational gradients 5 . With climate change, however, there is an increasing concern that this elevationally progressive pattern of vegetation activity may be altered [10][11][12] . Specifically, the warming rate is often faster at higher than at lower elevations, leading to a trend of temperature homogenization 13 . This temperature homogenization has also been observed along latitudinal gradients, where it is believed to contributed to species and ecosystem homogenization 14 . The question thus arises whether, analogous to latitudinal temperature and ecosystem homogenization, elevation-dependent warming would also lead to more homogenized or synchronized vegetation activity? The elevational pattern of vegetation activity changes over the past decades under the reported significant climate change has to date not been studied. Because a more homogenized or synchronized vegetation pattern may reduce both the aesthetic aspects and functional diversities, addressing these questions is important for understanding and predicting vegetation activity and its ecosystem functioning under climate change.
Here, we use satellite-derived vegetation activity data to investigate how the relationship between vegetation activity and elevation may have changed over the past three decades. We consider three vegetation activity indicators that are derived from the Advanced Very High Resolution Radiometer (AVHRR) GIMMS NDVI 3g dataset over 1982-2015 at the spatial resolution of 8 km × 8 km to represent vegetation activity strength and phenology: the maximum three-monthly normalized difference vegetation index (NDVI max3 ), the start date of the growing season (SOS), and the end date of the growing season (EOS). These three vegetation activity indices provide key information for understanding vegetation growth, carbon cycles, and vegetation phenology, and have been found to show clear elevational patterns in many regions 15,16 . Here the relationship between vegetation activity and elevation is quantified as an elevational gradient (EG), which is defined as the average change of a vegetation activity index per 100 m increase in elevation. For each 8 km × 8 km grid cell, we estimate multi-year EG of vegetation activities by regressing the multi-year averaged vegetation activity indicators against its elevation in a 9 × 9 moving window centered at the focused grid cell (see Methods). Similarly, we also estimate EG for each year with annual vegetation activity data, and then derive the temporal EG trend using least squares linear regressions between annual EG and the year (see Methods). Note in this study, we use the average value from four different phenology derivative algorithms as the estimated SOS or EOS (see Methods). To test the robustness of the results, we also use other moving window sizes, i.e., 5 × 5,7 × 7,11 × 11,13 × 13 (see Methods). The aim of this study is three-fold: (1) to investigate the EG of the above three vegetation activity indicators during 1982-2015 and their spatial patterns; (2) to test if the temporal evolution of these elevational gradients follows a hypothetic trend of increased homogenization or synchronization; and (3) to identify possible mechanisms underlying the observed spatio-temporal dynamics of vegetation activity EG.

Results
Multi-year averaged EG of vegetation activities. We first explored the spatial patterns of the multi-year averaged EG values for each of the three vegetation activity indices during 1982-2015 (Fig. 1). The results showed highly heterogeneous, rather than consistent patterns of EG values across the global land surface for all three vegetation activity indices. The EG of NDVI max3 (EG ndvimax ) was negative (decreased vegetation greenness at higher elevations) in about 38% of the study area, mostly in cold regions including the Tibetan Plateau and the Artic regions (Fig. 1a). By contrast, it was positive (increased vegetation greenness at higher elevations) in the tropics and in most temperate regions, such as south-eastern America, Central Europe, and Eastern China (Fig. 1a). Furthermore, EG ndvimax calculated with other moving window sizes also showed similar patterns ( Supplementary Fig. 1), proving the robustness of the spatial pattern of EG ndvimax shown in Fig. 1a.
The highly varied EG ndvimax pattern may be caused by regional differences in the dominant mechanisms governing vegetation greenness. In cold regions, vegetation growth is mostly limited by low temperature 17 , and decreased vegetation greenness with increasing elevation (negative EG ndvimax , Fig. 1a) is therefore probably related to the cooler temperature at higher elevations (negative EG in temperature (EG tem ), Fig. 2a and Supplementary  Fig. 2a). In the tropics, where vegetation often grows near its thermal optimum, temperature may not be such a dominant limiting factor for vegetation greenness 18,19 . Instead, in the tropics the lower temperatures at increasing elevation may actually reduce evapotranspiration demand and respiratory carbon losses, resulting in enhanced vegetation growth and greenness at higher elevations. Furthermore, in many temperate regions, the increase of precipitation ( Fig. 2b and Supplementary  Fig. 2b) may have contributed to the observed positive EG ndvimax (Fig. 1a). Numerous field studies 20,21 have suggested that increasing precipitation along the elevation facilitated vegetation growth at higher elevations. The positive EG ndvimax may also be partly attributed to the reduction of human disturbances at higher elevations (Fig. 2c, d). For instance, the impact of reduced human activities at higher elevations, including reduced thinning intensity and less pollution 22 , can be particularly important for determining the sign of EG ndvimax in densely populated regions like Eastern North America, Europe, Eastern and Southern China.
The spatial distribution of the elevation gradient of SOS (EG SOS ) was also very heterogeneous, although positive EG SOS (later start of growing season at higher elevations) dominated in most of the study area (about 68%, Fig. 1b). The largest EG SOS appeared in Northern Europe and in south-eastern America, where SOS dates were often delayed by more than 5 days when climbing up every 100 m (Fig. 1b). By contrast, negative EG SOS values were found in northeast North America and central Europe (Fig. 1b). In addition, analyses with each of the four methods in deriving SOS from satellite imagery ( Supplementary  Fig. 3), rather than using the four-method averaged value, all showed similar spatial patterns in EG SOS as those shown in Fig. 1b. The corroborated results with different methods further confirmed the robustness of our finding about the spatial distribution of EG SOS signs.
The positive EG SOS over most regions is very likely attributable to the decreased temperature with elevation 3 , i.e. negative EG tem . As shown in Supplementary Fig. 4a and 4b, spring temperature decreases with increasing elevation in most regions. In northeast North America, however, both WorldClim 23 and ERA Interim temperature datasets 24 show increasing spring temperature with elevation (positive spring EG tem ), which explains the observed negative EG SOS there.
Unlike EG SOS , the EG of EOS (EG EOS ) showed a more fragmented spatial pattern (Fig. 1c). Positive EG EOS values (later autumn senescence at higher elevations) were found in southeast Europe and western North America, while negative ones (earlier autumn senescence at higher elevations) were found in high latitude regions and eastern Asia. Also in contrast to EG SOS , the spatial patterns of EG EOS derived with different methods did not always corroborate ( Supplementary Fig. 5). In particular, the pattern derived with the piecewise logistic method was significantly different from that obtained with the other three methods (i.e., HANTS, polyfit, and double logistic, see Methods). Further analyses showed that autumn NDVI increased with elevation in most regions ( Supplementary Fig. 6), and that the spatial pattern of the elevational gradient of autumn NDVI (EG ndviau ) was consistent with that of EG EOS based on the algorithms HANTS, polyfit, and double logistic, being mostly (56%) positive across the study area. The evidence from EG ndviau therefore appears to support the elevational EOS changes derived with the algorithms of HANTS, polyfit, and double logistic.
Like spring temperature, autumn EG tem was also mostly negative ( Supplementary Fig. 7a, c), inconsistent with the highly mixed spatial distribution of EG EOS . This spatial mismatch between autumn EG tem and EG EOS suggests that autumn temperature may not be the main environmental factor controlling autumn phenology. It has been suggested that radiation, precipitation, winter warming and even spring phenology [25][26][27] can all impact autumn phenology. We find that autumn precipitation amplifies with elevation in most regions ( Supplementary Fig. 7b, d), which may contribute to the postponed autumn phenology and enhanced autumn NDVI at higher elevations in some regions. Furthermore, the temperature sensitivity of autumn phenology is smaller than that of spring phenology 17 , which may also help explain the more fragmented spatial pattern of EG EOS . Nonetheless, so far our knowledge on the major mechanisms shaping the spatio-temporal pattern of autumn phenology remains very limited [28][29][30] .
Temporal dynamics of EG of vegetation activities. Next, we examined the temporal trends of the annual EG of vegetation activities from 1982 to 2015. To facilitate hypothesis testing, we calculated for each grid cell both the multi-year averaged EG sign and the temporal change in EG. Along the elevational gradient, the multi-year averaged EG sign can either (1) indicate increasing vegetation activity (positive EG ndvimax , negative EG SOS , and positive EG EOS ) or (2) suggest decreasing vegetation activity (negative EG ndvimax , positive EG SOS , and negative EG EOS ). Over For vegetation greenness (EG ndvimax ), we observed that about 48% of the study area falls into categories (1A) and (2A) (Fig. 4a). Spatially, the temporal trend of EG ndvimax was very fragmented (Fig. 4a). Category (1A) was mainly found in southeast Europe, Eastern European highlands, southern Africa and the Rocky Mountains. By contrast, category (1B) was mostly observed in western India, Southeast China, eastern Australia, and western Europe. Category (2A) mainly occurred in the Siberian boreal region; while category (2B) was mainly found in Alaska and the Yukon region. Overall, this result of fragmented spatial distribution of the four categories of EG ndvimax trends suggests that the elevational difference in vegetation greenness is not prevalently reduced and that the prediction of increased elevational homogenization of vegetation greenness is not supported.
The geographic distributions of the temporal trends in elevational temperature and precipitation gradients (EG tem or the elevational lapse rate of temperature, and EG pre , respectively) were also mixed ( Supplementary Fig. 8a, b). However, the spatial pattern of the sign of the EG tem or EG pre trends did not well match the pattern of the EG ndvimax trends (comparing Fig. 4a with Supplementary Fig. 8a, b), indicating that changes in EG tem or EG pre may not have been the dominant reason for the observed EG ndvimax trends at the global scale. Nonetheless, EG tem or EG pre changes may have explained some EG ndvimax trends at the regional scale. For example, the reduced negative EG ndvimax trend in western Europe may have resulted from the associated negative EG tem trend, while the same EG ndvimax trend in eastern Australia may be attributable to the negative EG pre trend. Human activities may also have contributed considerably to the EG trends in some specific regions ( Supplementary Fig. 8c, d). Fast population growth can imply increased anthropogenic disturbances and more frequent thinning (reduced vegetation greenness). For example, in central China, faster population growth at higher elevations may be responsible for the observed EG ndvimax decrease, while in western Europe, faster population growth at lower elevations is likely one reason leading to the increased EG ndvimax (Supplementary Fig. 8c, d). In addition, anthropogenic activities such as afforestation may also enhance vegetation activities. For instance, substantial afforestation has been implemented in regions like eastern China since the 1980s 31,32 , mainly at lower elevations ( Supplementary Fig. 9). As a result, the increasing afforestation at lower elevations leads to a negative temporal trend of EG ndvimax .
We further investigated the changes of EG in spring and autumn phenology (SOS and EOS) using the same method. For both EG SOS and EG EOS we found that only about half of the study area fell into categories (1B) and (2B) that represent increased elevational synchronization of spring and autumn phenology (Fig. 4b, c). Therefore, these results also did not support the concept of global warming-induced synchronization at the global scale.
Interestingly, although most of the spatial pattern of EG SOS could be explained by that of spring EG tem , the trend in spring EG tem failed to explain the EG SOS trend ( Supplementary  Fig. 10a, b). Further analyses suggested that the EG SOS trend was also related to elevational gradient changes in the temperature   Fig. 10d). On the other hand, in only about 20% of the area where the trend in EG SOS was controlled by EG tem ( Supplementary  Fig. 10d). Hence, the spatial pattern of EG SOS trends was more dominated by changes in the temperature sensitivity of spring phenology instead of by changes in spring temperature. This important role of spring phenology sensitivity to temperature in shaping the spatial pattern of spring phenology was previously reported 33,34 .

Discussion
We examined the spatio-temporal patterns of the elevational gradients of three vegetation activity indices representing vegetation greenness and phenology since 1982, and found a clear geographical pattern in the EG of vegetation greenness and spring phenology, but not in the EG of autumn phenology. Elevational temperature lapse was the dominant factor determining the elevational changes of spring phenology. Temperature was also important in shaping the EG of vegetation greenness; however, its significance and even the regulating mechanisms for the EG of vegetation greenness in tropical vs. boreal mountains were very different. The lack of explicit geographical pattern for the elevational gradient of autumn phenology and the spatial mismatch between the elevational gradients of autumn temperature and that of autumn phenology suggested that temperature may not be the main controlling factor for the EG of autumn phenology at the global scale. Little is known about the mechanisms controlling autumn phenology, especially when compared to our understanding of spring phenology. Further work on understanding mountainous autumn phenology is thus critically needed. Overall, our findings on the spatial patterns of the EG of vegetation activities suggest that temperature is not the sole dominant factor shaping the vegetation-elevation relationship, and highlights the importance of understanding vegetation activities and their drivers in the spatial context. More importantly, we also found highly variable temporal changes for the elevational gradients of all vegetation activities (Fig. 4), in contrast to the hypothesized increased vegetation greenness homogenization or vegetation phenology synchronization along the elevational gradient in response to climate change. Furthermore, previous concerns of more homogenized/ synchronized vegetation activities were based on the findings of enhanced warming rate at higher elevations 13,35-37 and the assumption of the dominant role of temperature in determining vegetation activities 5,14,36 . Nonetheless, Zeng et al. 37 suggested that warming is not always occurring faster at higher elevations. While vegetation activities can be more sensitive to climate change at higher elevations 3,5,38 , here we also find that mechanisms governing the spatio-temporal dynamics of the elevational gradients of vegetation greenness and phenology are not globally universal. Elevation-dependent warming did indeed lead to a reduction in the difference of spring phenology timing across the elevation gradient in some regions like the Alps 39 , yet temperature is not the sole element regulating the vegetation-elevation relationship 40 Our results are subject to a few caveats. First, we used climate data from the Worldclim dataset, which was interpolated from meteorological station records. Considerable errors in deriving the gridded climate values in the dataset, particular in high latitude areas, may arise from the sparse distribution of meteorological stations. For example, the uniform pattern of EG pre and EG tem in high latitudes may be partly caused by the fact that climate in the surrounding grids was likely interpolated from the same station data (Fig. 2a, b). More meteorological stations are thus needed in data-sparse areas, especially at higher latitudes.
Second, the GIMMS NDVI 3g dataset used in this study is at a spatial resolution of 8 km × 8 km, which can be rough for mountainous areas. However, this GIMMS NDVI 3g dataset is still likely the best one we can currently use for a global-scale analysis. For example, although the MODIS NDVI or EVI datasets has a higher spatial resolution, they cover a relatively short time period only since 2000. Future development of higher-resolution satellite datasets may further improve the robustness of such global-scale analyses of vegetation activity dynamics. To test the robustness of our satellite-based findings on the spatio-temporal EG trends of vegetation phenology, we perform the same analysis with a pan European phenological database PEP725, the largest ground-based phenological monitoring network. We find that the spatial pattern of ground-based EG SOS is similar to that of satellite-based, i.e., later leaf unfolded date with increasing elevation (Fig. 5a, b); and the area fraction of each of the four different spatio-temporal EG trend categories for ground-based EG SOS is also consistent with that of satellite-based EG SOS (Fig. 5c, d). This again confirms the robustness of our finding based on large-scale satellite data analysis. However, for EOS, most regions show negative correlations between leaf coloration date and elevation (sooner start of leaf coloring at higher elevations) in station-based EG EOS , a pattern in weak agreement with that from satellite data (Fig. 6a, b). Such discrepancy between satellite-and ground-based EG EOS spatial patterns in central Europe may partly arise from the sparse distribution of ground-based phenological stations, where only about 14.4 ± 5.8 sites have records for leaf unfold date and 12.1 ± 4.4 sites for leaf coloration for one 9 × 9 moving window. Furthermore, EOS from ground-based measurements is defined as the date with 50% leaf coloration, which may not be the same as satellite-based EOS derived from different algorithms detecting the abrupt change of NDVI. Nevertheless, the ground-based result still supports our major finding of no clear sign for increased elevational autumn phenology synchronization (Fig. 6c, d).
While our results do not support a prevalent trend of the elevational homogenization/synchronization of vegetation activities, such trends of recently increased homogenization/synchronization   (Fig. 4). Vegetation homogenization may have critical implications for landscape-scale diversity and ecosystem function 46 . For example, the increased vegetation phenological synchronization may reshuffle species composition and reduce species functional diversity 47 . Mountain areas harbor the most diverse terrestrial ecosystems and provide key ecosystem services, such as aesthetic appreciation, carbon storage, water conservation, and so on 48,49 . Increased vegetation homogenization and reduced species functional diversity along elevational gradients may impair the capacity of mountain ecosystems in providing these key services. Hence, for regions that do experience elevational vegetation homogenization, there is a great need to adopt appropriate policy tools to mitigate the adverse ecosystem impacts of such vegetation homogenization. These policy tools should be based on findings of region specific mechanisms driving the trend of vegetation homogenization.
Nonetheless, for most regions, we did not find apparent homogenization/synchronization trends in the EG of vegetation activities, nor do we have a clear answer on what drives them. More importantly, where a general rule governing the vegetation EG dynamics is lacking, even for regions where the vegetation EG trend was found to be statistically correlated with one or another environmental factor during the study period of 1982-2015, it remains highly uncertain how the vegetation EG will change in the future. This uncertainty imposes great challenges in predicting future landscape vegetation homogenization or phenology synchronization. It subsequently also casts uncertainty in best options for managing ecosystems to mitigate possible adverse impacts of climate change on mountainous biological diversity and ecosystem services.

Methods
Satellite-sensed vegetation index data. Normalized Difference Vegetation Index (NDVI) data is widely used as a proxy of vegetation greenness and photosynthetic activity 50,51 . Here we used the third generation Global Inventory Monitoring and Modeling Studies (GIMMS) group NDVI dataset 52 (referred as NDVI 3g ). The dataset has a 15-day temporal frequency from 1982 to 2015 with a spatial resolution of 8 km. We removed cropland-dominated grids from the analysis for their intense human management. Specifically, grids with more than 80% coverage by croplands according to the MODIS MCD12C1 land cover classification 53 were excluded. Note that this study excluded regions with bare soil/sparse vegetation where multi-year average NDVI is below 0.1.
Climate dataset. Monthly average temperature and precipitation data used in this study was acquired from the Worldclim version 2 dataset 23 . The Worldclim dataset was generated by interpolating the average values of climatic records from 9000-60000 weather stations using thin-plate splines for the period 1970-2000, with a high spatial resolution of about 1 km × 1 km. We also used the reanalysis dataset ERA Interim dataset 24 from the European Centre for Medium-Range Weather Forecasts (ECMWF), which provides monthly air temperature and precipitation data at a spatial resolution of 12 km for the period 1979-2018.
Human activity dataset. The global datasets of population density and night light can represent the strength of human activities. Here the population density data were obtained from the Gridded Population of the World version 4 data (GPW v4) from the Socioeconomic Data and Applications Center at Columbia University. Population density for the years of 2000, 2005, 2010, and 2015 from GPW v4 was estimated at a spatial resolution of 1 km × 1 km, using an annualized growth rate extrapolation method based on the 2010 round of Population and Housing Censuses 54   Ground-based phenology observation dataset. The ground-based phenology observation data are obtained across Central Europe from the Pan European Phenological database (PEP725). We seleted the records of the first leaf unfolding date and 50% leaf coloration date from 2256 sites and focused on the sites with continuous records for 1982-2011.
Vegetation activity indicators. The maximum three-monthly NDVI max3 was calculated following the below steps: (i) calculating the mean monthly NDVI value for the period 1982-2015, and (ii) sorting the twelve average-monthly NDVI values to obtain the first three maximum values.
We also derived two phenological indices from the GIMMS NDVI 3g dataset, i.e. the SOS and the EOS. The SOS date was calculated as the average SOS date estimated by the below four methods: HANTS (harmonic analysis of time series), polyfit, timesat and spline [56][57][58] . The EOS date was obtained from the averaged EOS date estimated by the following four methods: HANTS, double logistic, piecewise logistic, and polyfit 15,57,[59][60][61] .
The calculation of elevation gradients. Before analyses, all datasets were resampled into a unified spatial resolution of 8 km (the same as the NDVI data). The elevation gradient of NDVI max3 (EG ndvimax ) was calculated as the linear regression slope of the 34-year-averaged NDVI max3 against elevation for each 9 × 9 moving window. Digital Elevation Model (DEM) data 62 was obtained from the Shuttle Radar Topographic Mission (SRTM) with a spatial resolution of 30 m. Similarly, the elevation gradient of SOS/EOS (EG SOS /EG EOS ) was calculated as the linear regression slope of the multi-year-averaged of SOS/EOS against elevation within each 9 × 9 moving window. The elevation gradients of temperature (EG tem ), total precipitation (EG pre ), population density (EG pd ), and night light index (EG NLI ) were also calculated using the same methods. A positive elevation gradient indicates that the value of focused index is higher at higher elevations, while a negative one suggests that the value is higher at lower elevations. The statistical significance of the elevation gradient was determined by 95% confidence intervals. Note that only 9 × 9 windows with more than 40 grid cells of valid vegetation index values and more than 50 m in elevation difference were chosen for the regression analysis. Similar analyses were also performed using 5 × 5, 7 × 7, 11 × 11 or 13 × 13 moving windows to test the robustness of the results. Note that only 9 × 9 windows with more than 40 grid cells of valid vegetation index values and more than 50 meters in elevation difference were chosen for the regression analysis. Based on the Pan European Phenological database (PEP725), We calculate the elevation gradient of SOS (EOS) as the linear regression slope of leaf unfolded date (leaf coloration date) for each 9 × 9 moving window with more than five stations the same way as we do for satellite-based analyses.
Driving factors of temporal trends in EG. To investigate possible effects of climate change and human activities on the temporal trends in vegetation activity EG, for each pixel we obtained the linear trends of EG in vegetation activity, climate, and night light index during 1982-2015 based on least-squares regressions. Note the data of night light index is available for 1982-2013; and because the number of available population density data for each grid during the study period is only 4, we therefore calculated the elevational gradient of the population density difference between 2015 and 2000 (EG pd_d ), rather than estimating the trend in EG pd .
Note EG of vegetation phenology can also be influenced by phenology temperature sensitivity 34 . We therefore also calculated the elevational gradient of the temperature sensitivity of vegetation phenology. The elevational gradient of the temperature sensitivity of SOS was calculated as the following steps. (i) For each pixel, the preseason was determined as the period before SOS for which the absolute value of the correlation coefficient between SOS and air temperature was highest 63 . Note that since the relationship between SOS and temperature is usually regarded as negative, respectively, we thus exclude those positive correlation coefficients for SOS. (ii) For each 9 × 9 moving spatial window, we used an identical preseason length, which is the median values of preseason length of each grid. (iii) The temperature sensitivity of SOS was calculated using a linear regression for the period 1982-2011 with temperature as the dependent variable and SOS as the independent variable. (iv) The elevational gradient of temperature sensitivity of SOS was estimated for each 9 × 9 moving window. To examine the dominant factor driving the temporal change of EG SOS , we classified all grid into another four   Table 1) based on the sign of EG SOS trends, spring EG tem trends, and elevation gradients of temperature sensitivity of SOS: (a) the EG SOS trend is co-determined by both EG tem and the temperature sensitivity of SOS, (b) the EG SOS trend is primarily driven by the temperature sensitivity of SOS, (c) the EG SOS trend is primarily driven by spring EG tem , (d) the EG SOS trend is possibly driven by other factors.
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.