Environmental Drivers of NDVI-Based Vegetation Phenology in Central Asia

Through the application and use of geospatial data, this study aimed to detect and characterize some of the key environmental driv ers contributing to landscape-scale vegetation response patterns in Central Asia. The o bj ctives of the study were to identify the variables driving the year-to-year vegetation d y amics in three regional landscapes (desert, steppe, and mountainous); and to determine if th identified environmental drivers can be used to explain the spatial-temporal variabi lity of these spatio-temporal dynamics over time. It was posed that patterns of change in t rrestrial phenology, derived from the 8 km bi-weekly time series of Normalized Difference V getation Index (NDVI) data acquired by the Advanced Very High Resolution Radio meter (AVHRR) satellites (1981–2008), can be explained through a multi-scale analysis of a suite of environmental drivers. Multiple linear stepwise regression analys es were used to test the hypotheses and address the objectives of the study. The annually c omputed phenological response variables or pheno-metrics time (season start, season length, and an NDVI-based pr oductivity metric) were modeled as a function of ten environme tal factors relating to soil, topography, and climate. Each of the three studied regional landscapes was shown to be governed by a distinctive suite of environmental dr ivers. The phenological responses of the steppe landscapes were affected by the year-to-year variation in temperature regimes. The phenology of the mountainous landscapes was influen ced primarily by the elevation gradient. The phenological responses of desert land scapes were demonstrated to have the greatest variability over time and seemed to be aff cted by soil carbon content and OPEN ACCESS Remote Sens. 2011, 3 204 year-to-year variation of both temperature regimes and winter precipitation patterns. Amounts and scales of observed phenological variabi lity over time (measured through coefficient of variation for each pheno-metric time) in each of the regional landscapes were interpreted in terms of their resistance and resili ence capacities under existing and projected environmental settings.


Introduction
Environmental drivers like climate, topography and soil properties affect vegetation dynamics at different spatial and temporal scales, ranging from instant to long-term and from local to regional scales.This makes the assessment of their response and ecological processes complex because due to the spatial and temporal scale dependence of drivers, the landscape-scale responses that are the most substantial for a specific place and time might be less important in other space or time scales.The development and advancement of the field of landscape ecology plays an important role in the attempts to unravel complex and interrelated landscape-scale terrestrial response dynamics, their drivers, and consequences.Various digital biophysical datasets (e.g., vegetation indices) represent an essential part of environmental assessments and have become particularly important when studying the pace and extent of landscape change dynamics across space and time [1][2][3][4][5].
While much attention has been given to predictive modeling of climate scenarios at regional to global scales, the scales at which the interplay of environmental factors will impact landscape-scale vegetation response across time and space has received far less emphasis.A more nuanced understanding of this matter will be an important contribution to the field of global change research.Notably, such research can provide insight and interpretation of climate and environmental change in Central Asia, a region that has had limited attention from the international scientific community [6].Therefore, through the application of geospatial data, this study proposed to detect, document, and explain the environmental drivers contributing to landscape-scale vegetation response patterns in Central Asia, an arid region that occupies about four million km 2 in the heart of Eurasia.Specifically, the objectives of this research were to develop a spatially explicit model for various regional landscapes of Central Asia, to (1) identify the key environmental variables driving annual (the year-to-year) vegetation dynamics in each regional landscape and (2) determine if the identified environmental drivers can be used to explain the variability-the tendency for deviation in the response-detected in vegetation dynamics across time.Terrestrial vegetation is often viewed as the most overt evidence of biological response to climatic and other environmental factors [7,8], and because of this notion, land surface phenological variables were used in this study as a measure of landscape-scale response to environmental drivers [1,9,10].
The main hypothesis of this study is that the patterns of change in land surface phenology can be identified through a multi-scale analysis of environmental drivers and can be expanded into two related hypotheses.First, there is a set of particular environmental drivers that consistently drive vegetation dynamics from year to year in a given region.Second, regional-scale landscapes will demonstrate unique yet predictable combinations of spatially explicit (e.g., in accordance with latitudinal and/or altitudinal gradients) phenological responses to these environmental drivers, enabling such drivers to be used as landscape-scale predictors of variability in vegetation response over time.
Variability in response patterns is often associated with the concepts of stability in a given ecosystem [11,12], and measures of variability are used to describe characteristic features of sensitivity to change of that ecosystem.Ecosystem sensitivity to change, as a response to short-term perturbations and long-term stressors, is a function of environmental factors at various temporal and spatial scales [13,14], and depends on two components that make up the concept of ecosystem stability: the capacity to withstand perturbation while maintaining regular functioning (resistance), and the ability to recover structural and functional attributes from perturbation (resilience) [12,15,16].Some systems may exhibit a higher degree of resistance than resilience; moreover, the systems that demonstrate high resistance and high resilience are expected to be the most stable systems (Table 1).Hence, it can be argued that ecosystems which have low resistance and high resilience capabilities are more able to adapt to long-term environmental change than ecosystems with high resistance and low resilience capabilities, as the latter may lose their recovery and functioning abilities with long-term stress.For instance, healthy woodland systems of the mountainous regional landscape of Central Asia might be comparatively resistant to fire due to the protective bark structure and relatively high moisture content of the soil [17,18].But if these systems catch fire and burn due to prolonged regional droughts, their recovery rate (resilience capability) might be relatively low and the systems may never return to their initial states.On the other hand, the grassland vegetation of the dry deserts and steppe areas of Central Asia might not be ostensibly resistant to fire [19], but it might recover quickly, exhibiting high resilience qualities [20].
Table 1.Conceptual framework: resistance and resilience capacities of some ecological systems.

Resilience
Low High of this regional landscape are a general absence of trees, and a continental climate [22] with maximum precipitation in summer and minimal precipitation in spring and fall seasons [23].The steppe regional landscape's vegetation response patterns are assumed to be governed mostly by the climate regimes that follow a very clear latitudinal gradient [24].
Figure 1.Map of study area extent with areas of the three regional landscapes vegetation responses that were assessed through modeling of their pheno-metrics: steppe regional landscape outlined by the green color; mountainous regional landscape outlined by the purple color; and desert regional landscape outlined by orange color.The map also demonstrates locations of the rainfed and irrigated agricultural areas excluded from the analysis.The area of the mountainous regional landscape incorporates the territories of Tajikistan and Kyrgyzstan (Figure 1), covering about 135 thousand km 2 .The mountainous regional landscape represents a key focal point in the Central Asian landscape that experiences adverse effects of the rapid rates of climate change [25][26][27].The snow pack and glaciers in the high mountains of the Tien Shan and Pamir, located in this area, are the origins for numerous rivers meandering through the Central Asian terrain and constitute about 70% and 21% of the total fresh water resources of the Central Asian region, respectively [28].Intensified melting of the glaciers and snow pack occurs under conditions of warming climate [25] is expected to lead to a decrease of glacier areas by about 20% in the next two decades [27], causing temporary increases in runoff and river flow, overflowing and flooding of existing mountain lakes, and the increasing loss of melt water due to evaporation caused by warmer temperatures.This regional landscape is characterized by the highest amount of rainfall distribution in Central Asia, which peaks in late winter and early spring and totals 800 to 2,000 mm per year [29].The vegetation dynamics of the mountainous regional landscape are expected to be driven mostly by temperature thresholds along the altitudinal gradients [30] and exhibit patterns of clinal variation-a gradual phenotypic and/or genetic variation over a geographical area-similar to those observed across the latitudinal gradient of the steppe regional landscape [31].Similarity, altitudinal and latitudinal patterns can elicit responses to the same environmental factors, i.e., spatial clines in climate and energy [32].

Mountainous
The semi-desert and desert regional landscape comprises the territories of Turkmenistan and Uzbekistan (Figure 1) and will be referred to as the desert regional landscape in the study hereafter.This regional landscape covers approximately one million km 2 and is predominantly represented by the sandy Karakum and Kyzylkum deserts [33].Major freshwater resources for this regional landscape come from the perennial rivers that are fed by seasonal melt water from the snow packs and glaciers of the mountainous regional landscape.The vegetation dynamics of the desert regional landscape are assumed to be mostly precipitation-driven [30], as this regional landscape is characterized by the lowest amount of precipitation in Central Asia [27].It receives about 75-100 mm precipitation per year peaking in the spring [23].

Datasets
Response Variables: Deriving Metrics of NDVI-Based Vegetation Dynamics Annual phenological metrics (pheno-metrics) for the study sites were derived from the 1981-2008 time-series of biweekly (15 days) composited 8-km Normalized Difference Vegetation Index (NDVI) values of the Global Inventory Modeling and Mapping Studies (GIMMS) project acquired by the Advanced Very High Resolution Radiometer (AVHRR) sensors 7, 9, 11, 14, 16 and 17 on board the National Oceanic and Atmospheric Administration (NOAA) satellite platforms; the NDVI data were calibrated and corrected for view geometry, volcanic aerosols, and other miscellaneous issues that were not related to vegetation response [34,35].The source for this data set was the Global Land Cover Facility (www.landcover.org).These NDVI data are widely used for time series trend analysis [36,37], because of the efforts by the producers to minimize inconsistencies [34,35].The NDVI time-series dataset spans 28 years, 24 images per year, resulting in 672 gridded images with observations of vegetation dynamics.TIMESAT time-series analysis software was used to extract metrics of vegetation dynamics (phenological metrics) for each year ).An adaptive Savitsky−Golay smoothing filter [38] was used because it maintains distinctive vegetation time-series curves and minimizes various atmospheric effects [39,40].For the purpose of this study, two timing and one greenness metrics were used: (1) the start of growing season, the time when base NDVI value has increased by 20% [40] of the distance between the pre−season minimum and the seasonal maximum; (2) the length of growing season, the difference between the start and end dates of the growing season; and (3) the integrated NDVI (large integral), for the season.We used a consistent start of growing season NDVI threshold (20%; other thresholds were used with less success) to optimize our interannual phenological retrieval rate for pixels with a seasonal response, while maintaining a start of season that relates to the beginning of the increase in photosynthetic activity.The use of this threshold could introduce a bias (an earlier SOS would be expected with lower threshold and a later season start with a higher threshold), but is not likely impacting the relationship between response and explanatory variables.The timing based pheno-metrics are important indices that can characterize seasonal photosynthetic thresholds and patterns among different vegetation types and along latitudinal and/or altitudinal gradients.These metrics are also quite informative when considering the capacity of vegetation types to assimilate carbon during their growing season [41,42].Metrics integrating NDVI values over the growing season are often used as a proxy measurement for seasonal vegetation productivity [43], and in this study the large integral is used as an NDVI-based productivity metric.Each of the selected metrics was accordingly expected to reveal unique spatio-temporal response patterns to environmental drivers in each of the studied regional landscapes.

Explanatory Variables: Obtaining Key Environmental Drivers
The suite of potential explanatory environmental variables hypothesized as drivers of the vegetation dynamics measured in this study included climate, soil, and topography characteristics.To address potential multicollinearity, a common issue when interpreting outputs of multiple regression models, interdependency among all the explanatory variables was assessed to ensure that their interrelationships had lower coefficients of determination than R 2 = 0.64 (or correlation coefficient R = ±0.8)[44,45].Consequently, ten environmental variables were used, among which were eight seasonal precipitation and temperature-based variables (out of 10 climate variables considered), one topography-based factor (of 3 considered), and one soil-based variable (of 7 considered) (Table 2).
Climate-based data included the mean monthly precipitation (mm) and temperature (°C) data available from the Global Land Data Assimilation System (GLDAS) modeled dataset with one degree of spatial resolution [46].The climate dataset [46] represents synthetic data from the land surface models (uncoupled from the atmospheric models) forced with precipitation gauge observations, satellite data, radar precipitation measurements, and output from numerical prediction models [47].From the GLDAS data, seasonal temperature and precipitation data were calculated from 1980 to 2008, to coincide with the available temporal extent of the NDVI-based metrics.Four seasons were created for (cumulative) precipitation and (average) temperature variables: the winter season included December, January, and February (DJF); the spring season included March, April, and May (MAM); the summer season consisted of June, July, and August (JJA), and the fall season included September, October, and November (SON).Because fall season variables demonstrated high multicollinearity to other seasonal variables, fall temperature and precipitation were not included in the phenological modeling phase.However, antecedent fall temperature and precipitation did not demonstrate clear collinear behavior and therefore were incorporated in the modeling effort to account for possible lags in phenological responses.The dynamics of vegetation in a given ecosystem are influenced not only by the weather but also by topography and soil properties of the system.Elevation values (m) were derived from the World Digital Elevation Model (DEM) file, included in the sample datasets of ENVI software.The DEM data were 0.1 degree of spatial resolution.The global gridded surfaces of selected soil characteristics from the Global Soil Data Task Group International Geosphere−Biosphere Programme (IGBP−DIS) dataset were used to represent soil variables within the study sites [48].The carbon content variable showed the highest spatial variability after being tested for multicollinearity among seven soil-based metrics (soil carbon content, total nitrogen content, field capacity, wilting point, profile available water capacity, thermal capacity, and bulk density) and was selected to represent a soil-based explanatory metric.Calculated quantities of the carbon content at a depth of 0-100 cm were extrapolated and interpolated from point locations to values representing areas by linkage to soil maps to represent gridded carbon content values (kg/m 2 ) at five arc-minutes (0.083333 degree) spatial resolution [48].
All explanatory variables were resampled to 8 km to match the spatial resolution of the phenological response metrics.Climate-based variables were temporally dynamic explanatory metrics that were used to account for changes in vegetation responses over space and time.Soil-and topography-based drivers, stable over time, were used to characterize spatial variation of the response patterns.All the pixels identified as either irrigated or rainfed cropland (Figure 1) were removed from the analysis to avoid the effect of these land cover types on the evaluation of vegetation dynamics in non-cropland areas, the primary focus of this research.Cropland was identified using the Global Irrigated Agriculture Map (GIAM) dataset [49] and the Global Map of Rainfed Cropland Areas (GMRCA) dataset [50] respectively, available from the International Water Management Institute web portal [51].

Phenological Modeling Methods: Analyses and Testing Hypothesis
What Environmental Factors Drive Spatial and Temporal Variation in Vegetation Dynamics?
To address the study's first hypothesis, the year-to-year response of the three selected pheno-metrics to all ten environmental explanatory variables was assessed.For each phenol-metric, stepwise multivariate regression models were run for each of the 28 years  in each of the three studied regional landscapes separately, totaling 252 models run for this stage of the analysis.The annual variation in each of the pheno-metrics was modeled as a function of environmental factors at yearly time steps: where x identifies the suite of both stable (s) and dynamic (d) environmental factors and a = 28 years (Table 2).Each response metric was analyzed using a unique suite of environmental variables that was most relevant to the metric's spatio−temporal variation.Model results were examined to determine the most frequently appearing yearly drivers of vegetation dynamics in each of the regional landscapes, and to assess the causes of variation in these dynamics across time.Additionally, measures of explanatory power (adjusted R 2 ) were obtained for individual explanatory variables (simple linear regression model) that contributed significantly (p <= 0.05) to each of the models: This statistic was then used for a separate level contribution assessment of each of the ten environmental variables.

What Environmental Drivers Explain Phenological Variability (COV) across Time (28 yrs)?
To address the second hypothesis determining whether the identified environmental factors explain detected variability in vegetation dynamics across time, the COV, the ratio of the standard deviation and the mean, for each phenol-metric was used as a measure of variability in the vegetation dynamics over time [52,53].Higher COV values are associated with higher susceptibility to the changes in environmental drivers of a given system, and often are used to represent the vulnerability and resilience properties of that system [52,[54][55][56].Antle et al. [54] simulated the effects of climate change on land use and net returns to grain production systems by using the COV as a measure of vulnerability to climate change with and without adaptation in agricultural production systems, and argues that adaptive capacity is associated with lower COV values.Hurd et al. [55] studied the impact of climate change projections on key aspects of water supply and used COV as a measure of vulnerability, and states that high variability (high COV) values are indicative of the areas that are closer to system thresholds and are vulnerable to substantial adverse effects of climate change.Luers [57] used COV values as a linked measure of sensitivity and exposure of the system to characterize relative vulnerability of a variable of concern (agricultural yield) to a set of forces of disturbance.
In this study, COV values are viewed as measures of temporal variability of annual phenological responses, to represent the given systems' stability: where x s are the stable (soil and elevation) environmental factors and a d x , are the seasonal climate-based explanatory variables, which were averaged across the 28 years to avoid pseudoreplication.To assess the relative variability in the pheno-metrics' responses among the regional landscapes, the COV values were calculated using the inputs from the entire extent of Central Asia, not for each regional landscape separately.The explanatory drivers included in these series of models were the variables that contributed explanatory power to the multilinear regression models during at least 14 out of 28 years.In total, nine regression models were run for this step of the analysis, i.e., a model for each of three pheno-metrics for three regional landscapes.The residuals for the three regional landscapes-the difference between actual COV values of the NDVI-based pheno-metrics and their modeled outputs from the nine regression models-were mapped to identify and compare spatial patterns in the variation that was not explained by the models.
In this study area, although it is expected that the mountainous regional landscape will demonstrate the least variability over time (lowest COV), these values will not necessarily be indicative of higher resilience capabilities of this regional landscape but rather of high resistance capabilities (Table 1), which might be due to the species richness and altitudinal heterogeneity of the system, e.g., the species longevity capacities that tend to follow the elevation gradient of plant growth forms ranging from herbs to trees [58][59][60].The temperature-driven variability patterns of the steppe regional landscape are also expected to be indicative of the sensitivity to change in this regional landscape: higher COV values mean that the systems are susceptible to change due to low resistance capabilities, while lower COV values mean higher resistance capabilities under current and projected climate variability impacts.The healthy systems of this regional landscape are expected to have higher resistance and higher resilience capabilities (Table 1) to natural perturbation events, as the undisturbed steppe vegetation communities culminate with high biodiversity, productivity and ecosystem stability concurrently [61].Finally, because phenological responses of the desert regional landscape are precipitation-driven, they are expected to have higher variability (higher COV) values than those of the steppe regional landscape.Higher variability values are assumed to be indicative of lower resistance, and also might be indicative of either high (if the system maintains its heterogeneity) or low (degraded biodiversity) resilience capabilities of the desert regional landscape (Table 1).It should be noted that COV data might be very sensitive to areas with low vegetation cover, such as desert landscapes, and the interannual land surface phenological variations detected in these regions could partially be due to changes in soil surface reflectance, as well as due to the affects of calibration and atmospheric correction [56].In areas with high vegetation cover, such as in the northeastern region of Brazil, variation in the COV values helped to detect seasonal oscillations in vegetation dynamics when high COV values were detected during the dry seasons [52].

Yearly Models of Yearly Spatial Variation in Vegetation Dynamics as a Function of Climate, Soil, and Elevation Factors
In examining the functional relationship between the three NDVI-based metrics and the ten environmental drivers to explain the year-to-year variation in vegetation dynamics across the three regional landscapes of Central Asia (Appendices A1-A9), each environmental driver explained a comparable amount of variation in the pheno-metrics for each of the regional landscapes (Figures 2-4).Generally, the number of explanatory variables contributing to the spatial variation of the NDVI-based phenological metrics for all three regional landscapes across years was higher for the productivity metric than for the start and length of season timing metrics (Figure 5).Among the ten explanatory variables, summer temperature consistently contributed the least explanatory power to the spatial variation in NDVI-based phenological responses in all three regional landscapes and all 28 years (Figure 5).All the outputs displayed and discussed in this study contributed significantly (p ≤ 0.05) to the regression models.The following sections describe the modeling results for each of the three regional landscapes.

Steppe Regional Landscape
The temperature regimes were expected to be the most obvious and important factors affecting vegetation dynamics in the steppe regional landscape.Generally, climate variables were most frequently used and contributed the highest explanatory power in this region's models of the NDVI-based vegetation dynamics (Appendices A1-A3 and Figure 2).
The spatial variations in start of season were explained more consistently across years by climate drivers than by topography and soil drivers (Appendix A1).Generally, temperature regimes had consistently the most variation explained in start of season from 1981 to 2008 (Appendix A1 and Figure 2(A)).After holding the effect of all other variables constant, higher spring temperature and lower antecedent fall and winter temperature regimes were most highly related to earlier growing season start dates (Appendix A1).Among all environmental drivers, summer temperature regime showed the lowest number of years impacting the start of season: the summer temperature variable was present only during 10 years among 28 years of observations (Appendix A1).In all 28 years, at least six of the ten environmental variables contributed to the explanation of variation in start of season (Appendix A1).Explained variation in start of season by the selected environmental drivers varied from year to year and was highest in 1986, reaching 44% or an adjusted R 2 of 0.44 (will be referred as a percentage of variation explained hereafter), and was the lowest at 8% in 1985 and 2000 (Appendix A1).Separately, spring temperature explained the most spatial variation in start of season dates, having on average about 8% explained variation (Figure 2(A)).
Spatial variations in the length of season for the steppe regional landscape (Appendix A2) were explained well by elevation and spring precipitation patterns, with both metrics contributing to the best-fit models during 27 of the 28 years of observations (Appendix A2), and relating to shorter season length with increased rainfall and higher altitude after the effect of all other explanatory variables was held constant.Summer temperature had the fewest occurrences in the models, being present during only 13 of the 28 years (Appendix A2), but when assessed separately, it contributed on average the highest (9%) amount of explained spatial variation in the season length (Figure 2(B)).At least five of the ten environmental drivers related to the variation in length of season each year (Appendix A2).Across years the explained spatial variation in length of season did not exceed 42% (1986) and was lowest at 5% (1993) (Appendix A2).The spatial variations in the NDVI-based productivity metric for the steppe regional landscape were most consistently explained across years by soil carbon content, which was present in all 28 years (Appendix A3) and individually explained about 10% of variation on average (Figure 2(C)).Higher carbon content of soil led to higher values of the NDVI-based productivity metric (Appendix A3) when the effect of all other explanatory variables was held constant.Antecedent fall temperature was also consistently present in the models, occurring in 27 years (Appendix A3).Elevation had the lowest explanatory power (on average less than 1%) for the NDVI-based productivity metric, while antecedent fall, current winter, and spring temperature variables individually explained about 15% of spatial variation on average and up to 35% in some years (Figure 2(C)).NDVI-based productivity models incorporated at least seven of ten environmental variables across the years (Appendix A3).Spatial variations in the NDVI-based productivity metric for the steppe regional landscape explained by the model and selected environmental factors ranged between years, from 8% (2008) to 48% (1983) (Appendix A3).

Mountainous Regional Landscape
Across years, the total number of environmental variables that explained spatial variation in NDVI-based vegetation dynamics in the mountainous regional landscape (Appendices A4-A6) was less than that for the steppe regional landscape (Appendices A1-A3).Elevation explained the highest amount of spatial variation (Figure 3) and was one of the most frequent drivers explaining the variation in vegetation dynamics of the mountainous regional landscape across years (Appendices A4-A6).The explanatory power of the climate variables was lower than that of elevation, which is often viewed as a proxy for climate when assessing vegetation dynamics along elevation gradients [45] as temperature decreases and precipitation often increases at higher altitudes [62].However, in this region, whereas temperature variables seemed to be overshadowed by the elevation variable in explaining spatial variation, precipitation variables contributed quite consistently to the best-fit models during the 28 years.
The spatial variation in start of season in the mountainous regional landscape was most consistently explained by elevation and precipitation (Appendix A4).Higher elevation related to later season start dates when all the other explanatory variables were taken into consideration.Furthermore, when individual contribution of each environmental variable was examined, elevation contributed the most explanatory power, which on average was about 10% and reached up to 26% (1993) (Figure 3(A)).Although elevation values were not considered as a biophysical parameter, they served as a good proxy for temperature and precipitation constrains that vary with altitude, especially given the limitation of coarser spatial resolution of climate variables.The number of environmental variables explaining spatial variation in season start was as low as one (elevation, in 1984) and did not exceed nine variables in any given year (Appendix A4).For any of the 28 years the explained spatial variation in start of season by the set of selected explanatory drivers did not exceed 29% (1992) and was the lowest at 7% (1986; Appendix A4).The spatial variation in the length of season for the mountainous regional landscape was best explained by the soil carbon content and precipitation patterns across years (Appendix A5).Higher carbon content values related to shorter season length other variables were held constant (Appendix A5), which can be explained by the fact that soil carbon content varied with altitude and represented a proxy for soil moisture content in this regional landscape.Winter temperature contributed the least number of times to all models explaining spatial variation in length of season, being present during only six of 28 years (Appendix A5).Soil carbon content and elevation had the highest explanatory power at the individual level, each having on average about 15% of explained variation (Figure 3(B)).The number of environmental variables explaining spatial variation in season length ranged from three to nine variables across years (Appendix A5).The whole-model explanatory power of environmental drivers among all yearly season length models ranged from 5% (1990 and 1995) to 38% in 1994 ( Appendix A5).
The spatial variations in the NDVI-based productivity metric for the mountainous regional landscape were most consistently explained by elevation, which was present in all 28 years (Appendix A6) and individually explained on average about 25% of variation across years (Figure 3(C)).Values of the NDVI-based productivity metric decreased with higher elevation when all other explanatory variables were held constant (Appendix A6), which is probably due to cooler temperature with higher altitude.Winter precipitation had the lowest separate explanatory power (on average about 2%) for the NDVI-based productivity metric (Figure 3(C)).Summer temperature occurred the fewest times in the multilinear regression models for the NDVI-based productivity metric, present in 12 of the 28 yearly models (Appendix A6).The number of explanatory factors of spatial variation in NDVI-based productivity values ranged from four to ten variables per year (Appendix A6).The power of studied environmental drivers to explain spatial variation in the NDVI-based productivity metric ranged from 22% in 1991 to 44% in 1986 (Appendix A6).

Desert Regional Landscape
Because deserts are generally pulse-driven systems, the selected precipitation variables were expected to be the most noticeable factors affecting vegetation dynamics in the desert regional landscape.However, temperature regimes were the most frequently used variables in the models (Appendices A7-A9), and explained the greatest amount of spatial variation in vegetation dynamics in the desert regional landscape (Figure 4).The interpretation of the desert results is not intuitive.The start of the season is driven primarily by the higher temperatures in the fall, winter and spring, and secondarily by winter and spring precipitation.The timing of both temperature trends and precipitation events are likely impacting desert plant response through multiple interactions.Higher fall temperatures might cause more evapotranspiration and consequently a later growing season unless spring precipitation in combination with higher spring temperatures causes an early start of the growing season.Higher spring temperatures provide for earlier starts of the growing season if enough precipitation accumulated during winter and spring.
The spatial variation in start of season was very consistently explained by antecedent fall and current winter temperatures (Appendix A7).Higher values of these variables related to later season start dates if other environmental variables were held constant.Summer temperature showed the fewest occurrences among all environmental drivers of start of season, being present during only 13 out of 28 years of models (Appendix A7).For all the years examined, best-fit models explaining spatial variation in start of season included at least five of the ten environmental factors (Appendix A7).The explained variation in the start of season by the suite of selected explanatory drivers from year to year was highest in 1983 reaching 44% and was the lowest at 7% in 1984 (Appendix A7).Winter temperature had the highest individual explanatory power for the start of season metric and on average, about 12% of spatial variation was explained by this climate metric, reaching 30% in some years (2001-2002) (Figure 4(A)).The spatial variations in the length of season in the desert regional landscape were well explained; the length of season tended to be longer with higher antecedent fall temperatures and lower winter rainfall when all other variables were held constant (Appendix A8).Notably, taking into account that higher antecedent fall temperatures related to later season start and longer season length, it can be surmised that the end of the growing season is substantially later with higher antecedent fall temperatures.The summer temperature again was used the least in the length of season models, taking place only during nine of the 28 years (Appendix A8).The number of environmental drivers relating to spatial variation in length of season ranged from two to nine of the ten explanatory variables considered (Appendix A8).The explained spatial variation in the length of season by the suite of selected drivers in any given year did not exceed 31% (2005) and was the lowest at 2% (1993) (Appendix A8).Winter temperature regimes had the highest explanatory power at a separate level contribution assessment, having on average about 5% of explained variation (Figure 4(B)).
The spatial variation in the NDVI-based productivity metric for the desert regional landscape was most consistently explained by soil carbon content, which was present in 27 out of 28 years (Appendix A9) and separately explained on average about 13% of the variation (Figure 4(C)).When holding other explanatory variables constant, higher carbon content values related to higher NDVI-based productivity values (Appendix A9).Winter temperature was also consistently present in the models, occurring in 26 years (Appendix A9).Summer temperature occurred the least frequently in the models across years (Appendix A9) and had the lowest explanatory power (on average less than 2%) for the NDVI-based productivity metric (Figure 4(C)).The number of environmental variables that related to variation in NDVI-based productivity metric was at least six out of ten per year (Appendix A9).Across all years, the spatial variation in the NDVI-based productivity metric in the desert regional landscape explained by the suite of selected explanatory factors ranged from 73% (2003) to 22% (1990 and 1997) (Appendix A9).
An assessment of the environmental variables relating to spatial variation across years and across regional landscapes reveals that the steppe and desert regional landscapes had similar relationships among NDVI-based metrics and assessed environmental drivers (Figure 5(A,C)).Generally, temperature regimes, with the exception of summer temperature, were the most consistent variables relating to the spatial variation in the vegetation dynamics of both regional landscapes (Figure 5(A,C); Appendices A1-A3 and A7-A9)).As expected, the impact of precipitation variables on the spatial variation in pheno-metrics differed between these two regional landscapes [23], with spring and summer precipitation being the most important for steppe regional landscape and winter precipitation for the desert regional landscape (Figure 5).

Patterns of Variability in Vegetation Dynamics across Landscapes
Overall, the spatial distribution of the temporal variability in the three NDVI-based metrics (represented by COV) demonstrated distinct latitudinal and altitudinal patterns for all three regional landscapes (Figure 6).The pheno-metrics of the steppe and the mountainous regional landscapes showed similar magnitude in temporal variability (COV) patterns (Figure 6).This expected similarity is in contrast to the assessment of drivers of the spatial variation in the NDVI-based metrics, where the steppe and desert regional landscapes had analogous patterns (Appendices A1-A9 and Figures 2-5).
The pheno-metrics for the desert regional landscape showed to be more variable (higher COV) over time than those for the steppe and the mountainous regional landscapes (Figure 6).This high variability in vegetation dynamics (Figure 6) was expected and characteristic of the "pulse stability" [19] of this regional landscape.Its ecosystems are adapted to the particular intensity and frequency of the natural perturbation through lower resistance and higher resilience properties [19].
The variability (COV) of pheno-metrics in the steppe regional landscape was relatively low, revealing that systems of this regional landscape might have higher resistance.The variability in the mountainous regional landscape showed as expected low COV values, revealing that its systems demonstrate higher resistance.Below, all the implications of the short-term perturbations and long term stressors are discussed with regards to the climate and other natural variability, i.e., processes that are not primarily governed by humans.Each of the nine models developed to assess the drivers of variability in pheno-metrics had a different set of associated explanatory environmental variables (Table 3).Winter, spring, and summer precipitation and elevation were included in all nine models, with spring precipitation and elevation being the most consistent variables that contributed to these models (Table 3).Summer temperature was almost consistently absent from the models (Table 3).The mountainous regional landscape had less explanatory variables included into the models than the steppe and desert regional landscapes, which both had the same representative set of environmental variables in the models for each of the response variables (Table 3).The explanatory variables for the steppe regional landscape had the largest explanatory power for each of the pheno-metrics (Table 3).Among the NDVI metrics, the variability in the length of season was explained the best for all three regional landscapes (Table 3), in contrast to the NDVI-based productivity metric that was explained best among the year-to-year models (Appendices A3, A6 and A9).The spatial patterns in the residuals, obtained for each of the three regional landscapes, demonstrated that the models for the desert regional landscape produced the least amount of residuals in the COV values of the phenological response variables (mapped residuals in Figure 7) compared to the steppe and mountainous regional landscapes.Table 3. Summary results for the coefficient estimates for a series of stepwise multilinear regressions.The regression models were run to assess how the temporal variability (COV) in three NDVI-based metrics (season start, season length and productivity) is related to the suite of explanatory environmental variables for the three regional landscapes.Explanatory variables for each of the metric of each regional landscape were identified as most frequently occurring (at least 50% of time) variables during the 28 years of the year-to-year models.The bottom row includes the explained variability values (Adj.R 2 *100).Gray colored cells identify variables that were not frequently present in the year-to-year models and were intentionally excluded from these models.Dashed cells identify variables that were not contributing statistical significance to the overall explanatory power of the model.Note: T°-temperature; PPT-precipitation; C content-soil carbon content.Interannual variability in the season start metric for the steppe regional landscape was mostly related to temperature regimes, when all the other explanatory variables were held constant.Higher winter and spring temperatures related to increased variability (high COV) in season start while higher fall temperatures related to lower variability (low COV) in season start (Table 3).Although the explanatory power of all the variables contributing to the model reached 63% (Table 3), the spatial patterns in the residuals revealed that the model slightly overestimated (negative residuals) the variability in start of the growing season (Figure 7(A)).

Mountainous Regional Landscape
Winter temperature had the largest effect on the length of season's variability, with higher temperature relating to higher variability of the season length (Table 3).Overall the explanatory power of the environmental drivers contributing to the COV season length model was the highest among all nine COV-based models, reaching 70% (Table 3).Furthermore, the modeled explanatory variables performed relatively better for this model (lower absolute values of the residuals) than those for the COV season start model, slightly overestimating COV values in the southern part and underestimating it in the northern part of the steppe landscape (Figure 7(B)).
Variability (COV) in the NDVI based productivity metric was also predominantly affected by the temperature regimes.Lower fall temperatures and higher winter and spring temperatures related to higher variability in productivity responses across time (Table 3).The explanatory power of the COV productivity model did not exceed 25% (Table 3).As it can be inferred from the mapped residuals (Figure 7(C)), the model enabled a prediction of the spatial variability in the productivity metric of this regional landscape.
The variability (COV) values for the three pheno-metrics for the steppe regional landscape were generally low, suggesting that the regional landscape exhibits good resistance but might not be as resilient to change (Table 2).Given that the temperature regimes have the largest affect among other assessed environmental variables, the current shifting climate regimes might have direct implications for agricultural and natural land cover types of the steppe regional landscape.Although the cereal production in the northern and eastern part of the regional landscape can benefit from the longer growing season and warmer winters [26], natural land cover types might experience reduction of overall vegetation productivity and irregularity in the growing season, which altogether could lead to phenological asynchrony across trophic levels due to temporal mismatches between resource availability and consumer demand [63].
Modeling Interannual Variability of Vegetation Dynamics as a Function of Environmental Drivers for the Mountainous Regional Landscape In the mountainous regional landscape more rainfall during the winter and spring seasons led to lower COV season start values and the overall model explained 12% of the variation (Table 3).Higher fall temperatures related to more variability in season length dates, and contributed to the overall 17% of explanatory power of the COV season start model (Table 3).Warmer spring temperatures, lower soil carbon content, and higher precipitation amounts related to higher variability in COV productivity .The explanatory variables of this model explained only 6% of total variation (Table 3).In spite of the fact that all three models of the mountainous regional landscape (for COV season start , COV season length , and COV productivity metrics) had the lowest explanatory power across the three regional landscapes (bottom row in Table 3), these models enables prediction of the variability of all three response variables.This can be inferred from the mapped residuals of the mountainous regional landscape (Figure 7(A-C)), which were close to zero.
The COV values of the mountainous phenological responses were low, suggesting low susceptibility to change due to its higher resistance capabilities.Since precipitation has the largest effect on the three pheno-metrics representing this regional landscape, lower rainfall coupled with higher temperatures may have a negative feedback on ecosystems of the mountainous regional landscape.The changing climate regimes negatively affect mountain-restricted species causing shifts in treeline and migration of species towards the summits [64].Furthermore, because these climate regime shifts have caused changes in the extent of glaciers and much faster rates of snowmelt and glacier retreat [27,[65][66][67][68][69], they could lead to changes in the vertical zonation of mountainous flora and fauna and, thus to alteration of their phenological parameters.
Modeling Interannual Variability of Vegetation Dynamics as a Function of Environmental Drivers for the Desert Regional Landscape In the desert regional landscape more variability in the season start dates was related to higher spring temperatures, lower fall temperatures, and lower carbon content in the soil (Table 3).While the overall model explained only 25% of total variation in the COV, the mapped residuals of the explanatory variables were very close to zero (Figure 7(A)).
COV season length was higher with lower fall and winter temperatures and higher spring temperatures.The goodness of fit of the model was 48% (Table 3) where, as it can be inferred from the mapped residuals (Figure 7B), modeled variables had sufficient prediction power.
COV productivity values were higher with lower winter and summer temperatures and higher spring temperatures.The explanatory variables of the model explained about one third of the total variation (Table 3) and with the residuals being close to zero (Figure 7(B)), they seemed to predict quite well the variability of all three of the response variables.
The observed high COV values for the desert pheno-metrics suggested that this regional landscape is very susceptible to change, particularly with lower resistance (Table 2).It is expected that desert landscapes have higher resilience, components of which include ecosystem elasticity (rate of restoration to a stable state after perturbation) and amplitude of stability (deformation extent from which an ecosystem will return to its initial state) [70].Although the systems of the desert regional landscape are assumed to be able to recover after short-term perturbation and long-term stress events, the long-term impacts of projected climate change could likely degrade these systems that are susceptible to change [26,71,72].

Conclusions
Through the application and analyses of digital biophysical time series data, this research provides new and valuable understanding of the impact of environmental drivers on the spatial, annual and interannual variation of vegetation dynamics in three regional landscapes of Central Asia.Additionally, this study links measures and scales of observed interannual variability in vegetation dynamics in steppe, mountain and desert landscapes to their sensitivity to change patterns, specifically to their resistance and resilience capacities under existing and projected environmental regimes of change and variability in these regional landscapes.
Although this study has provided important insights into regional landscape-scale vegetation dynamics, the results can likely be further enhanced by using more consistent and finer resolution spatial and temporal datasets.This matter is especially relevant to the climate datasets.The climate datasets are often compilations of records of climate data gathered from a variety of sources [73] and uncertainty might be introduced in climate datasets from the interpolation between weather station locations, elevation bias in the weather stations, elevation variation within grid cells, and through data partitioning and cross-validation [73].The uncertainty in climate data is generally the highest in mountainous and in poorly sampled areas [73].Additionally, the estimations of precipitation patterns might have higher uncertainty as they might be spatially less accurate than the temperature data due to nature of precipitation patterns and resolution of the rainfall data.Therefore, potential future research might include replicating the applied approach at a finer spatial scale with finer-resolution datasets.The use of minimum and maximum temperature regimes is worth considering as well [74,75].
Specifically, the research identifies potential drivers of spatio-temporal patterns in landscape-scale vegetation dynamics (phenology and productivity) in three regional landscapes (desert, steppe, and mountainous), which then were used to map the temporal variability of these dynamics.The results of this study suggest that a unique suite of specific environmental drivers governs the vegetation response of steppe, mountainous, and desert regional landscapes.Central Asia's landscape-scale vegetation dynamics showed strong relationships to environmental factors in across latitudinal and elevation gradients.The spatial and year-to-year variation in vegetation dynamics of the steppe regional landscape were demonstrated to be mostly affected by temperature.The temporal variability (COV values) of the response metrics in this region was also mostly affected by temperature, although temperature variables did not explain the observed spatial patterns of the COV values (mapped residuals) completely.Finally, relatively low temporal variability of vegetation dynamics suggested higher resistance to change in this regional landscape.
In the mountainous regional landscape the spatial variation in vegetation responses over 28 years were related primarily to climate variables, especially precipitation observed along the elevation gradients.The models that assessed temporal variability of the vegetation dynamics also suggested altitudinal clines of climate patterns as drivers.Although drivers of the temporal variability in vegetation dynamics had lower explanatory power (adj.R 2 values) than the steppe and desert environmental drivers, they still explained the spatial patterns of variability relatively well.Lower values of the COV of vegetation dynamics can be assumed to relate to higher resistance and lower susceptibility to short-term perturbations in this regional landscape.
Vegetation dynamics in the desert regional landscape were demonstrated to have the greatest variability (high COV) over time that seemed to be affected by the combination of temperature regimes, winter precipitation patterns, and soil carbon content.The high COV values of the desert regional landscape suggest its higher susceptibility to change and low resistance capabilities.
Given that climate precipitation and temperature variables showed the largest impact among other assessed environmental variables on vegetation dynamics across all three regional landscapes, projected climate changes could severely impact water use and agricultural production.Land use policy and decision makers in each of the regions will need to develop economically and environmentally effective climate change adaptation and mitigation plans to sustain natural resources and agriculture in Central Asia.Healthy ecosystem functions need to be maintained through processes that could confer adaptive capacities such as rapid recovery, high growth rates, rapid succession of ecosystem development stages, and/or flexible and opportunistic timing of vegetation growth responses.Because components of landscape susceptibility to change include resistance and resilience [70], the systems of the three regional landscape might be able to adapt to projected climate changes [25] if the future rates of change turn out to be less rapid and more gradual than is projected.As this is likely not the case, climate change implications in all three regional landscapes could involve both mitigation (e.g., increased water use efficiency, water storage) and adaptation (e.g., reduced water use and crop changes) efforts across the extent of these landscapes.
Appendix A Appendix A1.Steppe start of season summary table for series of stepwise multilinear regression models, run for each year from 1981 to 2008, where the start of season metric is a response variable to a suite of ten explanatory environmental variables (described in Table 1) and listed here (top row) in a descending order of total number of the year-to-year occurrences (bottom row) in the 1981 to 2008 models.The table presents coefficient estimates (β) for each of the environmental variables that contributed significantly (p ≤ 0.05) to the explanatory power (last column) of the regression models and were included in the total number of environmental X−variables (second column) present per year (first column).Note: T°-temperature; PPT-precipitation; Ant.Fall-antecedent fall.

Figure 2 .
Figure 2. Yearly adjusted R 2 values obtained for each of the response variables and their corresponding ten explanatory variables for the steppe regional landscape.(A) Start of the growing season.(B) Length of the growing season.(C) The NDVI-based productivity metric.

Figure 3 .
Figure 3. Yearly adjusted R 2 values obtained for each response variable and their corresponding ten explanatory variables for the mountainous regional landscape.(A) Start of the growing season.(B) Length of the growing season.(C) The NDVI-based productivity metric.

Figure 4 .
Figure 4. Yearly adjusted R 2 values for each of the response variables and their corresponding ten explanatory variables for the desert regional landscape.(A) Start of the growing season.(B) Length of the growing season.(C) The NDVI-based productivity metric.

Figure 5 .
Figure 5. Frequency of appearances (in percentage) of each of the ten environmental variables during the 28 years of modeling.(A) Steppe regional landscape.(B) Mountainous regional landscape.(C) Desert regional landscape.

Figure 6 .
Figure 6.Distribution of the coefficient of variation (COV) values for: (A) Start of the growing season.(B) Length of the growing season.(C) The NDVI-based productivity metric.

C
Drivers of Variability in Vegetation Dynamics for Each Regional Landscape

Figure 7 .
Figure 7. Distribution of the residual values for the modeled COV assessment for: (A) Start of the growing season.(B) Length of the growing season.(C) The NDVI-based productivity metric.

1 C
Modeling Interannual Variability of Vegetation Dynamics as a Function of Environmental Drivers for the Steppe Regional Landscape

Table 2 .
Response and both stable (s) and dynamic (d) explanatory variables to be used in the regression and correlation models.The unused explanatory variables are listed as well.

Steppe length of season summary
table for series of stepwise multilinear regression models, run for each year from 1981 to 2008, where the length of season metric is a response variable to a suite of ten explanatory environmental variables (described in Table1) and listed here (top row) in a descending order of total number of the year-to-year occurrences (bottom row) in the 1981 to 2008 models.The table presents coefficient estimates (β) for each of the environmental variables that contributed significantly (p ≤ 0.05) to the explanatory power (last column) of the regression models and were included in the total number of environmental X−variables (second column) present per year (first column).Note: T°-temperature; PPT-precipitation; Ant.Fall-antecedent fall.