Introduction

Vegetation activity often displays a clear elevational pattern1,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 regions3,4. It can also be related to vegetation growth and greenness5. 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 availability6 (relative to an optimum). Spatially, progressively changing vegetation activity along elevational gradients have been reported for many temperate and boreal regions7, with some cases also observed for the tropics8,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 gradients5. With climate change, however, there is an increasing concern that this elevationally progressive pattern of vegetation activity may be altered10,11,12. Specifically, the warming rate is often faster at higher than at lower elevations, leading to a trend of temperature homogenization13. This temperature homogenization has also been observed along latitudinal gradients, where it is believed to contributed to species and ecosystem homogenization14. 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 NDVI3g 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 (NDVImax3), 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 regions15,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 NDVImax3 (EGndvimax) 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, EGndvimax calculated with other moving window sizes also showed similar patterns (Supplementary Fig. 1), proving the robustness of the spatial pattern of EGndvimax shown in Fig. 1a.

Fig. 1
figure 1

The elevational gradient (EG) of vegetation activity indices. a The spatial pattern of the EG of NDVImax3 (EGndvimax) from 1982 to 2015. b the spatial pattern of the EG of SOS (EGSOS) from 1982–2011 over the Northern Hemisphere (north of 30°N). c The spatial pattern of the EG of EOS (EGEOS) from 1982 to 2011 over the Northern Hemisphere. A frequency distribution of the EG values is shown in the inset at the bottom-left for each panel. Regions marked with dots have statistically significant (p < 0.05 with the F-test) EG values

Full size image

The highly varied EGndvimax pattern may be caused by regional differences in the dominant mechanisms governing vegetation greenness. In cold regions, vegetation growth is mostly limited by low temperature17, and decreased vegetation greenness with increasing elevation (negative EGndvimax, Fig. 1a) is therefore probably related to the cooler temperature at higher elevations (negative EG in temperature (EGtem), 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 greenness18,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 EGndvimax (Fig. 1a). Numerous field studies20,21 have suggested that increasing precipitation along the elevation facilitated vegetation growth at higher elevations. The positive EGndvimax 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 pollution22, can be particularly important for determining the sign of EGndvimax in densely populated regions like Eastern North America, Europe, Eastern and Southern China.

Fig. 2
figure 2

The elevational gradient (EG) of climate and human activity variables. a The spatial pattern of the EG of mean annual temperature (EGtem). b The spatial pattern of the EG of total precipitation (EGpre). c The spatial pattern of the EG of population density (EGpd). d The spatial pattern of the EG of night light index (EGNLI). A frequency distribution of the EG values is shown in the inset at the bottom-left for each panel. Regions marked with dots have statistically significant (p < 0.05 with the F-test) EG values

Full size image

The spatial distribution of the elevation gradient of SOS (EGSOS) was also very heterogeneous, although positive EGSOS (later start of growing season at higher elevations) dominated in most of the study area (about 68%, Fig. 1b). The largest EGSOS 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 EGSOS 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 EGSOS as those shown in Fig. 1b. The corroborated results with different methods further confirmed the robustness of our finding about the spatial distribution of EGSOS signs.

The positive EGSOS over most regions is very likely attributable to the decreased temperature with elevation3, i.e. negative EGtem. As shown in Supplementary Fig. 4a and 4b, spring temperature decreases with increasing elevation in most regions. In northeast North America, however, both WorldClim23 and ERA Interim temperature datasets24 show increasing spring temperature with elevation (positive spring EGtem), which explains the observed negative EGSOS there.

Unlike EGSOS, the EG of EOS (EGEOS) showed a more fragmented spatial pattern (Fig. 1c). Positive EGEOS 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 EGSOS, the spatial patterns of EGEOS 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 (EGndviau) was consistent with that of EGEOS based on the algorithms HANTS, polyfit, and double logistic, being mostly (56%) positive across the study area. The evidence from EGndviau therefore appears to support the elevational EOS changes derived with the algorithms of HANTS, polyfit, and double logistic.

Like spring temperature, autumn EGtem was also mostly negative (Supplementary Fig. 7a, c), inconsistent with the highly mixed spatial distribution of EGEOS. This spatial mismatch between autumn EGtem and EGEOS 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 phenology25,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 phenology17, which may also help explain the more fragmented spatial pattern of EGEOS. Nonetheless, so far our knowledge on the major mechanisms shaping the spatio–temporal pattern of autumn phenology remains very limited28,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 EGndvimax, negative EGSOS, and positive EGEOS) or (2) suggest decreasing vegetation activity (negative EGndvimax, positive EGSOS, and negative EGEOS). Over time, the EG sign can be either (A) amplified or (B) reduced. The combined changes in EG along the elevational gradient and over time are illuminated in Fig. 3. Here (1A) and (2A) indicate increased elevational vegetation dissimilarity; while (1B) and (2B) imply reduced vegetation dissimilarity, i.e. increased elevational vegetation homogenization/synchronization.

Fig. 3
figure 3

A diagram of changes in vegetation homogenization along elevational gradients. The vegetation activity here is showcased with the example of the maximum three-monthly normalized difference vegetation index (NDVImax3). Along the elevation, the multi-year averaged elevational gradient (EG) sign can be either (1) increasing vegetation activity (positive EGndvimax) or (2) decreasing vegetation activity (negative EGndvimax). Over time, the EG sign can either (A) amplified or (B) reduced. The combined changes along the elevation and over time for EG dynamics are illuminated in Fig. 3. Here (1A, a) and (2A, c) indicate decreased elevational vegetation homogenization; while (1B, b) and (2B, d) imply increased vegetation homogenization

Full size image

For vegetation greenness (EGndvimax), we observed that about 48% of the study area falls into categories (1A) and (2A) (Fig. 4a). Spatially, the temporal trend of EGndvimax 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 EGndvimax 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.

Fig. 4
figure 4

Spatial patterns of different temporal trends of the elevational gradient (EG) of vegetation activities. a The spatial pattern of different categories of the temporal trends in EGndvimax during 1982–2015. b The spatial pattern of different categories of the temporal trends in EGSOS during 1982–2011 over the Northern Hemisphere (north of 30 °N). c The spatial pattern of different categories of the temporal trends in EGEOS during 1982–2011 over the Northern Hemisphere (north of 30 °N). Category colors in this figure correspond to that shown in Fig. 3. Regions marked with dots have statistically significant (p < 0.05 with the F-test) EG values. The area fraction of each of the four different categories in the temporal trends of EGndvimax/EGSOS/EGEOS is shown in the inset at bottom-left for each panel

Full size image

The geographic distributions of the temporal trends in elevational temperature and precipitation gradients (EGtem or the elevational lapse rate of temperature, and EGpre, respectively) were also mixed (Supplementary Fig. 8a, b). However, the spatial pattern of the sign of the EGtem or EGpre trends did not well match the pattern of the EGndvimax trends (comparing Fig. 4a with Supplementary Fig. 8a, b), indicating that changes in EGtem or EGpre may not have been the dominant reason for the observed EGndvimax trends at the global scale. Nonetheless, EGtem or EGpre changes may have explained some EGndvimax trends at the regional scale. For example, the reduced negative EGndvimax trend in western Europe may have resulted from the associated negative EGtem trend, while the same EGndvimax trend in eastern Australia may be attributable to the negative EGpre 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 EGndvimax decrease, while in western Europe, faster population growth at lower elevations is likely one reason leading to the increased EGndvimax (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 1980s31,32, mainly at lower elevations (Supplementary Fig. 9). As a result, the increasing afforestation at lower elevations leads to a negative temporal trend of EGndvimax.

We further investigated the changes of EG in spring and autumn phenology (SOS and EOS) using the same method. For both EGSOS and EGEOS 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 EGSOS could be explained by that of spring EGtem, the trend in spring EGtem failed to explain the EGSOS trend (Supplementary Fig. 10a, b). Further analyses suggested that the EGSOS trend was also related to elevational gradient changes in the temperature sensitivity (St) of spring phenology (Supplementary Fig. 10c, see Methods). Based on the signs of the EGSOS trend, the spring EGtem trend, and the elevational gradient changes in St (EGSt) of spring phenology, we also grouped all grids into another four categories according to the dominance of EGtem and EGSt in controlling the EGSOS trend (see Methods). We found that for about 31% of the Northern Hemisphere (NH) land area, the trend in EGSOS was controlled by EGSt, and also for about 31% of the NH, it was controlled by both EGSt and EGtem (Supplementary Fig. 10d). On the other hand, in only about 20% of the area where the trend in EGSOS was controlled by EGtem (Supplementary Fig. 10d). Hence, the spatial pattern of EGSOS 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 reported33,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 elevations13,35,36,37 and the assumption of the dominant role of temperature in determining vegetation activities5,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 elevations3,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 Alps39, yet temperature is not the sole element regulating the vegetation-elevation relationship40,41. Precipitation and human activities, for instance, have both played important roles in determining the elevational gradients of vegetation activities and their changes40,42,43,44,45, with various contributions in different areas. Furthermore, these multiple factors that control the vegetation-elevation relationship, as well as their recent dynamics, are highly heterogeneous across the global land surface. Such heterogeneities probably underpin the highly variable responses of the vegetation-elevation relationship to recent climate change. These variable vegetation-elevation responses under climate change need to be further investigated for better understanding and improved prediction of future vegetation activities under different climate scenarios.

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 EGpre and EGtem 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 NDVI3g 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 NDVI3g 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 EGSOS 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 EGSOS is also consistent with that of satellite-based EGSOS (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 EGEOS, a pattern in weak agreement with that from satellite data (Fig. 6a, b). Such discrepancy between satellite- and ground-based EGEOS 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).

Fig. 5
figure 5

Spatial patterns of the elevational gradient of SOS (EGSOS) and their temporal trends. a The spatial pattern of ground-based EGSOS from 1982 to 2011. b The spatial pattern of satellite-based EGSOS from 1982 to 2011. Frequency distribution of the EGSOS values is shown in the inset at the bottom-right for each panel. c The spatial pattern of different categories of the temporal trends in ground-based EGSOS during 1982–2011. d The spatial pattern of different categories of the temporal trends in satellite-based EGSOS during 1982–2011. The area fraction of each of the four different categories in the temporal trends of EGSOS is shown in the inset at bottom-right for each panel. The location of the regions with European phenology stations is shown in the inset at top-right for each panel

Full size image
Fig. 6
figure 6

Spatial patterns of the elevational gradient of EOS (EGEOS) and their temporal trends. a The spatial pattern of ground-based EGEOS from 1982 to 2011. b The spatial pattern of satellite-based EGEOS from 1982 to 2011. Frequency distribution of the EGEOS values is shown in the inset at the bottom-right for each panel. c The spatial pattern of different categories of the temporal trends in ground-based EGEOS during 1982–2011. d The spatial pattern of different categories of the temporal trends in satellite-based EGEOS during 1982–2011. The area fraction of each of the four different categories in the temporal trends of EGEOS is shown in the inset at bottom-right for each panel. The location of the regions with European phenology stations is shown in the inset at top-right for each panel

Full size image

While our results do not support a prevalent trend of the elevational homogenization/synchronization of vegetation activities, such trends of recently increased homogenization/synchronization do appear in some regions (Fig. 4). Vegetation homogenization may have critical implications for landscape-scale diversity and ecosystem function46. For example, the increased vegetation phenological synchronization may reshuffle species composition and reduce species functional diversity47. Mountain areas harbor the most diverse terrestrial ecosystems and provide key ecosystem services, such as aesthetic appreciation, carbon storage, water conservation, and so on48,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 activity50,51. Here we used the third generation Global Inventory Monitoring and Modeling Studies (GIMMS) group NDVI dataset52 (referred as NDVI3g). 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 classification53 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 dataset23. 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 dataset24 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 Censuses54. Night light data with a spatial resolution of 30 arc-seconds (approximately 1 km) for the period 1992–2013, covering regions from 65 °S to 75 °N, were downloaded from the Defense Meteorological Satellite Program—Operational Linescan System (DMSP-OLS) night light dataset55, which was derived from the National Oceanic and Atmospheric Administration (NOAA) National Geographic Data Center (NGDC).

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 NDVImax3 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 NDVI3g 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 spline56,57,58. The EOS date was obtained from the averaged EOS date estimated by the following four methods: HANTS, double logistic, piecewise logistic, and polyfit15,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 NDVImax3 (EGndvimax) was calculated as the linear regression slope of the 34-year-averaged NDVImax3 against elevation for each 9 × 9 moving window. Digital Elevation Model (DEM) data62 was obtained from the Shuttle Radar Topographic Mission (SRTM) with a spatial resolution of 30 m. Similarly, the elevation gradient of SOS/EOS (EGSOS/EGEOS) 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 (EGtem), total precipitation (EGpre), population density (EGpd), and night light index (EGNLI) 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 (EGpd_d), rather than estimating the trend in EGpd.

Note EG of vegetation phenology can also be influenced by phenology temperature sensitivity34. 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 highest63. 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 EGSOS, we classified all grid into another four categories (categories a-d in Supplementary Table 1) based on the sign of EGSOS trends, spring EGtem trends, and elevation gradients of temperature sensitivity of SOS: (a) the EGSOS trend is co-determined by both EGtem and the temperature sensitivity of SOS, (b) the EGSOS trend is primarily driven by the temperature sensitivity of SOS, (c) the EGSOS trend is primarily driven by spring EGtem, (d) the EGSOS 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.