Climatic Regulation of Vegetation Phenology in Protected Areas along Western South America

: Using 19 years of remotely sensed Enhanced Vegetation Index (EVI), we examined the effects of climatic variability on terrestrial vegetation of six protected areas along southwestern South America, from the semiarid edge of the Atacama desert to southern Patagonia (30 ◦ S–51 ◦ S). The relationship between satellite phenology and climate indices, namely MEI (Multivariate ENSO Index), PDO (Paciﬁc Decadal Oscillation) and SAM (Southern Annular Mode) were established using statistical analyses for non-stationary patterns. The annual mode of phenological activity ﬂuctuated in strength through time from the semiarid region to the border of southern Patagonia. Concomitantly, enhanced synchrony between EVI and climatic oscillations appeared over interannual cycles. Cross correlations revealed that variability in MEI was the lead predictor of EVI ﬂuctuations over scales shorter than 4 months at lower latitudes and for the most poleward study site. The PDO was correlated with EVI over lags longer than 4 months at low latitude sites, while the SAM showed relationships with EVI only for sites located around 40 ◦ S. Our results indicate that the long-term phenological variability of the vegetation within protected areas along southwestern South America is controlled by processes linked to climate indices and that their inﬂuence varies latitudinally. Further studies over longer time scales will be needed to improve our understanding the impacts of climate change on vegetation condition and its effect over phenological variability.


Introduction
Climatic variability has been identified as the lead driver of observed changes in phenological cycles worldwide [1][2][3]. Phenological changes have consequences that cascade across the entire ecosystems: from altering the timing of flower emergence, and their mismatch with pollinators life cycles, to an earlier onset of greening and delayed senescence of forest canopies [4][5][6][7]. Climate-driven drought disturbs the phenology of forests in boreal, temperate, and tropical ecosystems, impacting tree growth and mortality [8,9]. On the other hand, the strong warming-induced trend for earlier spring and later autumn has changed the photosynthesis-mediated carbon uptake in the temperate forests of eastern North America during the past two decades [10]. Although changes in phenological patterns have been reported mostly for the northern hemisphere [11,12], the impacts of climate variability on the dynamics of vegetation phenology have a global reach and are expected to display a complex structure over space and time [13].
Understanding the association between climatic change and vegetation phenology requires long-term evaluation of the spatiotemporal interplay between oceanic conditions and land vegetation characteristics. Capturing vegetation characteristics solely based on in situ measures is challenging due to the large spatial heterogeneity across the scales [14]. To circumvent this shortcoming, spectral indices derived from remotely sensed platforms have been employed as a proxy to monitor annual and interannual variability from the field [15] to regional [16] and global scale [17]. With the proliferation of satellites and radiometric measurements, several vegetation indices have been developed among which the Enhanced Vegetation Index (EVI) and the Normalized Difference Vegetation Index (NDVI) have been more highlighted for their enhanced performance in capturing plant phenology [18,19].
In hydro-climatological studies, oceanic conditions and oscillations have been presented in the form of large-scale climate indices that capture patterns of heat and mass distribution in the atmosphere and oceans [20,21]. Among climate indices, ENSO (El Niño-Southern Oscillation) has been recognized as one of the leading modes of global climate variability, and other climatic modes, such as the Pacific Decadal Oscillation (PDO) or the Atlantic Multidecadal Oscillation (AMO) interact with ENSO to modulate phenological patterns of vegetation at a global scale [22,23].
The temporal dynamics of climate indices characterize the state of the coupled oceanatmosphere systems across the entire Pacific ocean, which is tightly coupled to global-scale patterns of climatic variability, from interannual to decadal scales [24].
Using the spectral indices researchers have studied climate indices and vegetation phenology. For example the authors of Zhao et al. [25] have shown ENSO (El Niño-Southern Oscillation) as a key driver of interannual variability of vegetation phenology using NDVI. Besides climatic effects, the influence of precipitation and temperature on phenological activity has been studied in different regions (e.g., arid, semiarid, tropical) [26,27]. In the semi-arid climate of the middle east, the agricultural crops are severely stressed by heat and limited rainfall [28] while in Mediterranean areas, vegetation phenology is stressed during summer due to soil moisture depletion [29] and in southern Africa and South America the growth of vegetation is limited by water availability [30]. In temperate forest, the increase in air temperature has also altered vegetation phenological patterns (e.g., phenological advances of approximately 3-8 days for each • C increase of 1), vegetation activity and regional carbon cycling [10,31].
The impact of climate change is particularly important and more pronounced over biodiversity-rich areas [32]. The coast of Chile, extending along western south America from the Atacama desert to Patagonia region, spans 12 vegetation formations and 127 vegetation types [33]. This broad latitudinal gradient provides a model system to evaluate the impacts of climate variability on the spatio-temporal dynamics of vegetation phenology. Ecological studies have recognized the importance of climatic variability as a driver of changes in the productivity and composition of the different vegetation types found along this latitudinal gradient [34,35]. For example, long-term monitoring in semiarid central Chile has shown that vegetation dynamics are tightly linked to interannual changes in the ENSO-controlled precipitation regime [36].
Moreover, an ENSO-driven drought now extending for over a decade (2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020) in central Chile has imposed different trends on terrestrial productivity depending on latitude [37]. Similarly, the coupling of positive ENSO and SAM anomalies seem to underpin a pattern of increased forest fires across central and southern Chile over the same period during the last decade [21,[38][39][40]. Together, these studies suggest that long-term climatic variability is impacting seems to modulate vegetation dynamics along western south America and that their effects are not homogeneous in space neither in space nor in time. However, the temporal climatic dynamics driving observed changes in vegetation dynamics have not been investigated to date.
Satellite-derived vegetation indices have been used before to assess the temporal variability of different vegetation ecosystems of Chile [29,30]. However, their ability to capture phenological heterogeneity and their linkages with large-scale climatic oscillations are currently lacking.
In this study, we evaluated the sensitivity of EVI to climatic variability along an extensive region characterized by heterogeneous biomass distribution and hypothesized that climatic oscillations have an asymmetric influence on the spatiotemporal phenological variability across the study region. To this end, the aim of this study was to comprehensively examine the effects of climatic variability on the EVI over multiple temporal scales, with an emphasis on the impacts that the time-varying climate has on the phenological patterns along the meridional extent. Section 2 introduces the study area and satellite data. Section 3 describes the main results evaluating the terrestrial plant phenology response to climatic oscillations along the latitudinal gradient. The discussion is presented in Section 4, emphasizing the differences in phenological responses to climate change. Finally, the main conclusions are highlighted in Section 5.

Study Area
The central-southern coast of Chile along the western Pacific coast of South America is characterized by a narrow semiarid region, where the Atacama desert, the driest on the planet, changes into a Mediterranean climate (30 • S-37 • S). While humidity in the semiarid region is provided mostly by seasonal coastal fog, the Mediterranean region is characterized by seasonal precipitation (>200 mm) concentrated during the austral winter months (June, July, August), which is subject to considerable interannual variability. On the other hand, the southern-austral region experiences a temperate climate (38 • S-53 • S), with year-round rainfall that can exceed 2000 mm, particularly south of 40 • S [41][42][43].
The vegetation of the lower-latitude part of the coastal temperate sector, between 38 and 43 • S, corresponds to Valdivian rainforest, dominated by species with wide leaves adapted to high humidity [44]. The Sub-Antarctic forest poleward of 47 • S is dominated by deciduous Nothofagus species such as Nothofagus pumilio and Nothofagus antarctica, among others [44,45]. Due to the high endemism and its isolation, the Valdivian and Sub-Antarctic forests integrate a global biodiversity hotspot [46,47], which has been affected by human disturbances (selective logging and livestock grazing) and natural processes [48].
We focused our study on forested areas located within the boundaries of National Parks to tease out the potential effects of long-term human activities on vegetation (i.e., changes in land use and/or irrigation) from climate-driven phenological dynamics. In Chile, the conservation strategy that defines the legal boundaries and management plans for these protected areas follows that stated by the International Union for Conservation of Nature (IUCN) as: "a clearly defined geographical space, recognized, dedicated and managed, through legal or other effective means, to achieve the long-term conservation of nature with associated ecosystem services and cultural values". These areas generally include heterogeneous vegetation and land surface cover, for example rocky or glaciated areas, and are managed under the National System of Protected Wild Areas (SNASPE), which is administrated by the Chilean National Forestry Corporation (CONAF) [33,49]. We selected six national parks that spanned the vegetation gradient from the semi-arid edge of the Mediterranean sector to the Cold Patagonian tundra, using cartography provided by the SNASPE to define the units of study (Figure 1). The selected parks were: Fray Jorge (30 • S), La Campana (32 • S), Alerce Costero (40 • S), Tantauco (43 • S), Laguna San Rafael (47 • S), and Torres del Paine (51 • S) ( Table 1).

Phenology and Climatic Data
We used satellite images of the Enhanced Vegetation Index (EVI; MOD13Q1 product, version 6) from the Moderate Resolution Imaging Spectroradiometer (MODIS) sensor onboard the Terra satellite. EVI is an L3 product with a resolution of 250 m, generated every 16 days from daily images. The EVI algorithm choose the best available pixel value (i.e., low clouds, low view angle, and the highest NDVI/EVI value) from all the acquisitions from the 16 days period. We downloaded and analyzed a time series of images from February 2000 to December 2018. Using the polygon of each site, we generated an internal 250-m buffer for each park and the corresponding EVI time series (e.g., [50]). Our analysis aimed to capture biospheric responses to climate, and hence they included the totality of the park area and the different kinds of vegetation within (see Table 1). Zonal statistics per observation (scene) were calculated for each park, at each time. This results in a temporal series of one data point every 16 days, for each park (the mean). That data were aggregated by month using a weighted mean, in order to match the temporal resolution of the climatic data, with the weight consisting in the number of valid pixels per scene at each time and park. When there was only one scene per month, that observation is used. When there was two, the weighted mean is computed. All data extracting and processing were implemented in Python and the MODIS imagery re-projection and cropping were carried out using the Modis Reprojection Tool.
In order to identify the long-term relationships between climatic oscillations and phenological variability, we examined three climatic indices: MEI, PDO and SAM, which were obtained from the National Weather Service website (www.cpc.ncep.noaa.gov, accessed on 1 September 2020) of the National Oceanic and Atmospheric Administration (NOAA). These three monthly climatic indices strongly influence climate along western south America, and their variability modulates changes in precipitation and streamflow along the study region [51,52].

Data Analysis
We performed a wavelet analysis on MODIS phenological data to explore the interannual variability in phenological indices. Additionally, a bootstrapping method was used to evaluate the significance of wavelet spectra (more details are found in [53]). The relative importance of frequencies for each time step is represented in the time-frequency plane to form the local wavelet power spectrum (WPS). Here, we used the Morlet wavelet to detect the characteristics in phenology fluctuations, which according to Cazelles et al. [54] is defined as: where exp(−iω 0 t) is the product of a complex sinusoidal, π −1/4 is a normalization factor to ensure unit variance and ω 0 represent the central angular frequency of the wavelet [54].
One of the key advantages of a continuous wavelet function is the mathematical relationship between the wavelet scale and its frequency a and its frequency f . We also computed the global WPS where all wavelet coefficients of the same frequency f were averaged. The global WPS is defined by Cazelles et al. [54] as: where σ 2 x is the variance of time-series (EVI), T is the time duration of the time series (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018) and W x ( f , τ) is the wavelet transform of the x(t) . All analyses were performed using the MATLAB wavelet scripts (www.biologie.ens.fr/~cazelles/bernard/Research.html, accessed on 6 April 2021). To quantify the co-variation between EVI and climate indices we computed the wavelet coherence, which measures the correlation between the spectra of two series [55] and is defined as: where <> indicates a smoothing operator in time and scale, W x,y ( f , τ) is the wavelet co-spectrum of the two time series x(t) and y(t), respectively, W x,x ( f , τ) and W y,y ( f , τ) are the wavelet transform of the x(t) and y(t). The wavelet coherence (R x,y ( f , τ)) is equal to 1 when there is a perfect linear relation between the frequency of two temporal signals, and is equal to 0 when signals x(t) and y(t) are independent with no common frequencies at a particular temporal scale [54,56]. Finally, we computed cross-correlations between the EVI signal and the climatic indices for each park with time lags from 0 to 7 months. The semi-annual lag was chosen as the maximum delay where the impact of climate-driven seasonal variability could be detected on the phenological cycle.

Results
Vegetation productivity, using EVI as proxy, showed a significant interannual variability in the intensity of the 1-year period-the phenological band. In general, WPS of EVI revealed a unique and significant annual signal, that persisted throughout the study period for all evaluated parks (Figure 2 (Figure 6). In several portions of the analysis, intra-anual (0.5-year) modes were significant, yet they should be interpreted with caution due to their short duration and the coarse temporal resolution of our data.     [56]). Figure 7. Climate indices lead the analysis and positive correlations mean that the corresponding climate index and EVI are in phase (e.g., positive MEI before leads to positive EVI after). At Fray Jorge Park, EVI showed a quick and weakly significant positive relationship with MEI that persisted for 3 months. For PDO, the response was negative at lag-0 and switched to a positive relationship over lags between 6 and 7 months. No significant cross-correlations were observed with the SAM at this site (Figure 7a). A similar pattern was observed in La Campana National Park, but the positive cross-correlation between EVI and MEI persisted for 5 months and the positive significant lags with PDO appeared at 6 and 7 months (Figure 7b). Alerce Costero and Tantauco National Parks, located ca. 10 degrees south of La Campana, did not share the lagged correlation structure with the Mediterranean region. No significant correlations with MEI were observed, but there were significant negative associations with PDO between 2-and 5-month lags for Alerce Costero and 4-month lag for Tantauco, which also showed a significant and positive zero-lag correlation with SAM (Figure 7c,d, respectively). In the case of San Rafael National Park, significant and positive correlations were observed at lag-0 for both PDO and SAM, and also for SAM at 1 month (Figure 7e), while for Torres del Paine, the most austral National Park, EVI showed significant positive correlations with MEI at 3-4 months and negative at 0 and 4-5 months with PDO (Figure 7f).

Discussion
Our work shows a persistent and time-varying association between phenological patterns and different climatic oscillations over a broad range of vegetational formations spread along the western coast of south America. These results provide new insights on how the spatio-temporal variable natural climate forcing modulates the physiological status in terrestrial ecosystems over global, regional and local scales (see references in [58,59]) and influences the provision of essential environmental variables.
Due to the critical role that vegetation plays in the global biogeochemical and hydrological cycles, and the services we derive from them, it is relevant to understand how plant dynamics respond to global climatic oscillations over multiples scales [60]. Our study highlights the impacts that tropical ocean-atmosphere dynamics have on seasonal vegetation activity, ranging from the semi-arid Mediterranean scrub down to the sub Antarctic tundra of Patagonia. Previous studies evidence the strong impact of ENSO events (e.g., extreme MEI anomalies) on vegetation of the Southern Hemisphere [34,51]. The warm phase of ENSO (El Niño) covaries with severe drought events in some regions (e.g., Australia, southeast Asia, northeastern South America), whereas the cold phase (La Niña) is associated with droughts in western south America [22,61]. In concomitance with ENSO, PDO also appears as a driver of the phenological patterns along western south America, mediated by the meridional displacement of the subtropical gyre, changes in the pattern of the winds, and SST anomalies [43,62,63].
Changes in the precipitation pattern associated with negative phases of PDO are related to a drying trend in central Chile (30 • S-40 • S) over the last 30 years [42,43,64]. The detrimental effects of the reduced rainfall are evidenced as changes in the vegetation growing-season during the so-called megadrought period (2010-2018) in central Chile [37,64]. The uninterrupted sequence of dry years, i.e., years with a rainfall deficit, coincided with a negative EVI signal over pasture lands and areas with native forests cover. When coupled with climate effects (e.g., warming), these intense drought periods could promote fire activity (frequency and severity), mainly in areas with extensive plantations of non-native forests, native sclerophyll forests, and scrublands, vegetation types that largely determine the primary productivity patterns in central Chile [37,65]. While in northern Patagonia, both coastal ocean and land biological processes are forced by climatic variability over multiple temporal scales [66].
In the Southern Hemisphere, where instrumental records are sparse and relatively short, the opportunities for understanding the large-scale dynamical mechanisms driving climate variability are limited [67]. According to Marshall [68], the SAM oscillation is defined by the gradient between the high-pressure belt at medium latitudes and the lower pressures around coastal Antarctica. As such, it is calculated as the zonally averaged pressure difference between 40 • S and 65 • S. SAM oscillations are phase-dependent with ENSO, e.g., when a positive phase of ENSO (El Niño) coincide with a negative phase of SAM ( Figure A2), climate anomalies dominate south of 40 • S [69], generating extreme conditions due to high solar irradiation and reduced precipitation [66,70]. Such phaselocking could drive part of the observed changes in phenological activity, e.g., weak amplitude (see Figure 7).
Despite the low but significant lagged-effect of ENSO and PDO on the phenological pattern for the 43 • S-52 • S zone, our results strengthen the conclusion that the teleconnection between ocean-atmosphere oscillations are evident in the coast of South America. Multiple studies have associated ENSO with spring rainfall in southern China [71] and across western North America [72]. Although these climatic drivers have not been formally established, following a lack of in situ data, our results can be used to assess the impacts of climatic teleconnection processes with local phenology patterns.
EVI provides a consistent measure of greenness and improves sensitivity over dense vegetation conditions closely related to the seasonal variation of biological activity and their sensitivity to climatic oscillations over multiple time-lag relationships [73]. Throughout the study period (2000-2018), we observed significant time-lagged effects of the local vegetation to global climatic oscillations. However, in Fray Jorge National Park, several long-term studies that investigated the effects of climatic forces on population dynamics of wild rodents suggest that pulses of higher rainfall, induced by ENSO, impact ecological processes over multiple scales of organizations (e.g., [74][75][76]), including human communities and land use. Garreaud et al. [37] showed the effects of climate variability on rainfall along a latitudinal gradient spanning most of Chile and including all types of vegetation. It showed that, along the central region (30-38 • S), a large negative rainfall anomaly (nicknamed Megadrought) had a major impact on the gross plant productivity (as measured by the EVI index) during the 2010-2015 period, disproportionally impacted the semi-arid shrublands of north-central Chile (30-33 • S) and the coastal region south of 30 • S, following the curtailed rainfall during this period. In contrast, positive but patchy anomalies in EVI were recorded between 34 • S and 37 • S, which were associated to forest plantations involving non-native species with a rapid initial growth in leaf area [77]. Thus, the complex interplay among climate-driven variability in air temperature, rainfall, and other variables in the latitudinal gradient of Chile, in addition to the different vegetation types (Table 1), should result in marked differences in phenological activity (evidenced in the power spectra for different parks; Figures 2 and A1) in our study. Therefore, it is essential to couple the impact of climate variability with the patterns of air temperature, precipitations, and hydrological regimes to further inspect the contrasting impacts of climate events on the plant productivity and phenology along western South America.
Remote sensing vegetation indices can be used as a robust tool to assist decisionmaking for the conservation and management of areas with high vegetation endemism. Time-lagged responses provide us with insights into how earlier conditions impact vegetation activity. In most open shrubs, at middle and low latitudes, climatic forcing reveals their significant effect on vegetation growth [73]. This study highlights the complexity in biogeographic climate-vegetation biomass effects that should be considered in predictive climate models of phenological variations along western South-America.

Summary
The Pacific coast of western South America has a wide latitudinal range of climate variability, with consequent change in vegetation communities. In recent years there has been increasing interest regarding the effects of climatic oscillations on phenological patterns. Our study demonstrated that the response of EVI to capture the annual cycle of vegetative communities in the temperate climate region minimizes canopy-soil variations and improves sensitivity over dense vegetation conditions. On the other hand, the MEI and SAM indices showed an effect on multiple lags along the latitudinal gradient. Coincidentally, PDO also presents an effect (4-6 month lags) along the study area. Future studies should carefully analyze the PDO as the ENSO/SAM variability only partially determines the phenological patterns along our broad study region. To this end, our results provide a first step towards the design of conservation strategies that integrate climatic information into the managing and monitoring of essential environmental variables, chiefly the phenology terrestrial productivity, along the Pacific coast of western South America.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A
In this contribution, we also showed the NDVI (Normalized Difference Vegetation Index)'s temporal trends ( Figure A1). As EVI, the NDVI also serves for detecting temporal changes in vegetation phenological activity. EVI have a greater dynamic range than NDVI to capture the phenological cycle, due to the fact that NDVI is mainly saturated in the red band region of the electromagnetic region where the energy is strongly absorbed by pigments [78,79] and also these values strongly tend to decrease in the presence of high background amounts of dry biomass. Moreover, NDVI varies depending the chlorophyllthe concentration and amount of vegetation biomass saturating when both chlorophyll and biomass present high levels [78,80]. For example, in the Chiloé Island Ecosystem, recognized as a biodiversity "hotspot", NDVI is susceptible to uncertainty of capturing the annual photosynthetic cycle [79]. EVI perform well under high aerosol loads (which may add significant noise to the signal), biomass burning conditions and topographic effect [78]. However, the temporal co-evolution of EVI and NDVI showed in general complementary patterns across the study region. In terms of the effects of climatic variability over both EVI and NDVI, our results showed the asymmetrical response of the phenological signal.