Spatial Patterns of Vegetation Activity Related to ENSO in Northern South America

Interannual variability of vegetation activity (i.e., photosynthesis) is strongly correlated with El Niño Southern Oscillation (ENSO). Globally, a reduction in carbon uptake by terrestrial ecosystems has been observed during the ENSO warm phase (El Niño) and the opposite during the cold phase (La Niña). However, this global perspective obscures the heterogeneous impacts of ENSO at regional scales. Particularly, ENSO has contrasting impacts on climate in northern South America (NSA) depending on the ENSO phase and geographical location, which in turn affect the activity of vegetation. Furthermore, changes of vegetation activity during multiple ENSO events are not well understood yet. In this study, we investigated the spatial and temporal differences in vegetation activity associated with ENSO variability and its three phases (El Niño, La Niña, Neutral) to identify hotspots of ENSO impacts in NSA, a region dominated by rainforest and savannas. To achieve this, we investigated time series of vegetation variables from 2001 to 2014 at moderate spatial resolution (0.0083°). Data were aggregated through dimensionality reduction analysis (i.e., Global Principal Component Analysis). The leading principal component served as a proxy of vegetation activity (VAC). We calculated the cross‐correlation between VAC and the multivariate ENSO index separately for each ENSO phase. Our results show that El Niño phase has a stronger impact on vegetation activity both in intensity and duration than La Niña phase. Moreover, seasonally dry ecoregions were more susceptible to El Niño impacts on vegetation activity. Understanding these differences is key for regional adaptation and differentiated management of ecosystems.

, with important impacts on the productivity of vegetation within this region (Bastos et al., 2018;Köhler et al., 2018;Liu et al., 2017;Luo et al., 2018;Salas et al., 2020;van Schaik et al., 2018).ENSO is a coupled oceanic-atmospheric event in the Pacific ocean that can lead to excesses or deficits of rainfall for different parts of NSA depending on the season and geographical location.ENSO can be considered an excellent example to investigate extreme weather conditions, and can contribute significantly to assess the implications of extreme climate on vegetation at multiple spatial scales.Increasing our regional understanding of ENSO impacts on climate and vegetation variability, is of paramount importance to predict more accurately the impacts of global change on natural ecosystems, and, in consequence, its further implications to the human populations that inhabit this region.
Hydro-climatological studies have shown important changes in precipitation and temperature associated with ENSO phases (i.e., El Niño, La Niña and Neutral) both at regional and subregional scales (Bolaños et al., 2021;Cai et al., 2020;Poveda et al., 2001Poveda et al., , 2006;;Waylen & Poveda, 2002).These studies have shown that each ENSO phase does not manifest homogeneously across the entire NSA domain, but rather at different intensities and with contrasting patterns of precipitation.The complex topography imposed by the Andean mountains as well as the diversity of ecosystem types characteristic of the tropical rain and dry forests, savannas and wetlands of NSA pose big challenges to understanding the impacts of different ENSO phases on climate and vegetation activity.In particular, there is limited knowledge on how variability among ENSO phases affects vegetation activity taking into account the diverse topography and ecosystems in NSA.
In general, the impacts of ENSO on vegetation are studied for individual events (e.g., El Niño 2015-2016) from local to global scales.Local studies have used ground data and permanent plots to assess plant functional traits, and changes in floristic biodiversity in regions affected by ENSO (González-M et al., 2021;Muenchow et al., 2013;Muenchow et al., 2020).Regional and global analyses have assessed changes on vegetation productivity and photosynthesis (Bastos et al., 2018;Liu et al., 2017;Luo et al., 2018;Patra et al., 2017) and teleconnections between gross primary productivity (GPP) and ENSO (Le et al., 2021).Specifically, for tropical South America, research has focused extensively on productivity in the Amazon biome.For example, van Schaik et al. (2018) compared GPP anomalies among the Amazon sub-basins for El Niño 2015Niño -2016 and found that the northeastern sub-basin hosted the largest GPP loss (56% from October to March).In contrast to analyses based on El Niño events, analyses related to La Niña are scarce and usually carried out at global scales.For instance, Bastos et al. (2013) found that La Niña of 2011 explained more than 40% of the variance in global net primary productivity.For the same La Niña event, Pandey et al. (2017) reported an increase in methane emissions by wetlands (5%).Analyses of individual events contribute to explain yearly changes in global vegetation activity but they are limited in: (a) detecting systematic/concurrent hotspots of ENSO across time; and, (b) identifying differences in ecosystem responses due to local conditions.An alternative to evaluating individual ENSO events is to implement time series analyses of climatological and ecosystem-productivity variables that cover multiple ENSO events.This approach can be very powerful in assessing and comparing changes over time that could be associated with lagged effects during transitions between ENSO phases.Currently, multiple proxies of vegetation activity are available (e.g., greenness, productivity).Some of these variables are either calculated from direct satellite retrievals such as vegetation indices (Huete et al., 1997(Huete et al., , 2002;;Tucker & Sellers, 1986) or data products obtained from model-data fusion methods or data-driven models (Jiang & Ryu, 2016;Jung et al., 2020;Running et al., 2015).Nevertheless, tropical regions such as NSA usually face challenges related to data availability and quality due to cloud cover and lack of ground-truthing (Estupinan-Suarez et al., 2017;Hilker et al., 2012).
Another important challenge in analyzing vegetation activity metrics is that multiple data products are available that may provide contrasting results concerning the relationships between climatic and ecosystem variables.It is often difficult to determine whether one single vegetation variable is preferable over others, but data analysis methods such as dimensionality reduction contributes to capturing a clearer vegetation signal and to reducing noise (Estupinan-Suarez et al., 2021;Kraemer et al., 2020).In particular, these methods can combine the information from multiple data streams into a small set of indicator variables retaining most of the shared information in the data, which offers an opportunity to aggregate different metrics of vegetation activity without a major loss of information.
In this study, we investigated the spatial and temporal differences in vegetation activity associated with ENSO phases (i.e., El Niño, La Niña) in NSA from 2001 to 2014.In particular, we address three main research questions; The study area is encompassed from latitude 14°N to 14°S and longitude 83°W to 60°W, and it is in good agreement with the Northwestern South America region suggested by the IPCC AR6 report (Arias et al., 2021).The highest elevation is 6,500 m asl in the Peruvian Andes, and the wettest region is the Biographic Choco on the Pacific coast with more than 11,000 mm mean annual precipitation.NSA is characterized by regional water cycling and nutrient feedbacks between the Amazon and the Andes (Poveda et al., 2006) as well as atmospheric-vegetation feedbacks (Zemp et al., 2017).It offers a unique opportunity to analyze different topographic and moisture gradients, represented in diverse climatic zones and ecoregions within the tropics.
In this region, the impacts of ENSO on weather are contrasting depending on the season (i.e., austral summer and winter), ENSO phase, and geographical location (NOAA, 2021a).The main climatic impacts in the Caribbean and Pacific coast, and the Amazon are summarized on Table A2.

Vegetation Variables
As a proxy of vegetation activity, we selected two vegetation indices and two data sets derived from models.These variables are the normalized difference vegetation index (NDVI) (Tucker & Sellers, 1986), the enhanced vegetation index (EVI) (Huete et al., 1997(Huete et al., , 2002) ) and the fraction of absorbed photosynthetically active radiation (FPAR) (Knyazikhin et al., 1998;Sellers et al., 1997) from the Moderate Resolution Imaging Spectroradiometer (MODIS) Terra.Specifically, we worked with version 6 of MOD13A2 (layer 1 and 2 for NDVI and EVI respectively) and MOD15A2H (layer 1 for FPAR).GPP was obtained from the Breathing Earth System Simulator (Ryu et al., 2011).Data sets were acquired from the Regional Earth System Data Lab (RegESDL) where they are available in a harmonized grid at 0.0083° spatial resolution and 8-daily temporal resolution (Estupinan-Suarez et al., 2021).Quality flags were previously implemented by the RegESDL for the MODIS products (for more details see Table S2 in Estupinan-Suarez et al. (2021)).Data spans from 2001 to 2014, a period when all four variables overlapped.Detailed information on how to access the RegESDL is in the Supporting Information S1.

Global Principal Component Analysis
Each variable was normalized pixel-wise to a global mean of zero and a variance of one.To extract the main features of vegetation activity, we used dimensionality reduction analysis.Specifically, we used an online version of the principal component analysis (PCA) developed to deal with very large data sets and spherical distortions (Kraemer et al., 2020), and refer to it as Global PCA.Within this framework, the entire data set was projected into a unified principal component space regardless of the time and space features.The Global PCA steps are: (a) to create covariance matrices for each time step; (b) to combine the covariance matrices of each time step into a global covariance matrix; (c) to calculate the loadings of the PCA from the global covariance matrix; and (d) to project data pixel by pixel to the new principal component (PC) space.
In addition, the first PC (PC 1 ) trajectory was cross-checked with the one from the input variables.This is because the PCA orthogonal rotation is arbitrary and not sensitive to the variables directions.In our case, the PC 1 trajectory was opposite to the one from the input variables, then PC 1 was flipped with a simple multiplication by −1.Moreover, the Global PCA handles no data values, and also it assigns a weighted factor to each covariance matrix per time step.Therefore data was not gap-filled.To account for pixel size variations and accurately estimate the final covariance matrix we used the WeightedOnlineStats.jl package (https://doi.org/10.5281/zenodo.6494412)that efficiently handles large data sets.
Furthermore, we conducted a separate Global PCA for each of the dominant natural land cover classes, and referred to it as the land cover PCA for simplicity (see Appendix).With this approach, we compared the fraction of variance captured for PCs when doing a separate analysis for broadleaved evergreen forests, shrublands and grasslands.
Subsequent analyses were performed with the PC 1 that captures the largest fraction of variation, and holds the main information of vegetation activity (see Section 3.1).Henceforth, we refer to PC 1 as the leading vegetation activity component (VAC).

Time-Lagged Correlations Between ENSO Phases and the Vegetation Activity Component
We calculated the Spearman correlation (ρ) between VAC and the bi-monthly Multivariate ENSO index v2.0(MEI) on a pixel-by-pixel level.MEI is calculated based on five atmospheric and oceanic variables (i.e., sea level pressure, sea surface temperature, surface zonal winds, surface meridional winds, and outgoing longwave radiation) between 30°S-30°N and 100°E-70°W over the tropical Pacific, and it takes into account seasonal and sub-seasonal variability (NOAA, 2021b).Time series of monthly MEI are available from the National Oceanic and Atmospheric Administration (NOAA) (https://psl.noaa.gov/enso/mei/,last access 1 September 2021).To match the variable's temporal resolution with the one of MEI, the index was interpolated to 8-daily time steps using B-splines.Based on MEI thresholds by NOAA, the index values were classified as follows: ≤−0.5 correspond to La Niña whereas values ≥0.5 to El Niño, values from −0.5 to 0.5 are Neutral (NOAA, 2021b).
The Spearman correlation was computed separately for El Niño, La Niña and the Neutral phases.Furthermore, we calculated cross-correlations between MEI and VAC to assess time-lagged correlations from zero to 6 months.Due to seasonal variability, we only considered the months that have data from all three ENSO phases (January-line to compare VAC during El Niño and La Niña.Thus, we subtracted the correlation coefficient difference between El Niño (La Niña) and the Neutral phase pixel-wise.Henceforth, we refer to this difference as the correlation anomalies, and interpreted it as a measure of the effect of the ENSO phases on variation in vegetation activity.

Hotspots of Correlations Between ENSO and Vegetation Activity by Ecoregions
To investigate regions with the strongest association between ENSO phases and vegetation activity, we evaluated the study area by terrestrial ecoregions (Figure 1).The ecoregion map integrates biogeographical principles with spatial distribution models of species and communities (Olson et al., 2001).Within these ecoregions, we computed the median Spearman correlation between ENSO phases and the vegetation component.

Global PCA
We calculated a unified principal component space for the entire study area.VAC, the leading vegetation activity component, captured a fraction of 0.45 of the variance.The second and third PC explained a fraction of 0.21 and 0.19 of total variance, respectively.Vegetation indices such as EVI and NDVI were the main variables contributing to the leading component (Table 1).In general, VAC had a good agreement with the input variables in grasslands and shrublands, with correlation coefficients ρ > 0.7 (Figure A1).Nevertheless, the agreement between VAC and GPP decreases when GPP ≥6 gCm 2 /day, highlighting that VAC is not sensitive when productivity is high.Regarding broadleaf evergreen forests, the correlation between VAC and the input variables dropped dramatically ρ ≤ 0.5, except for EVI.In fact, EVI has a nearly linear relationship with VAC within each of the assessed land cover classes, which contrasts with the observed pattern for NDVI.In this sense, we found that the correlation with NDVI highly varies among land cover classes with a ρ ∼ 0.3 in broadleaved evergreen forests in comparison with ρ ∼ 0.7 in savannas.
When conducting the Global PCA separately by land cover class (i.e., land cover PCA), we found that VAC accounted for 0.68 and 0.71 of the variance in shrublands and savannas, respectively, and decreased to 0.33 in broadleaved evergreen forests (Table A3).Additionally, EVI, NDVI and GPP had similar loadings for savannas and shrublands (from −0.51 to −0.53), whereas EVI was the main variable in broadleaved evergreen forest (−0.67 vs. ≥−0.42)(Table A4).

Vegetation Activity During Different ENSO Phases Along Lags
We compared the variability of VAC-MEI correlations among ENSO phases.The correlation anomalies map (Figure 2) shows the differences of the Spearman correlation coefficients between the El Niño (La Niña) and the Neutral phase.The differences of correlation coefficients at each pixel and for each time lag illustrate how the strength of these correlations change spatially and temporally.Specifically, areas near the Pacific coast of Ecuador and Peru, the Caribbean coast, and northern inland savannas (i.e., Llanos and surroundings) are the areas with the highest variability in vegetation activity associated to El Niño and La Niña.
In general, we observed contrary patterns between areas north and south from the equator during all ENSO phases.In the northern region, correlation anomalies were negative during El Niño while positive during La Niña.Additionally, the correlation anomalies were stronger during El Niño than La Niña either in intensity (ρ) and duration (along lags).For example, in the Caribbean coast of Venezuela, the correlation anomalies were  stronger for all lags during El Niño and only from 4 to 6 months during La Niña.In the Sinu valley, the correlation anomalies was stronger during El Niño in the first lags (0-2 months), whereas for La Niña occurred from 3 to 6 months lags.
The opposite was reported for the southern regions near the Ecuadorian and northern Peruvian coast, it was an increase on the correlation anomalies during El Niño that also became stronger along time (2-6 months lag).
During La Niña, the correlation anomalies turned negative and stronger from 3 to 6 months lag.The Amazon biome included in the study area, that is, the upper western Amazon basin, showed a more homogeneous pattern during La Niña than during El Niño.Mostly, during La Niña the correlation diminished from 2 to 6-month lag.
The strongest period was at 3-month lag, then it slightly increased in the following months.Contrary, there was not a unified response of vegetation activity during El Niño.Around 10°S, the correlations were negative whereas in the northern part of the Amazon they were positive (Figure 2a).This heterogeneity is more discernible between 2 and 6 months lag, and highlights the variety of responses within the Amazon biome.
Despite these strong spatial differences, these results consistently show that time lags are more pronounced during El Niño than during La Niña phase.In brief, variability of the vegetation activity is higher during El Niño in the mountain ranges near the southern Pacific coast, Caribbean coast and the Llanos along all lags.On the contrary, changes in vegetation activity related to La Niña are narrowed to specific lags.

Hotspots of ENSO Impacts on Vegetation Activity
We also calculated the lagged Spearman correlation between VAC and MEI by ecoregions for each ENSO phase.We found that El Niño had the strongest correlation with VAC, and followed the same trajectory that the Neutral phase which had weaker values (Figure 3).On the contrary, the trajectory observed during La Niña is opposite to the ones from other ENSO phases, and also had weak correlations.
Contrasting correlation patterns emerged from the ecoregions along the latitude gradient.In fact, we clearly observed that in the northern region, the VAC-MEI correlations were stronger and negative during El Niño while they were weakly positive during La Niña.Conversely, we found a positive (negative) correlation with El Niño (La Niña) in the southern region.Interestingly, El Niño showed the strongest correlation in ecoregions located north and south from the equator, despite the opposite signs of the correlation that is, a positive correlation in the southern region and a negative one in the north.
In particular for El Niño, the strongest correlation, and smallest standard deviation, was found in ecoregions closer to the coasts.Specifically, ecoregions closer to the Pacific ocean (i.e., Central Andean Puna and Central Andean wet Puna) and the Caribbean coast (i.e., La costa xeric shrublands) have the strongest correlation |ρ| > 0.6.Continental ecoregions such as Guianan savanna, Llanos and Beni savanna had |ρ| ∼ 0.5 and a larger standard deviation.Overall, we observed a higher variability of VAC associated to ENSO in dry ecosystems; from the 16 ecoregions with the strongest correlation (|ρ| > 0.4), 69% are either savannas, dry forests, or xeric shrubs/scrubs.Figure A2 illustrates the gradient between ρ and precipitation of the driest quarter, this suggests that the correlation becomes weaker in wetter ecoregions.
To understand the climatic drivers of these ecoregional changes on vegetation, we explored time series of precipitation and air temperature.In general, plots of MEI against the climate variables showed similar trajectories that the ones observed with VAC (Figure A3).As expected, rainfall increases during La Niña in the norhthern ecoregions and decreases during El Niño in the southern ecoregions, highlighting the importance of regional heterogeneity.Otherwise, the trajectories regarding the MEI-temperature correlation were not uniform along lags and often overlapped (Figure A4).

Discussion
We analyzed time series of four vegetation variables to assess the effects of inter ENSO variability across time and space in NSA.To do this, we implemented the Global PCA approach (Kraemer et al., 2020) with the aim to capture a clearer vegetation activity signal that incorporated information from different sources.The derived leading component (i.e., VAC) was correlated to MEI within the different ENSO phases.Our findings show: (a) VAC captures the largest amount of variation from vegetation variables but its performance varies among land cover classes.(b) There are opposite patterns of vegetation activity along time lags during El Niño and La Niña.In the following we discuss these findings in detail.

The Vegetation Activity Component From the Global PCA
Our understanding of vegetation activity in the tropics is still limited.Recent studies are clarifying the climatic drivers at seasonal scales in different ecosystems (Chen et al., 2021;Hashimoto et al., 2021;Li et al., 2021;Uribe et al., 2021), but its variability at longer temporal scales is still in a very early phase.One of the major challenges is that clear-sky observations are limited due to high cloud cover, and occasionally aerosols from fire (Hilker et al., 2012).In addition, ground data for models calibration is limited in the tropics reducing their regional predictability (Jung et al., 2020).Using the Global PCA approach, we captured the main characteristics of vegetation variables related to greenness and productivity.The derived vegetation activity proxy defined as VAC captured the largest amount of variation in seasonally dry climate ecosystems such as savannas and grasslands, but it was limited in the rainforest regions (Figure A1).This was expected due mostly to two main reasons.First, areas such as the western Amazon or the Biogeographic Chocó, where annual precipitation could reach levels above 2,500 mm and 11,000 mm respectively (Poveda & Mesa, 2000), present large seasonal data gaps (Hilker et al., 2012) offering limited yearly data.Second, there are known drawbacks of optical satellites to detect changes in vegetation activity in ecosystem with very dense canopies such as the rainforest (Asner & Alencar, 2010;Huete et al., 1997;Köhler et al., 2018).
In general, our results show that VAC is dominated by the EVI signal.As a consequence, the information content in VAC is closer to vegetation greenness than to productivity.When looking at the dominant land cover classes, we found certain characteristics.For example, all variables show similar linear trends and correlations in grasslands with the exception of a nonlinear relationship between VAC and GPP (Figure A1).This clearly delineates the limitations of VAC regarding its ability to represent spatio-temporal variation of photosynthetic activity.With respect to broadleaved evergreen forests, EVI is the variable with the highest Spearman correlation (ρ = 0.67), and NDVI has the weakest correlation (ρ = 0.3).This was unexpected considering that both are measuring greenness, but reinforces the idea that EVI is more suitable for broadleaved evergreen vegetation because its lower sensitivity to atmospheric effects as well as less canopy saturation (Huete et al., 1997(Huete et al., , 2002)).

Lagged Effects of ENSO in Vegetation Activity in NSA
Regions near the coast of southern Ecuador and northern Peru, including La Sierra, experienced a progressive increase on the correlations anomalies from lag zero to six during El Niño (Figure 2).In this area, the dominant land cover is grasslands and shrublands that might benefit from an increase in rainfall associated to El Niño.A similar response was observed in the Caribbean coast and neighbor inland areas such as the Sinu Valley and Orinoquia savannas (Llanos), but for La Niña phase and from lag three to six.Conversely, these ecoregions showed a negative correlation between VAC and El Niño across all lags, which can be related to a stronger water deficit during the dry season.
In comparison to studies that pointed to a reduction of vegetation activity during El Niño episodes in the Amazon basin (Hilker et al., 2014;Liu et al., 2017;Luo et al., 2018;Patra et al., 2017), we observed a more heterogeneous pattern.During the first 3 months, the western Amazon basin, covered by our study area, had a boost on the correlation anomalies.This could be explained by the fact that solar radiation is a limiting factor in tropical rainforests (Graham et al., 2003;Nemani et al., 2003), which benefit from less clouds during El Niño (Moura et al., 2019).Furthermore, this area is part of the wettest sub-basin in the Amazon (van Schaik et al., 2018), therefore plants are not exposed to water stress in the early stages.However, at a 6 month lag, the correlation of vegetation activity with MEI became predominantly negative.This reveals contrasting vegetation responses related to a lagged effect, and it could be indicative of a link between prolonged solar radiation and the crossing of other environmental thresholds in the Amazon (Brando et al., 2010).
In general, we found negative correlations in the Amazon region during both ENSO phases at lag six, although the spread of these correlations cover a larger area during La Niña phase.Taking into account that higher precipitation and less solar radiation is expected during La Niña, the vegetation activity should drop (Moura et al., 2019).This is consistent with observations by Graham et al. (2003) at the site level in the Panama rainforest, who reported light as the main limiting factor for forest carbon uptake.Nevertheless, this result has to be carefully interpreted when working with optical satellite data due to quality issues in tropical forests and saturation of vegetation indices.
Overall, we summarize the interplay between ENSO phases and VAC as follows.There is a clear and contrasting response between the Ecuadorian and Peruvian regions near the coast on the one hand, and the Caribbean coast and Orinoquia savannas on the other hand, during both the El Niño and La Niña phases.This highlights the importance of spatial heterogeneity within the region, with a clear separation of the northern and southern regions of the study area.This variability is associated to the complex climatology of the region.For instance, Salas et al. (2020) showed that anomalies in precipitation and streamflow during ENSO phases are interlinked with moisture advection of regional jets affecting the anomalies in magnitude and spatial distribution.Thus, the Caribbean Jet, Choco Jet, and Orinoco Jet caused different anomaly patterns in the Caribbean, northern Pacific coast (Choco) and the Amazon, respectively.
Finally, the VAC-MEI correlation anomalies are stronger during El Niño than La Niña both in terms of intensity and delayed effects in the southern Pacific and Caribbean coasts as well as in the Orinoquia savannas.Elucidating these differences in ENSO phases is key for assessing changes, not only for the natural ecosystems but also in terms of further impacts for the population depending on them.

ENSO Hotspots by Ecoregions
Of significant importance is our finding of a high correlation between vegetation activity and MEI during El Niño in dry ecoregions.About 67% of ecoregions with ρ > |±0.4| are considered arid or semi-arid ecosystems.Currently, there is limited understanding on the physiological and climatological mechanisms that may lead to this response in such ecosystems, but some authors have highlighted some implications.For example, González-M et al. ( 2021) reported a negative net biomass balance in tropical dry forests during El Niño of 2015 for NSA.For similarly dry ecosystems such as the Cerrado in Brazil, Zanella De Arruda et al. ( 2016) found that this grassland dominated ecosystem is a source of CO 2 during the drier years.In general, research on tropical ecosystems with seasonal dry climate is still limited (Pennington et al., 2018), and even scarcer regarding their response to ENSO.Taking into account that semi-arid ecosystems are considered the main driver of interannual variability of GPP (Ahlström et al., 2015), and that ENSO is the main climatic driver at the same temporal scale, it is crucial to increase our understanding in this context.Moreover, research focused on tropical savannas and dry forests is extremely relevant in light of decreased rainfall scenarios in places such as the Amazon due to the amplified effects of deforestation, drought, and climate change (Hilker et al., 2014;Parsons et al., 2018;Shiogama et al., 2011;Vilanova et al., 2021;Zemp et al., 2017).
In addition, our results also show that ecoregions with the highest correlation anomalies (ρ > |±0.6|) were in relative proximity to the coast, that is, La Costa xeric shrublands and Central Andean puna, followed by continental ecoregions such as Guianan savanna, Llanos and Benni savanna.This suggests that ecoregions near the coast are slightly more sensitive to MEI than the continental ones independent of their latitudinal location.Another geographical remark is the clear inverse response of VAC in ecoregion's located North and South from the equator during the same ENSO phase.In summary, while the sign correlation between VAC and MEI is positive in the North, it becomes negative in the South and conversely.This pattern can be explained by the opposite effects that ENSO phases have in the different geographical locations in NSA, and it is aligned with the diverse changes on temperature and rainfall previously described for the region (NOAA, 2021a; Salas et al., 2020).Overall, our results suggest that due to the warm ENSO phase, the identified ecoregions have a higher sensitivity to changes in climate.Therefore, we considered them as the hotspots of ENSO effects on vegetation activity.
Lastly, we highlight the opposite trajectory of the lagged correlations during La Niña in comparison to the Neutral phase (Figure 3).Although the correlations between VAC and MEI are weaker during La Niña (from −0.23 to 0.09), the trend consistently has the opposite direction than the Neutral phase along lags.We repeated the same lagged correlation analysis but replacing VAC by precipitation (Figure A3).Our results showed the same trajectories along lags for each ENSO phase and are consistent with the trajectories observed with VAC.These results were unexpected because La Niña is defined as an enhanced Neutral phase from the atmospheric point of view (Bureau of Meteorology, 2012;Philander, 2018).These findings require further investigation to assess the local climate conditions during La Niña in the identified ecoregions.Currently, the assessment of land vegetation variability during Niña and the Neutral phase is limited.But studies such as Pandey et al. (2017) reporting a global methane increase (6-9 Tg CH 4 yr −1 ) during La Niña 2011 reveal the large impact and relevance of the ENSO cold phase.Regionally, important questions emerged after the extensive flooding by La Niña 2011.For example,; how did the ecosystems recover after nutrients were recharged by the river sediments?For how long were soils water-logged?And, What are the consequences for the carbon cycle?Unfortunately, assessments of such impacts of La Niña are still pending.

Limitations and Opportunities
The Global PCA approach aggregates the largest vegetation variability from the entire region into a unified PCA space.This has the advantage that all pixels are projected using the same loadings, and prevents opposite loading signs among neighboring pixels in a pixel-wise analysis.Nevertheless, PCA relies on orthogonal linear transformations and is limited for analyzing nonlinear processes.There is an opportunity to use nonlinear dimensionality reduction methods such as Isomap (Tenenbaum et al., 2000) or locally linear embeddings (Roweis & Saul, 2000) in future studies.Until now, these methods have not been computationally implemented for capturing simultaneously the variability in multiple dimensions (e.g., variables, space, time).
This study was carried out in a limited temporal and spatial domain that was defined by the overlap in the time spans of the assessed variables, together with consistency requirements for their spatial resolution.Nevertheless, there are opportunities to expand this analysis to cover larger areas such as the entire tropical domain and include longer time windows as more information becomes available.Moreover, our method can be used to assess different ENSO indices in addition to the MEI used here.Preliminary studies assessing covariability between ENSO and vegetation variables for watersheds in NSA found similar results regardless of the ENSO index (Estupinan-Suarez et al., 2020), but this needs to be investigated further, and may not be the case for other regions.Otherwise, big hopes are in forthcoming data sets such as the Advanced Baseline Imager (ABI), a geostationary satellite for NSA that will increase data availability of vegetation greenness (Hashimoto et al., 2021).
Better estimates of the vegetation activity will contribute significantly to investigating the lagged correlations shifts in the Amazon, and also causal relations in tandem with other variables such as soil moisture, vegetation type, among others.
Finally, in the last decade, we have witnessed an increase in extreme events.Droughts, floods and compound events are occurring with greater frequency and intensity (AghaKouchak et al., 2020;Beniston et al., 2007;Myhre et al., 2019;Papalexiou & Montanari, 2019;Rahmstorf & Coumou, 2011).Likewise, extreme weather events have also occurred in the past during ENSO years.Government strategies and environmental interventions (e.g., policies) to address the climate and ecological crisis could benefit from a deeper regional understanding of ENSO-induced extreme weather events and associated changes in vegetation activity.This can be seen from two perspectives.First, due to the expected increase of ENSO events in terms of frequency and intensity and its causal effect in the tropics (Cai et al., 2014;Fasullo et al., 2018;Yun et al., 2021), it is necessary to pay more attention to tropical regions taking into account the multiple facets of ENSO.Second, although predictions of El Niño and La Niña teleconnections remain largely uncertain for the extratropics (Johnson et al., 2022;Singh et al., 2022;Williams et al., 2023), identifying hotspots of high vegetation variability during multiple ENSO events could elucidate vegetation characteristics related to repetitive exposure to climate extremes in the past, and may inform what to expect during similar conditions in the future.

Conclusion
Our results show that in NSA, a region with high heterogeneity in terms of climate and biophysical characteristics, vegetation activity presents opposite patterns within the same ENSO phase.Thus, when the correlation between vegetation activity and ENSO is positive during El Niño south of the equator, it is negative toward the north, and vice versa during La Niña.Understanding these differences is key for regional adaptation to more frequent and prolonged droughts or floods.Moreover, this could also contribute to the development of differentiated management plans for different ecosystem types that respond differently to the extreme phases of ENSO.
According to the lagged effects of ENSO, the variability of vegetation activity is stronger in intensity and more prolonged during El Niño than during La Niña.Interestingly, this is independent of the sign of the correlation, highlighting that positive and negative correlations are stronger during El Niño events.In addition, seasonally dry ecosystems require attention, considering that they are the most sensitive to ENSO.Advancing knowledge of these ecosystems is key to understanding the interannual variability of vegetation activity at regional and global scales.It is also important to promote research efforts in this ecosystem, since in NSA they are mostly focused on the tropical rainforest.(See Appendix Figure A4)

Global Research Collaboration
This research is the result of an ongoing collaboration between international scientists that began at the PEACE meeting in Colombia in 2016 (www.peacecolombia.org).A large number of the participants were scientists from Colombia, as well as national scientists working abroad; we thank them all.In addition, our study benefited from the RegESDL, a data infrastructure that facilitates big data analysis in the region, launched in 2021 and also an initiative that started during the PEACE meeting.In addition, the contribution and in-person participation of GP was made possible thanks to the DAAD grant

Figure 1 .
Figure 1.Map of northern South America.Gray polygons are selected ecoregions from Olson et al. (2001).The numbers indicate the ecoregions name.The black lines are the country borders.

Figure 2 .
Figure 2. Map of correlation anomalies.The color map represents the difference in lagged Spearman correlations between MEI and the leading vegetation activity component (VAC) for (a) the El Niño and Neutral phase, and (b) the La Niña and Neutral phase.The black lines are country borders.

Figure 3 .
Figure 3. Median values of the lagged Spearman correlation between VAC and MEI by ecoregions.Ecoregions are in rows.From left to right: first column shows the ecoregion location.Second column represents the median correlations within the ecoregion by ENSO phases.The ribbon correspond to the standard deviation.Data distribution from El Niño, La Niña and Neutral phases by lags is showed in columns three to five respectively.Showed ecoregions have the highest correlation and an area larger than 50,000 km 2 .

Figure A4 .
Figure A4.Median values of the lagged Spearman correlation MEI and air temperature at 2 m above the ground (source: source: ERA-Interim 80 km) by ecoregions.Ecoregions are in rows.From left to right: first column shows the ecoregion location.Second column is the median correlation values within the ecoregion by ENSO phases.The ribbon correspond to the standard deviation.Data distribution from El Niño, La Niña and Neutral phases by lags is shown in columns three to five respectively.Shown ecoregions have the highest correlation and an area larger than 50,000 km 2 .

Table 1
Loadings From the Global PCA by (Normalized) Variable

Table A1 Area
Calculated per Land Cover Classes for NSA in 2014 Based on the Global Land Cover Map by ESA (Source: ESA (2017)) Note.*: Anomalies not reported.Source: NOAA (2021a).

Table A3 Explained
Fraction of Variance by Principal Components (PC) From the Global PCA and the Land Cover PCA (See Appendix Table A4)

Table A4 First
Principal Component Loadings From the Global PCA and the Land Cover PCA Research Stays for University Academics and Scientists 2019 (57440915).LMES acknowledges the support of DAAD and its Graduate School Scholarship Program (57395813), as well as the International Max Planck Research School for Global Biogeochemical Cycles.All authors acknowledge support from Ulrich Weber and Fabian Gans for data management and preprocessing.The Regional Earth System Data Lab was supported by ESA via the ESDL (now DeepESDL) project.Authors thank the reviewers for their valuable comments and feedback.The Max Planck Society paid for the processing costs and open access to this publication.Open Access funding enabled and organized by Projekt DEAL.