Article Phenological Characterization of Desert Sky Island Vegetation Communities with Remotely Sensed and Climate Time Series Data

Climate change and variability are expected to impact the synchronicity and interactions between the Sonoran Desert and the forested sky islands which represent steep biological and environmental gradients. The main objectives were to examine how well satellite greenness time series data and derived phenological metrics (e.g., season start, peak greenness) can characterize specific vegetation communities across an elevation gradient, and to examine the interactions between climate and phenological metrics for each vegetation community. We found that representative vegetation types (11), varying between desert scrub, mesquite, grassland, mixed oak, juniper and pine, often had unique seasonal and interannual phenological trajectories and spatial patterns. Satellite derived land surface phenometrics (11) for each of the vegetation communities along the cline showed numerous distinct significant relationships in response to temperature (4) and precipitation (7) metrics. Satellite-derived sky island vegetation phenology can help assess and monitor vegetation dynamics and provide unique indicators of climate variability and patterns of change.


Introduction
Climate patterns from intra-seasonal to decadal and century scales directly influence the timing, magnitude (productivity), and spatial patterns of vegetation growth cycles, or phenology [1,2].Vegetation phenology can be detected by satellites and other remote sensing systems due to the unique seasonal and spectral reflectance and transmittance characteristics of canopy, plants and leaves [1,3].Compared to other land surface components (e.g., soils), red and NIR reflectance values for green vegetation are low and high respectively.The Normalized Difference Vegetation Index (NDVI = (ρ NIR -ρ red )/(ρ NIR + ρ red )) [4] is one of the most commonly used spectral vegetation indices for agricultural and natural resource monitoring applications and has been related to many biophysical properties, such as fractional green cover, chlorophyll density and leaf-area index (LAI), the fraction of absorbed photosynthetically active radiation (fAPAR), biomass and net primary production [5][6][7][8].Remotely sensed vegetation phenological data have been used in global climate change studies, showing trends and responses such as earlier start of the growing season, later end of season, and higher seasonal productivity [9,10].In addition, remotely sensed vegetation phenology information is used to determine land cover classes [11,12] and vegetation dynamics in response to climate [13] and disturbances and stressors like fire [14] and drought [15].Metrics created from spectral vegetation indices, such as start of season, peak index or greenness value, and seasonal integrated greenness are used in many models and studies to describe phenology and to represent vegetation interactions with climate-based [5,10,16] and anthropogenic [17,18] factors.
Most remote sensing phenology studies use satellite-based vegetation index (VI) time series as a source to derive phenological metrics (phenometrics) to identify critical measures (e.g., timing of start, peak, length of the growing season, amounts of greenness) in the annual seasonal vegetation index signatures that match phenological phases of vegetation development.Many algorithms have been designed to retrieve phenometrics [19].Jönsson and Eklundh [20,21] designed a phenometrics retrieval approach that is based on both threshold and curve derivative methods.This approach reduces residual atmospheric noise in the time series data, and provides a user-friendly output in either a point-based or image format.

Nature of the Study Area
Sky islands are biogeographically isolated mountains, which are connected by lowland valleys and terrain that act as both barriers and bridges to dispersion and colonization by new species from one mountain or valley to another mountain or valley.Many sky islands throughout the world have resulted from isolation caused by erosive and tectonic forces over time and the interplay with varied bedrock formations [22].Sky islands contain a large number of diverse landcover types [23] and contain a large part of the world's biodiversity [24].Drought stress and wildfire disturbance [25,26] threaten the rare and endemic species found there as well as the ecosystem functions of sky islands [27,28].Sky islands also fulfill an important role as safe havens and food resources along bird migratory routes [29,30].Sky islands have been shown to be regionally important as ecological regulators of water and CO 2 budgets [31].
This research focused on one of the sky islands, which is part of the Madrean Archipelago, a mid-latitude inland sky island complex between the Rocky Mountains and the Sierra Madre Occidental, in the midst of the Chihuahuan and Sonoran Deserts (Figure 1).Like all sky islands, those of the Madrean Archipelago are characterized by steep ecological gradients [32][33][34].Southeastern Arizona's sky islands have a relatively small spatial footprint compared to the surrounding desert, with steep gradients and highly variable distribution of vegetation communities.In contrast, remotely sensed phenology data often have coarse spatial resolutions with pixels of 1 km to 1 degree in size [35].However, the 250 m resolution NDVI time series data from the Moderate Resolution Imaging Spectroradiometer (MODIS) [36] have been used successfully for monitoring natural resources [37] and modeling vegetation response after wildfires using land surface phenological metrics [14,38,39].Coarse global change and land classification models often exclude entire arid biomes due to the weak phenological signatures and low productivity of these areas [40].Coarsely resolved phenological signatures may overlook small-scale patterns that indicate important ecosystem processes [41].Thus, in research of remotely sensed vegetation phenology there is a need to further explore climate-phenology interactions to improve our understanding of the impact of climate [42] on arid lands and finer scale landscape features like sky islands.

Study Objectives
The objectives of this study were to (1) characterize and examine the variation in sky island vegetation growth and phenology patterns using satellite time series of NDVI data, (2) derive and examine vegetation type-specific phenological metrics, and (3) examine the relationship between phenometrics and climate variables along an elevation cline.We hypothesized that the timing and the magnitude of the growth cycles of vegetation communities vary according to precipitation and temperature patterns, considering that these vegetation types are mostly limited by water and insolation regimes [43].Vegetation communities might have unique responses to these drivers that might be visible via phenometrics.Specifically, we extracted and analyzed phenological and climatic metrics for the Santa Rita Mountains, a sky island in southeastern Arizona, and (1) examined vegetation phenological characteristics for different vegetation types across a steep elevation gradient; (2) visualized spatial and temporal patterns of phenometrics; (3) examined the relationships between climate metrics (e.g., annual precipitation) and phenometrics; and (4) explored which of these phenoclimatic relationships were distinctive for specific vegetation communities.

Study Area
The study site extent was a 4,400 km 2 area south of Tucson, Arizona, which encompasses the sky island of the Santa Rita Mountains (Figure 1).The high NDVI values show the sky islands to be "refuges" in the desert "seas" (Figure 1), with a diversity of land cover types (Figure 2 [44]).The study area includes vegetation communities ranging from lowland palo verde-mixed cactus desert and mesquite grasslands to oak woodland and coniferous forest along an elevation gradient of about 2,000 m (Figure 2c).We assessed the vegetation dynamics in this area over the time span of mid-2000 through mid-2007, which allowed us to examine the growing seasons of 2001-2006.

Research Design: Sampling Land Cover, NDVI and Climate Time Series Data
Both image-based and point-based NDVI time series were examined with respect to their growth and phenology patterns.However, only areas with relatively homogeneous landcover types were selected to examine the variation in phenology and their response to climate variables.To sample the study area, several criteria were taken into account to ascertain representation of the major landcover types and avoid mixed or heterogeneous landcover patterns.The NDVI time series data were used to derive land surface phenometrics-the response variables in our modeling approach.Precipitation and temperature time series were the explanatory variables and were available at 4 km resolution.The land cover data were available at 30 m resolution and were used to stratify and select our samples that represented homogeneous areas within the same land cover type.Since the spatial resolution of the different data sets ranged from 30 m to 4,000 m, analysis was performed at the intermediate resolution of 250 m.

SWReGAP Landcover and DEM
To select sampling sites from which to extract time series data, we first derived vegetation community classes from the Southwestern Regional Gap Analysis Project (SWReGAP) dataset [44] (Figure 2a).The SWReGAP data were re-sampled to 250 m to match the MODIS NDVI product.The re-sampled SWReGAP product represented the majority vegetation community per MODIS pixel.We then filtered the data by analyzing the amount of contiguous 30 m SWReGAP pixels within a 250 m MODIS pixel, with at least 80% homogeneity in vegetation communities (by comparison with the original dataset).The resulting dataset included vegetation communities having at least three sites where four or more 80% homogeneous pixels were available for analysis (Figure 2b).Areas that were either non-natural or urban, and areas with high heterogeneity were not used in the analyses.The final dataset included eleven vegetation types, the full extent of which represented over 98% of the areal extent of the study site and corresponding vegetation clines, which are characteristic of many of the Sky Islands of the Madrean Archipelago.We sampled 41 sites (four pixels-25 ha-per site) throughout the study area: three to five sites per vegetation type, including six weather stations in six different vegetation types (Figure 2c).

Figure 2. (a)
The vegetation communities based on the SWReGAP vegetation cover map are dominated by the 11 (b) communities used in our analysis (see also Table 1).The vegetation community data set was re-sampled to 250 m and filtered for homogeneous sites.White areas/pixels denote either developed urban and agricultural areas, or areas without enough contiguous homogeneous pixels for the analyses.(c) The locations of homogeneous sample points and weather station are distributed over the elevation gradient.

MODIS NDVI Time Series
NDVI data at 250 m resolution and16-day composite intervals, acquired by MODIS on the Terra platform (MOD13Q1), were used for the NDVI time series and phenological analysis.Both the MODIS NDVI data and associated quality assurance (QA) data [36,45] were used in this research.This resulted in 23 composite NDVI and QA images (periods) per year, providing a total time series of about 161 NDVI/QA images.The MODIS NDVI images were reprojected into a Lambert Azimuth equal area projection suitable for analysis.For each of the selected homogeneous sites we extracted MODIS NDVI time series data from the 12th period of 2000 through the 11th period of 2007.Seasonal phenological metrics were then extracted for 2001 through 2006.The mean of NDVI time series data were extracted for all areas of interest excluding pixels that were affected by clouds or snow as indicated by the QA data.These data provided the basis for examining the trends, seasonality and anomalies in NDVI and the associated phenology metrics, and their relationships to climate variables, for the different land cover types.

Precipitation and Temperature
Climate data (Table 1) were extracted from Oregon State University's Parameter Regression on Independent Slopes Model (PRISM) dataset, at 4 km spatial resolution and monthly temporal resolution [46].The beginning of vegetation green-up or start of the growing season, marked by Start t is defined as the point in time for which the NDVI value (Start NDVI ) has increased by 20% of the distance between the left minimum NDVI value and the Peak NDVI value (Figure 3).When vegetation senesces and less vegetation cover is observed on the land surface at the end of the growing season, End t and End NDVI can be estimated similarly.The 20% threshold accounts for noise in the NDVI data and prevents the identification of false starts or ends of the growing seasons.The peak of a growing season is estimated by the timing of Peak t for which the NDVI value of reaches the maximum NDVI value (Peak NDVI ).The amplitude of the season is obtained as the difference between the Peak NDVI value and the average of the Start NDVI has End NDVI values.The length of season is computed as the difference between the start and end of the growing season, (End t -Start t ).Two separate integrals over the growing season are computed to estimate the production of the seasonally dominant vegetation type.The first integral (Small Integral-SI), given by the area of the region between the fitted function and the average of Start NDVI and End NDVI values, called the Base NDVI , represents the seasonally active vegetation, which may be large for herbaceous vegetation cover and small for evergreen vegetation cover.The second integral (Large Integral-LI), given by the area between the fitted function and the zero NDVI value bounded by the Start t and End t , represents the total vegetation stand and is a proxy for vegetation production.The left derivative or rate of increase in NDVI during the beginning of the growing season is related to the rate of green-up and plant growth and can be estimated by calculating the ratio between the amplitude and the time difference between the start and midpoint of the growing season.The right derivative is related to the rate of senescence, but is also governed by snow cover, which can obscure green vegetation from being observed by the satellite sensor.Different vegetation cover types are likely to have different rates of green-up and senescence [48].

Data Analysis
The PRISM precipitation and temperature data were compared to data sets from six ground-based weather stations within the study sites to examine the consistency of the modeled climate data.These data were used in subsequent analyses of the relationships between phenometrics and climate variables for dominant vegetation types along the elevation cline.These stations provided monthly precipitation data from 1989 through 2006 for a total of 425 measurements.Three of the weather stations also provided measurements from 2000 through 2006, for minimum temperature (216 measurements) and maximum temperature (219 measurements).Weather station measurements were well correlated (α = 0.05) in all three datasets (R 2 = 0.696, 0.928, and 0.904 for precipitation, maximum temperature, and minimum temperature, respectively).
A qualitative evaluation of the phenological characteristics for the sky island vegetation types was conducted by examining the seasonal NDVI trajectories.The spatial and temporal variation in NDVI among and within vegetation types were examined by comparing spatial patterns between wet and dry year exemplars.The difference between the wet and dry year NDVI images were used to evaluate the inter-annual spatial variability of sky island vegetation growth patterns.The inter-annual variation between the yearly phenometrics (e.g., Start NDVI Start t and LI) was quantified by the coefficient of variation (COV = standard deviation/mean).Images of COV phenometrics were computed based on each phenometric for each year between 2001 and 2006.Zero-or non-values were excluded from COV calculations.The relationships between the annually derived phenometrics and seasonal and annual climate metrics were explored using the following statistical techniques.To assess the variability of phenometrics with respect to climate and vegetation types we used correlation, MANOVA and stepwise regression analysis techniques.We used a Bonferroni adjustment to account for multiple statistical comparisons.Data were left out in instances where seasons were not detected by TIMESAT or where climate data were not available.JMP statistical software was the main tool used in our analysis (SAS Institute, 2006).

NDVI and Climate Time Series Patterns among Vegetation Types
NDVI images after the spring rainfall (May) and monsoon rainfall (September) seasons show the inter-annual, seasonal and spatial vegetation dynamics for a relatively "wet" (2001) and "dry" (2002) year (Figure 4).Year 2001 showed a somewhat higher NDVI compared to 2002 for both May and September.However, the spatial patterns observed in the difference from yearly average (2000-2006) NDVI values for May and September of 2002 showed that overall greenness for May 2002 was much lower than average, particularly for the higher elevation areas.Generally, for September 2002 the difference from average NDVI values increased for lower elevation lands and decreased for mid-elevation lands.The riparian corridor formed by the Santa Cruz River showed a variable response to the 2002 spring drought, which might relate to a smaller effect of rainfall events.
Averaged time series data for each of the sky island vegetation types showed that the phenological signatures were generally not coinciding with each other in terms of the shape and magnitude of the NDVI response (Figure 5).Differences in the magnitudes of the high and low NDVI values as well as in the timing and shape of seasons can be seen among (Figure 5) and within (Figure 6) vegetation communities.Time series data also showed a high degree of variability across seasons (Figure 5).The seasonal NDVI patterns were similar for some communities and increasingly different for vegetation types that are further apart in terms of elevation range.The NDVI values generally increased with elevation.There was also strong seasonal and inter-annual NDVI variability within vegetation types.The cross-site and inter-annual heterogeneity of the NDVI time series data highlights both the steep gradient of the Santa Rita Mountains and the variation to be expected in any vegetation community [41].When NDVI time series were averaged for each vegetation type and compared to one another, differences in phenology could be seen that often suggested different biological or ecological processes associated with each of the vegetation types (Figure 5).For example, the phenological pattern for the pine-oak woodland showed overall higher base and peak NDVI values, less obvious seasonality and longer seasons; this could be expected in an evergreen vegetation community where seasons are often only punctuated by snow cover [1,45].The bimodal growing season is not visible at higher elevations as the woody vegetation is rooted much deeper and responds more to temperature related cues.Figure 6a shows an example of the variability in NDVI for sampling sites that were all classified as a mesquite community.The spring and monsoon rainfall periods resulted sometimes in a bimodal growing season, generally corresponding to weak increases in NDVI values in the spring and strong increases in NDVI values during the monsoon season.Increasing minimum night time temperatures and decreasing precipitation trends (Figure 6b) were observed to be statistically insignificant (p > 0.01) but have been predicted in the literature [27].These climate trends, warmer summer temperatures in particular, have resulted in species flowering range changes in the Santa Catalina Mountain range, a sky island just north of the Santa Rita Mountains [50].

Spatio-Temporal Patterns of Phenometrics
Assessment of image-based spatial and temporal patterns of phenometrics revealed both spatial and temporal patterns (Figures 7-9).The start of the growing season was comparatively spatially asynchronous in 2001 (Figure 7a) and synchronous in 2002 (Figure 7b).These patterns in synchronicities could be related to a bimodal growing season in 2001and a uni-modal growing season (skipped the spring due to the lack of precipitation) in 2002 for most vegetation communities.Annual integration of the NDVI values show that vegetation production was higher in 2001 (Figure 7c) than in 2002 (Figure 7d), and suggests the impact of inter-annual climate variability.
Interestingly, we observe that the seasonally integrated or small integral (Figure 8) is lower for the higher elevation vegetation cover types than the middle elevation cover types, suggesting more seasonally dynamic responses at these intermediate elevation levels.This might be attributed to the broadleaf deciduous woody plants (e.g., oak) at these elevations that senesce in the fall and winter, while higher-elevation coniferous vegetation communities keep their green plant matter.The observed spatial patterns in the COV for the Start of the season and Small Integral (Figure 9) suggest that the inter-annual variability in the start of season (Start t ) and seasonal small integral (SI) responses are somewhat different across vegetation types and among phenometrics.

Start t SI (b) (a)
The effect of a forest fire in 2005 is also clearly visible at the top of the mountain range; this area was excluded from further analysis.Distinguishing phenological characteristics of vegetation types, including both their mean values and their variability, have been exploited in some land cover classifications [51,52], albeit at coarser biological and spatial scales than this study.Patterns shown in our phenometrics images suggest that classification might be a possibility at this spatial scale as well.

Phenometrics and Climate across Vegetation Communities
Phenometrics were highly correlated to many climate metrics in our analysis across vegetation types (Figure 10).Significant positive correlations to precipitation and negative correlations to temperature metrics were seen for many of the phenometrics describing the magnitude of seasons, such as Base NDVI , amplitude, and large integral.Similar patterns have also been seen in global studies of NDVI as related to precipitation and temperature [53].Linderman et al. [54] found that phenometrics varied across the African continent due to climate more than due to land use change.Correlations of timing metrics (start, peak, and end of season) with climate were less strong than correlations of many of the magnitude metrics, suggesting that there might be different causes than temperature and precipitation for variability in these phenometrics (e.g., soil and vegetation type, growing degree days) [2,55].Across vegetation types an earlier start of season was correlated with lower temperatures and higher precipitation.Start of season showed negative correlations for precipitation and positive correlations for temperature (i.e., an earlier season start correlates with more rain and lower temperatures), while peak and end of season were shown to be later in the season with higher precipitation and lower temperatures.As the Southwest was in a drought for most of our study period (2002 through 2005 and ongoing), accompanied by anomalously warm temperatures at least in the worst drought year (2002) [16], relationships showing dependence on water and cooler temperatures are not surprising [56].Increased precipitation could delay plant senescence and therefore the peak and end of the growing season.Lower temperatures might allow plants to continue to photosynthesize due to lessened vapor pressure deficit or other factors affecting drought stress [16,40].Also, these patterns might highlight the dominance of the lower-elevation desert vegetation types in our study area, communities which experience less low-temperature-limitations, except in the winter.
Stepwise regression for each metric likewise showed salient relationships between metrics and climate (Table 2).Each phenometric showed moderate to highly significant relationships with climate metrics, and the combination of relationships was unique per metric.The amount of variation in phenometrics that was explained by climate metrics varied from less than seven percent to over 74% (adjusted R 2 ; Table 2).Base NDVI , peak NDVI and large integral were best explained by climate metrics, while peak and end dates and rates of green-up and senescence were not well explained by climate metrics across all sample sites.Annual precipitation explained variation in almost all phenometrics.

Phenometrics and Climate by Vegetation Type
Analysis of phenometrics by vegetation type (11) showed correlations between both timing-and magnitude-based phenometrics and climate metrics.Correlations of a high and low elevation vegetation community are shown in Table 3 and Table 4 respectively, while an overview of all significant relationships for all 11 vegetation types are provided in Appendix A. After a Bonferroni adjustment for multiple comparisons, 121 tests (11 phenometrics × 11 climate variables) and a mean correlation of 0.62 [57], lowering the alpha from 0.05 to 0.0081 for each test, 128 pheno-climate relationships were significantly correlated (p ≤ 0.0081).Of these 128 relationships between phenometrics and climate metrics, 88 were related to precipitation metrics and 40 to temperature metrics.Another 198 pheno-climate relationships were significantly correlated at 0.0081 < p ≤ 0.05 (without Bonferroni adjustments; data not shown).Of these 198 relationships between phenometrics and climate metrics, 118 were related to precipitation metrics and 70 to temperature metrics.The pheno-climate relationships are dominated by precipitation metrics (66%), although temperature metrics also plays a considerable role (34%) among these.Negative correlations were found for 124 out of the 326 total pheno-climate relationships which were evenly distributed between temperature and precipitation metrics, but about 65% of these were related to timing based phenometrics.An overview of the number of significant pheno-climate relationships across all vegetation types by phenometric and climate metric is provided in Appendix B. This overview suggests that some of the phenometrics, like the end of the season and rate of green-up are not as often correlated to climate metrics as some of the other phenometrics.It also suggests that precipitation-based metrics are relatively more important than temperature-based metrics for this gradient, while the winter precipitation (4 mo.) has the fewest number of significant relationships.
As in our analysis across vegetation types, magnitude-based metrics were often positively correlated with precipitation and negatively correlated with temperature.Significant relationships were found in all of the vegetation communities.Correlation matrices suggested unique relationships between phenometrics and climate metrics for each community along the elevation gradient (data shown in Appendix A).For example, in the pine-oak woodland vegetation community, most phenometrics were correlated with precipitation metrics and not with any of the temperature metrics (Table 3).The base, peak and integrated NDVI metrics were positively correlated with several precipitation metrics (i.e., winter, pre-monsoon and annual precipitation), while most of the timing based phenometrics were negatively correlated with several of the precipitation metrics.For example, an earlier start of the growing season was a function of both pre-and post-monsoon precipitation, indicating that soil moisture accumulation during fall and winter will result in an earlier start of the next spring or summer growing season.The rate of senescence is slower when we have winter precipitation.In contrast, in the grassland community (Table 4) a later start of season was correlated with extreme maximum and minimum temperatures, while an earlier start of the season was correlated with more winter and pre-monsoon precipitation.The peak of season is earlier with more pre-monsoon, annual and winter precipitation and lower maximum temperatures.The end of season seems to be later with higher average minimum and maximum temperatures, suggesting that these temperatures extend the growing season.It is less evident why less winter precipitation would result in a later end of the growing season.More complex temperature, precipitation and environmental feedback mechanisms and interactions would need to be explored further.Most magnitude based metrics like SI, LI, NDVI peak , and NDVI base increase with precipitation metrics (e.g., summer, annual, and winter precipitation), and decrease with extreme low or high temperatures.
While some of the magnitude-based phenometrics correlated similarly to climate metrics for many vegetation types, correlations of timing-based phenometrics seemed more inconsistent between vegetation types.It could be that the timing of plant growth cycles is more variable than magnitude on a landscape scale.DeFries et al. [51] found that vegetation types were distinguished from one another using unique metrics for each separation; our study suggests the same.Heumann et al. [58] found that different TIMESAT-derived phenometrics responded to climate factors in different regions; similar patterns are illustrated in this research by the ability to distinguish vegetation communities across this steep climatic gradient that represents many ecoregions.As individual plants respond directly to soil moisture patterns rather than rain or snow per se [59], highly variable timing might be related to soil type and texture, as well as plant functional type [60], topography, or antecedent precipitation [59,61].And finally, with the studied landscape also experiencing drought it might be expected to display vegetation phenology that is highly sensitive to precipitation pulses and therefore potentially quite spatially heterogeneous [56,59,62].In addition to limitations in spatial resolution, the temporal resolution of the standard 16-day NDVI data products could potentially also preclude a more precise capture of vegetation community phenology; certain growth phases can be overlooked within a 16-day period [1,63], although most timing-based metrics tend to be within 3 days when based on 16-day composited MODIS time series data [64].Use of more finely resolved and cloud-free time series data might provide more precise phenology and could positively or negatively affect the relationships between phenology and other environmental factors.Finally, the algorithms or the settings used in TIMESAT could be optimized for each vegetation type, possibly resulting in slightly different metrics of the phenology for each of these communities.However, tracking the optimal setting for each pixel might not be practical and could lead to less spatially and temporally consistent phenometrics.

Conclusions
Use of phenometrics to describe vegetation response to climate trends and variability across a sky island, with the finer spatial resolution of the Moderate Resolution Imaging Spectroradiometer (MODIS) Normalized Difference Vegetation Index (NDVI) at 250 m and 16-day composite intervals, allows for a more detailed look at the vegetation communities in and around steep biological gradients than current phenology studies typically employ [13,35,65].The temporal and spatial heterogeneity in land surface phenological patterns along this climate and elevation gradient were well characterized and depicted by phenometric images and examples of COV images for the start of season and small NDVI integral.Examples of spatial patterns in land surface phenometrics showed that the start of the growing season and annual productivity for the 11 vegetation communities are highly variable between years.Example COV images provided an important spatial representation of the distinct interannual variability among phenometrics.
This phenoclimatological study provides the groundwork and support for our hypotheses.Remotely sensed phenological signatures and phenometrics differed across the elevation gradient and were unique for each vegetation type.Most vegetation types (11) and associated phenometrics (11) showed distinct relationships and significant correlations (128 relationships, p ≤ 0.0081, Appendix A; 326 relationships, p ≤ 0.05, Appendix B) for many of the seven precipitation and four temperature metrics.The climate metrics were less successful in explaining variability in each of the phenometrics across all sample sites at once (representing all 11 vegetation community types).This suggests that some phenological similarities exist across the elevation gradient, but that most of the variability among vegetation type communities could not be fully captured by the generalized pheno-climate models created in this study (Table 2).
Low-elevation desert sites (e.g., creosote, palo verde and mesquite communities) had seasonal shapes that included sharp seasonal peaks in NDVI and slight to strong bimodality.These features might point to the potentially high degree of sensitivity of these communities to precipitation events [1,49].Additionally, inter-annual variation might be related to temperature; low temperatures might postpone the start of the growing season in the spring, especially at higher elevations, while high temperatures might result in high evapotranspiration rates, possibly postponing the onset of the growing season or reducing the vigor and possibly the length of the growing season, especially in the lower elevations where shallow rooted herbaceous plants and shrubs are more water limited than plants at higher elevations.
Spatial and temporal resolution of the MODIS and PRISM data might have hindered more precise relationships, but our study demonstrates the value of these datasets in exploring landscape-scale phenology across steep elevation and linked biological gradients.The number of observations could be augmented and the results made more ecologically robust by including more samples from other Sky Islands.Environmental factors beside climate, such as topography (slope and aspect) and soil physical and chemical properties (water permeability, surface rock fragments, soil depth, carbon, and nutrient availability) might contribute the phenological changes as well and could be a focus of future research.Impacts of pulse driven and lagged land surface phenological response to climate variables can also be further explored, especially if we examine sky island resilience to changes.
Use of phenometrics in this manner could be helpful in updating land cover classifications over seasons [66].Monitoring ecotones and transition zones of sky islands could be practical for management and conservation purposes, and further development of the presented methods could result in tools for such monitoring.Additionally, sky islands could be used as indicators and barometers of stress [28].Sky islands provide the ability to do this, as they represent a latitudinal gradient as well as many of the US Southwest's dominant vegetation types within a relatively small area.Continuing research into the spatio-temporal phenological patterns of sky islands could assist in the conservation of these areas as biodiversity hotspots.Additionally, mapping of vegetation communities based on phenology could allow for the accounting of interannual and seasonal landscape-scale community responses to interacting environmental factors [67]

Figure 1 .
Figure 1.NDVI image of Arizona with a focus on the Santa Rita Mountains (pointed to in black box) and surrounding sky islands (circled), which are part of the Madrean Archipelago.The study site is centered on the Santa Rita Mountains, in southeastern Arizona, USA.The higher NDVI (Date: June 2007) values make the sky islands contrast with the Sonoran Desert lowlands.Urban areas are found within the magenta polygons.

Figure 4 .
Figure 4. Pre-monsoon (May; 4a-4c) and post-monsoon (September; 4d-4f) NDVI images for a wet year (2001) and a dry year (2002, especially in the spring).The difference from a multiyear NDVI average for May and September highlight below-average NDVI values (blue) for most of the study area in May, 2002.Some of the lower elevation areas show increased greenness patterns during the 2002 monsoon season (September 2002; orange patterns).Curvilinear spatial patterns represent riparian areas.May, 2001 May, 2002

Figure 5 .
Figure 5. Average NDVI trajectories for each vegetation type, 2000-2007, Santa Rita Mountains, AZ, USA.Trends in magnitude, timing, and shape of NDVI are different for each vegetation type and each year.A high elevation forest fire in 2005 is the reason for the marked decreases in NDVI values for Conifer-oak forest and the somewhat lower NDVI values for the Pine-oak woodlands.

Figure 6 .
Figure 6.(a) Mean and standard deviation of NDVI and (b) precipitation, maximum and minimum temperatures for the mesquite upland scrub vegetation type across seasons 2000-2007, Santa Rita Mountains, AZ.Effects of climate on the NDVI trajectories can be seen in a relatively wet year (2001; high NDVI) and the dry spring season (2002; low NDVI).(a)

Figure 7 .
Figure 7. Start of season (Start t ; Day of Year) and for (a) 2001 and (b) 2002, and Large Integral for (c) 2001 and (d) 2002 across the Santa Rita Mountains, AZ.Gray areas are agricultural or industrial land cover types and are excluded from analysis while white areas represent incomplete phenological retrievals for that year.

Figure 8 .
Figure 8. Small NDVI integral for 2001 (a) and 2002 (b) shows both the seasonal response and spatial patterns for the vegetation cline.The mid elevation vegetation types have larger small NDVI integral values than the high elevation vegetation types and the desert vegetation types.Gray areas are agricultural or industrial land cover types and are excluded from analysis while white areas represent incomplete phenological retrievals.(a) (b)

Figure 9 .
Figure 9. Spatial patterns of the coefficient of variation (COV) for the start of the season (SOS) (a) and the small NDVI integral (SI) (b) for 2001 through 2006.

Table 2 .
Results of stepwise regression of climate variables for each phenometric, for sites across the Santa Rita Mountains, AZ (n = 238; degrees of freedom listed with test statistics; R 2 values are adjusted to account for multiple variables).Negative associations are indicated by italics.* = data were log transformed; ** = data were square-root transformed.Both transforms were made to improve the normal distribution of the data.

Table 1 .
Summary of average elevation, temperature and precipitation values for each of the natural vegetation types along the Desert-Mountain gradient as described by the SWReGAP land cover project.

Table 3 .
Significant correlations (p ≤ 0.0081) between phenometrics and climate metrics for pine-oak woodland.Relationships that have negative R values are italicized.(Elevation ~2,204 m).

Table 4 .
Significant correlations (p ≤ 0.0081) between phenometrics and climate metrics for semi-desert grassland.Relationships that have negative R values are italicized.(Elevation ~1,357 m).