Mapping Temperate Vegetation Climate Adaptation Variability Using Normalized Land Surface Phenology

Climate influences geographic differences of vegetation phenology through both contemporary and historical variability. The latter effect is embodied in vegetation heterogeneity underlain by spatially varied genotype and species compositions tied to climatic adaptation. Such long-term climatic effects are difficult to map and therefore often neglected in evaluating spatially explicit phenological responses to climate change. In this study we demonstrate a way to indirectly infer the portion of land surface phenology variation that is potentially contributed by underlying genotypic differences across space. The method undertaken normalized remotely sensed vegetation start-of-season (or greenup onset) with a cloned plants-based phenological model. As the geography of phenological model prediction (first leaf) represents the instantaneous effect of contemporary climate, the normalized land surface phenology potentially reveals vegetation heterogeneity that is related to climatic adaptation. The study was done at the continental scale for the conterminous U.S., with a focus on the eastern humid temperate domain. Our findings suggest that, in an analogous scenario, if a uniform contemporary climate existed everywhere, spring vegetation greenup would occur earlier in the north than in the south. This is in accordance with known species-level clinal variations—for many temperate plant species, populations adapted to colder climates require less thermal forcing to initiate growth than those in warmer climates. This study, for the first time, shows that such geographic adaption relationships are supported at the ecosystem level. Mapping large-scale vegetation climate adaptation patterns contributes to our ability to better track geographically varied phenological responses to climate change.


Introduction
Temperate vegetation phenology demonstrates large spatiotemporal variations that are controlled simultaneously by current climatic variability and underlying plant heterogeneity.The latter is shaped by genetic differentiation of plants via adaptation to past long-term climate gradients.While much attention has been drawn to temporal variations of phenology in response to climate change [1,2], we have not been able to effectively account for spatial variations of phenological response that are caused by intra-and inter-specific differences [3,4].At the continental scale, geographic gradients of spring land surface greening that show gradually delaying timing towards colder locations are well-understood [5,6].However, within this observed variability, the portion determined by underlying genotypic difference (vs.immediate temperature effect) remains difficult to quantify [7,8].
Specifically, phenology of diverse plant populations not only varies phenotypically and temporally in response to current climate, but also genotypically and spatially due to long-term climate impacts, as reflected by differing climatic requirements of different genotypes for growth and development, such as varied amounts of thermal forcing (e.g., warm temperature accumulation or growing degree days) needed to trigger bud burst [9,10].From common garden studies, many woody species have shown lower forcing and higher chilling requirements for spring phenology within populations adapted to colder climates [11][12][13][14][15].A certain amount of chilling is necessary for many tree species to break winter dormancy [16,17].Grass species, however, appear to demonstrate the same forcing-based clinal variation without obvious chilling requirement [18].For certain woody species, lengthening photoperiod can partly compensate for lack of chilling or forcing [19][20][21].But the effect of photoperiod on spring phenology, in addition to being limited to a few species, is supplemental to temperature [22].Except with very mild winters, plants growing in their natural environments (vs. in common gardens or growth chambers) usually have their chilling requirements fulfilled early in the year, allowing thermal forcing to be the most effectual predictor of spring phenological timing.At the interspecific level, later flushing of species transplanted from warmer locations implies similar forcing-based geographic patterns [23].
Such cline-based geographic variations allow plants to maximize the use of available heat in varied climate conditions to initiate growth while minimizing frost risks in spring.Therefore we termed this phenomenon "growth efficiency" in a previous study, corroborating the findings from common gardens by comparing naturally adapted plants to a clonal reference [24].However, now realizing this name may create confusion regarding the processes being described, we are changing the name to the hopefully more straightforward "climatic adaptation" moving forward.These species-level relationships are the basis and rationale for us to explore corresponding geographic vegetation patterns at the ecosystem level.Given the challenges in directly measuring genotypic variations underlying vegetation heterogeneity, we employed an indirect approach to remove the effects of contemporary climate, by normalizing satellite-observed land surface phenology with a cloned plants-based phenological model.We assume that geographic patterns of normalized land surface phenology are reflective of genetic differences of vegetation, given that the influences of additional environmental factors (mainly soil types) are likely negligible compared to the overriding climatic effect [25], especially at large spatial scales.

Data Overview
The datasets used span a 30-year period from 1982 to 2011 and cover the conterminous U.S. Land surface phenology (LSP) start-of-season (SOS or greenup onset) dates were derived from satellite data collected by Advanced Very High Resolution Radiometer (AVHRR) and Moderate-Resolution Imaging Spectroradiometer (MODIS).A cloned plants-based Spring Indices-extended (SI-x) first leaf (FL) model was used to approximate genetically uniform phenological responses to climate.The specifics of these datasets are provided below.
The AVHRR and MODIS EVI2 time series were first smoothed using an algorithm based on normal development of the annual vegetation growth cycle [32][33][34].Temporal EVI2 trajectory was then reconstructed using a Hybrid Piecewise Logistic Model [35].SOS dates were identified according to the rate of change in curvature of reconstructed EVI2 temporal trajectories during the spring to summer period.Specifically, an SOS date corresponded to the day of year (DOY) on which curvature change rate in a fitted EVI2 curve exhibited a local maxima [36].In addition, if the annual variation in EVI2 was minimal (<0.04), such as in evergreen forests or deserts, SOS dates were not retrieved.Moreover, if the annual EVI2 time series was very flat (as in some shrublands), a very early SOS date may be produced.In such cases, a 5% threshold was used instead to determine SOS.

Modeled Phenology Based on Cloned Plants
To simulate a genetically-uniform response pattern of phenology that represent the portion of spatial variation caused solely by contemporary climate, we used first leaf predictions from a cloned plants-based Spring Indices-extended (SI-x) model [37].SI-x models are an extension of the original Spring Indices (SI-o) [38], which include three sub-models built respectively on phenologies of clone lilac Syringa chinensis "Red Rothomagensis", clone honeysuckle Lonicera tatarica "Arnold Red", and clone honeysuckle L. korolkowii "Zabeli".SI-o models predict both first leaf and first bloom using a multi-step regression driven by growing degree hours during specific time windows and high-energy synoptic events after the date of chilling satisfaction [39].The SI-x models were developed with the same methods as SI-o with the exception that chilling is not calculated, and instead, the starting date of thermal accumulation (leading to first leaf) is always set as 1 January for all stations/years [40].This adjustment allowed model predictions to be extended into the subtropics to cover southern parts of the continental U.S. The SI-x outputs showed minimal difference when compared with the SI-o models and with lilac phenological data over the region for which both models produce outputs.
While SI models do not reproduce all the details of multi-species phenological data at any site, or the specific phenology of some types of plants, these models integrate changes in spring temperatures into indices directly related to the growth and development of many temperature-responsive (but not water-limited) tree, shrub, and herbaceous plant species [37,38,41].Spatial patterns outlined by SI model predictions reflect general bioclimatic differences rather than being tied to the specific cloned plants used.Therefore, the SI-x models have been widely adopted as standardized measures of the start of the growing season (similar to what the Palmer Drought Index provides for evaluating moisture stress).They are employed as one of the "Climate Change Indicators in the United States" by the U.S. Environmental Protection Agency [42], and proposed as a part of the "Pilot National Climate Indicators System" under development by the National Climate Assessment within the U.S. Global Change Research Program [43].
Along with LSP SOS, we produced SI-x first leaf (FL) dates from 1982 to 2011, and at the same spatial resolution of 0.05 degree.For each year, temperature data from 678 to 771 weather stations (the total number of stations available varied by year) for the conterminous U.S. were used to run the model.FL dates predicted for all station locations were then interpolated to the entire region using an Empirical Bayesian kriging method (ArcGIS, ESRI, Redlands, CA, USA).An LSP SOS normal and an SI-x FL normal were calculated by averaging the two 30-year data layer time series, respectively.

Rationale and Algorithms
LSP SOS shows spatial patterns of vegetation phenology that are dependent on both climate variation and plant heterogeneity, the latter contributed by both intra-and inter-specific differences of plants.SI-x FL dates are a function of climate variables (mainly thermal forcing) only, and therefore insensitive to genetic differences of vegetation.Hence we inferred the portion of variation induced by vegetation heterogeneity through normalizing the observed LSP SOS with modeled SI-x FL.The normalized LSP is represented by a phenological index that approximates genetically-controlled (vs.climate-triggered) spatial variations in vegetation phenology.We term this phenological index "climatic adaptation index" given that the geographic patterns it depicts primarily result from adaptation of plants to long-term climate gradients.
Three algorithms were used to derive the climatic adaptation index (CAI).These algorithms are based on the same rationale of separating instantaneous climate effect (FL) from observed patterns (SOS), but offering different mathematical solutions to normalization.We assumed that the instantaneous climate and genotypes contribute to observed patterns in ways of being additive, interactive, or both, from which the three algorithms are derived, respectively.Allowing these alternatives provides opportunities to ascertain the pros and cons of each in estimating the climatic adaptation patterns.Therefore at the risk of redundancy, given that the current work is limited to developing the initial methodology of CAI, we employed the three algorithms in tandem, in the hope that an initial comparison can be made, and future studies will shed more light on the utilities of the indices, when they are tested in additional applications.
Assuming an additive relationship between climate and genotypes, the first type of CAI (CAI DIF ) was developed by subtracting SI-x FL from LSP SOS (1).This transformation retains the units of SOS and FL in days.A larger value of CAI DIF at a given geographic location corresponds to longer time and more forcing needed for greenup onset at that location.It should be noted that because SI-x FL is based on plants different than that underlying vegetation observed by LSP SOS, absolute values of CAI do not carry an explicit meaning.But the spatial pattern does reflect normalized phenological differences; therefore it is the relative difference of CAI values over space that matters.
Further, interaction between climate and plant genotypes may also be approximated with a multiplicative relationship.Hence, a second type of CAI (CAI DIV ) is calculated by dividing SOS with FL (2).Resultant values carry similar utility and limitation as those from the subtraction approach, but they are dimensionless quantities.A larger CAI DIV value corresponds to a relatively later phenology and a higher forcing requirement of vegetation for greenup onset.
A third algorithm combines additive and interactive (multiplicative) relationships to produce a CAI (CAI ND ) with a normalized difference equation.The difference between SOS and FL is normalized by dividing the sum of SOS and FL (3).Thus, resultant values are ratios in the range of ´1 to 1.This approach maximizes contrast between the two measurements and increases the overall consistency by reducing local base variations across a study region.CAI ND shares the similar connotations of the two previous types in that a higher value at a location corresponds to a relatively later phenology caused by a higher forcing requirement.
The algorithms of CAI derivation were implemented to 30-year (1982-2011) normals of LSP SOS and SI-x FL for the conterminous U.S.But large areas in the west and southwest of the U.S. were omitted for lack of valid predictions, given that SI-x FL and LSP SOS are mostly applicable to temperature-sensitive and moist ecosystems.Hence, we further extracted latitudinal gradients of CAIs for only the humid temperate domain in the eastern U.S. [44].The latitudinal gradients were calculated through averaging CAI values in that region by latitude at 0.05-degree intervals.

Results
The 30-year normal LSP SOS transitions from earlier dates (smaller DOY) in the south to later dates (larger DOY) in the north (Figure 1a).This well-known latitudinal gradient is most prevalent in the eastern U.S. and disrupted by altitudinal effects towards the western mountains.The Appalachians and agricultural fields in the Midwest are identifiable by regional variations they contribute.In the arid Southwest where vegetation is sparse and leaf flush is primarily rain triggered, very few valid SOS values are present.The 30-year normal of predicted SI-x FL follows the same geographic patterns (Figure 1b), but demonstrates a much smoother and more continuous spatial pattern compared to that of the LSP SOS.This reflects SI-x FL representing phenological variation induced by current climate only and not including contributions from the underlying vegetation heterogeneity.In addition, LSP SOS is systematically later than SI-x FL at most locations.This has led to predominantly positive values for CAI DIF and CAI ND , and values > 1 for CAI DIV , but does not affect the geographic variations they represent, as only the relative differences across space matter.
Climate 2016, 4, 24 5 of 11 eastern U.S. and disrupted by altitudinal effects towards the western mountains.The Appalachians and agricultural fields in the Midwest are identifiable by regional variations they contribute.In the arid Southwest where vegetation is sparse and leaf flush is primarily rain triggered, very few valid SOS values are present.The 30-year normal of predicted SI-x FL follows the same geographic patterns (Figure 1b), but demonstrates a much smoother and more continuous spatial pattern compared to that of the LSP SOS.This reflects SI-x FL representing phenological variation induced by current climate only and not including contributions from the underlying vegetation heterogeneity.In addition, LSP SOS is systematically later than SI-x FL at most locations.This has led to predominantly positive values for CAIDIF and CAIND, and values > 1 for CAIDIV, but does not affect the geographic variations they represent, as only the relative differences across space matter.The three CAIs show very similar geographic patterns that are traced by latitudinal and altitudinal variations (Figure 1c-e).Notably, the direction of CAI gradients is opposite to that of LSP SOS and SI-x FL, i.e., higher CAI values are found in southern and warmer areas, and lower values at colder locations of higher latitude or higher elevation.The latitudinal gradients of CAIs are most obvious, especially in the eastern humid temperate domain (Figure 1f).Lower CAI values in the north reflect relatively earlier SOS in reference to FL, and correspondingly lower forcing requirement of vegetation than in the south.More specifically, the continuous latitudinal gradients in the eastern U.S. are interrupted by slightly higher CAI values in the central Midwest region, affected by wide spread agricultural fields.Towards the western mountains, lower CAI values are common for vegetated areas at higher elevations.In the relatively low-lying Appalachians, less distinct (in The three CAIs show very similar geographic patterns that are traced by latitudinal and altitudinal variations (Figure 1c-e).Notably, the direction of CAI gradients is opposite to that of LSP SOS and SI-x FL, i.e., higher CAI values are found in southern and warmer areas, and lower values at colder locations of higher latitude or higher elevation.The latitudinal gradients of CAIs are most obvious, especially in the eastern humid temperate domain (Figure 1f).Lower CAI values in the north reflect relatively earlier SOS in reference to FL, and correspondingly lower forcing requirement of vegetation than in the south.More specifically, the continuous latitudinal gradients in the eastern U.S. are interrupted by slightly higher CAI values in the central Midwest region, affected by wide spread agricultural fields.Towards the western mountains, lower CAI values are common for vegetated areas at higher elevations.In the relatively low-lying Appalachians, less distinct (in comparison to western mountains) orographically-induced CAI variations are blended into the more dominant latitudinal gradients.
Extracted north-south (latitudinal) gradients of CAIs within the humid temperate domain in the eastern U.S. (cf. Figure 1f) are evident (Figure 2).All CAIs demonstrated gradually decreasing values with increasing latitude.Further, CAI decline with latitude is relatively rapid up to 38 ˝N; afterwards the change becomes small.Among the three CAIs, CAI DIF is most noisy and shows a sharper transition across the 38 ˝N parallel.CAI DIV is most smooth and resembles an exponential decline process that can be modeled with a logistic function [45].CAI ND integrates features of the previous types and additionally provides a wider range of CAI variability as shown by a slightly steeper slope.The CAI ND gradient also shows a piece-wise response divided by approximately the 38 ˝N parallel.Interestingly, this latitude coincides with the boundary of mesothermal (above-freezing winter) and microthermal (sub-freezing winter) climates [46], and is around a transitional zone between areas with different long-term phenological trends [47].Extracted north-south (latitudinal) gradients of CAIs within the humid temperate domain in the eastern U.S. (Cf. Figure 1f) are evident (Figure 2).All CAIs demonstrated gradually decreasing values with increasing latitude.Further, CAI decline with latitude is relatively rapid up to 38°N; afterwards the change becomes small.Among the three CAIs, CAIDIF is most noisy and shows a sharper transition across the 38°N parallel.CAIDIV is most smooth and resembles an exponential decline process that can be modeled with a logistic function [45].CAIND integrates features of the previous types and additionally provides a wider range of CAI variability as shown by a slightly steeper slope.The CAIND gradient also shows a piece-wise response divided by approximately the 38°N parallel.Interestingly, this latitude coincides with the boundary of mesothermal (above-freezing winter) and microthermal (sub-freezing winter) climates [46], and is around a transitional zone between areas with different long-term phenological trends [47].

Discussion
All three CAIs with their algorithms and results were reported to show alternative ways of normalizing land surface phenology, albeit the similarity of patterns revealed.Both additive and interactive (multiplicative) relationships appeared to have mathematically separated current climate effects from the primarily genotypic influences on spatial variation in vegetation phenology.The difference-based CAI (CAI DIF ) retains the unit in days, which more directly conveys ecological information as of the number of days the normalized phenology is different from one location to another.While the CAI based on (CAI DIV ) yielded a smoother transitional pattern than CAI DIF , the normalized difference CAI (CAI ND ) appeared to include information from both previous CAIs.Therefore, we generally recommend CAI ND for potential future use as far as the units are not of concern for intended applications.However, this does not rule out the utility of CAI DIF and CAI DIV , which may turn out to be more straightforward in particular cases.We leave the decision to potential users in consideration of their specific research questions or needs.
In this study, we derived CAIs primarily based on thermal forcing, given that SI-x models do not take into account photoperiod effects and adopt a uniform chilling fulfillment date (i.e., 1 January) for all locations.The smaller latitudinal gradient beyond the 38 ˝N parallel may indicate relatively higher chilling requirements of the northern vegetation which could have "held the reins" of spring phenology progression through counteracting with thermal forcing [48,49].This hypothesis, however, requires further testing with additional data and models.Future CAI development can be improved by explicitly incorporating these additional climatic factors.Nevertheless, thermal forcing is the dominant factor affecting spring phenological timing for the majority of tree species [22,23], and even more so for grass and understory species [18,50].Hence, under the current climate, neglecting occasional mild winters that could hamper chilling requirement fulfillment in limited regions, the forcing-based approach likely captures the first-order variability in temperate vegetation phenology that is caused by climate adaptation.
The influence of climate on vegetation operates at two very different temporal scales: The contemporary variability that is mostly accommodated by phenotypic plasticity, and long-term averaged conditions that have shaped the genotypes through adaptation [51,52].In spite of the similar spatial gradients, the short-term and long-term climatic effects on vegetation phenology take two distinct pathways (i.e., environmentally and genetically) and should be assessed separately.It is the latter that is often neglected in describing and accounting for phenological patterns over space and is here approximated using the derived climatic adaptation indices.Earlier studies comparing weather-predicted phenology and satellite-derived phenology noticed the fundamental difference between the two types of measurements [53,54], which is primarily that current ecosystem-scale phenological models (e.g., Spring Indices) assume a genetically-uniform phenological response to climate, while satellite observations capture the integrated effect of both climatic and genetic differences.
Normalized land surface phenology revealed geographic patterns that reflect inherent differences of vegetation resulting from adaptation to past climate.Given that the normalization removed spatial variation caused solely by phenotypic plasticity (as in cloned plants), the derived CAI maps can be analogously understood as showing how greenup onset date of vegetation would vary across the country if a uniform contemporary climate existed everywhere, mimicking the situation in a common garden, only applied more comprehensive (i.e., targeting vegetation as a whole vs. particular species).Under this analogy, due to lower thermal requirements, northern vegetation would initiate growth earlier than vegetation in the south, opposite to the annually observed latitudinal gradient of vegetation phenology.As widely documented in common garden studies [23,55,56], northern populations of a species may show either earlier or later spring phenology than southern populations [10,12,13,18,48].The former case is referred to as a countergradient variation in evolutionary biology, as the gradient of common garden (or normalized) phenology is opposite to that of in situ (or locally-observed) phenology; and the latter case a cogradient variation [57,58].A countergradient variation (north being earlier, as in the Northern Hemisphere) is reflective of different thermal forcing requirements of populations, while a cogradient variation (south being earlier) of different chilling requirements [59].Given temperature changes simulated by common garden experiments via transferring seeds across large distances from their source locations are often substantial (e.g., 10-20 ˝C), chilling-controlled gradient could thus be manifested.But the magnitude of current climate change (e.g., <2 ˝C) is unlikely to cause widespread chilling deficits for temperate plant phenology, hence a forcing-based countergradient variation is likely dominant in the present and upcoming decades.This study for the first time mapped this underlying countergradient pattern of temperate vegetation phenology at the ecosystem level.
Such ecosystem level patterns indicate that temperate vegetation grown in colder climates tends to require less thermal forcing for greenup onset than that adapted to warmer climates.A crucial implication is that temperate vegetation at different latitudes may respond to the same degree of climate change differently.Mean bud burst dates of multiple temperate U.S. tree species have been predicted to advance more in the north than the south under projected climate scenarios [60].An underlying cause of this predicted latitudinal difference of phenological response is likely the climatic adaptation gradient unveiled in this study.Hence, development of climatic adaptation indices and mapping of the intrinsic variability in climatic requirements of vegetation phenology is significant in terms of enabling monitoring and modeling of geographically-explicit responses of vegetation phenology to climate change.As demonstrated in this study, such applications can be potentially done at large spatial scales using satellite remote sensing.This makes it operationally feasible to build continental scale phenoclimatic models that are able to account for vegetation heterogeneity in tracking climate change impact on ecosystems.
Furthermore, climate-normalized phenology is additionally useful for tracking vegetation changes under other environmental stresses, such as forest fire and insect outbreak.This is because a climate-normalization, when implemented on an annual or decadal basis, helps remove spatial variations in observed vegetation phenology caused by short-term weather fluctuations.The derived CAIs, therefore, may more directly reflect variations in inherent vegetation properties, and when used in change detection, can better reveal underlying ecosystem changes.The U.S. Forest Service has used multi-year land surface phenology-based phenoregions to track anomalous phenological changes in relation to forest health [61][62][63][64].By implementing normalized land surface phenology with related climatic adaptation indices in monitoring tasks as such, we could better target vegetation changes related to processes such as forest succession or ecological disturbances.
Finally, the current study covered the conterminous U.S., with a focus on the eastern humid temperate ecological domain.Future work is needed to test the CAIs in other land regions around the world.We expect that similar patterns can be found, especially for vegetation in the humid temperate climates.Therefore similar applications of CAIs for tracking geographically-explicit impacts of climate change on vegetation phenology and other ecological transitions may be promising in additional temperate land regions.

Conclusions
In this study, we derived normalized land surface phenology over the conterminous U.S. that is independent of current climatic influence and reflects underlying vegetation heterogeneity as shaped by long-term climatic adaptation.The resultant climatic adaptation indices (CAIs) revealed geographic adaptation patterns of temperate phenology at the ecosystem level.A countergradient variation in spring phenology by latitude is particularly evident for temperate vegetation in the eastern U.S., suggesting that colder-climate adapted vegetation requires less thermal forcing for greenup onset than warmer-climate adapted vegetation, embodying the dominant clinal relationships within and among species at the ecosystem level.Such intrinsic geographic variability of vegetation phenology due to climatic adaptation has not been mapped previously.Normalized land surface phenology with climatic adaptation indices is potentially useful for tracking spatially-explicit phenological responses to climate change and additional underlying changes in vegetation.

Figure 1 .
Figure 1.Phenological variations in the conterminous U.S. (a) Land Surface Phenology (LSP) Start-of-Season (SOS) in day of year (DOY) based on 30-year (1982-2011) normal; (b) Spring Indices-extended (SI-x) first leaf (FL) in DOY based on 30-year (1982-2011) normal; (c) Climatic adaptation index (CAI) based on the difference between the SOS and FL (CAI DIF ); (d) CAI based on the division of SOS and FL (CAI DIV ); (e) CAI based on the normalized difference between SOS and FL (CAI ND ); (f) U.S. Forest Service ecological domains.
mountains) orographically-induced CAI variations are blended into the more dominant latitudinal gradients.

Figure 2 .
Figure 2. Latitudinal gradients of averaged climatic adaptation indices (CAIs) for the eastern Humid Temperate Domain in the U.S. (Cf. Figure 1f).The three CAIs are respectively developed through difference (DIF), division (DIV) and normalized difference (ND) of observed and modeled phenologies.

Figure 2 .
Figure 2. Latitudinal gradients of averaged climatic adaptation indices (CAIs) for the eastern Humid Temperate Domain in the U.S. (cf. Figure 1f).The three CAIs are respectively developed through difference (DIF), division (DIV) and normalized difference (ND) of observed and modeled phenologies.