Monitoring and predictive mapping of floristic biodiversity along a climatic gradient in ENSO’s terrestrial core region, NW Peru

1 –––––––––––––––––––––––––––––––––––––––– © 2020 The Authors. Ecography published by John Wiley & Sons Ltd on behalf of Nordic Society Oikos This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. Subject Editor: Bethany Bradley Editor-in-Chief: Miguel Araújo Accepted 6 August 2020 43: 1–13, 2020 doi: 10.1111/ecog.05091 43 1–13


Introduction
Studying the relationship between environmental gradients and species turnover has become a cornerstone of biogeography and community ecology since Humbold's and Bonpland's travels to South America more than 200 years ago (Whittaker 1967, MacArthur 1984. Naturally, water availability drives floristic composition and biomass production in (semi-)arid regions such as NW Peru (Muenchow et al. 2013a, c), which is the area under study in this manuscript. In fact, water availability is probably the most important limiting abiotic factor in NW Peru's tropical dry forests, where it controls plant recruitment, establishment and survival, seedling production, vegetation cover as well as species diversity (Balvanera et al. 2011, Espinosa et al. 2011, Salazar et al. 2020). However, water availability does not only vary in space but also in time, especially in regions heavily affected by the El Niño Southern Oscillation (ENSO).
ENSO is a recurrent climate phenomenon causing anomalies around the globe (Capotondi et al. 2015). ENSO consists of alternating episodes. El Niño represents the warm and La Niño the cold episode of this cycle. ENSO neutral refers to periods when neither El Niño nor La Niña episodes are present. The study area is part of the core region where the ENSO phenomenon exerts its greatest terrestrial influence (Richter and Ise 2005). Generally, El Niño episodes bring more rain than usual and La Niña episodes enforce the already dry conditions in the study area (< 50 mm of annual precipitation along the coast of NW Peru; Fig. 1d; Rollenbeck et al. 2015).
Especially El Niño episodes provoke drastic ecological responses (Kogan and Guo 2017) but occur irregularly (five times over the last 20 years in the study area), and are highly variable in terms of intensity. Nevertheless, the ecoclimatic impact of La Niña episodes is comparable in magnitude with that of El Niño episodes. For instance, La Niña episodes were most likely the main driver for the hiatus in global warming at the beginning of the 21st century (Kosaka and Xie 2013). Therefore, one objective of this manuscript is to study the effect of different ENSO episodes, including La Niña episodes, on the floristic composition, comprised of all observed vascular plants, of the same ecosystem. In any case, ENSO drives a large amount of the interannual variability of woody plant growth across the tropics (Rifai et al. 2018) though the specific response to El Niño episodes may vary from region to region (Phillips et al. 2009). Overall, the effect of extreme events is even more pronounced in (semi)-arid ecosystems (Holmgren et al. 2006), and is especially true for ENSO's terrestrial core region in NW Peru and SW Ecuador (Muenchow et al. 2013a, Rollenbeck et al. 2015.
NW Peru is also home to the tropical dry forest (TDF) ecosystems. TDFs are characterized by a high rate of endemism and feature an outstanding biodiversity which is almost as high as that of humid tropical forests (Espinosa et al. 2011, Pennington et al. 2018. At the same time, they are the most threatened tropical ecosystem (Miles et al. 2006), and in general chronically understudied (Escribano-Avila et al. 2017, Muenchow et al. 2018. Despite their unique biodiversity and threatened status, the TDFs of NW Peru are vastly understudied, and the spatio-temporal effect of extreme events like El Niño or La Niña episodes on the biodiversity patterns in this ecosystem remains largely unknown. Statistical learning has a long tradition in biodiversity research for detecting patterns in ecological data and predicting from these patterns (Guisan and Zimmermann 2000). It is especially useful for quantifying species-environment relationships and for predicting the occurrence of species or species richness (predictive mapping), and therefore is an indispensable instrument for biodiversity management (Fremout et al. 2020). Furthermore, predictive mapping can reveal ecological changes in space and time which otherwise might be hard or impossible to detect (Lovelace et al. 2019). Therefore, it is a crucial tool for gaining detailed insights into ecosystem dynamics, which other more descriptive methods such as global biodiversity measures and ordination techniques cannot capture. We aim at spatially predicting, for the first time, how different ENSO episodes might change and reshape the spatial distribution of the floristic composition in the study area.
Species richness is strongly associated with precipitation in the tropics (Esquivel-Muelbert et al. 2017a, b). This is especially true for dry tropical ecosystems where water shortage is the main constraint to plant growth. Still, edaphic and topographic properties are sometimes equally important in structuring vegetation of tropical ecosystems (Soethe et al. 2008, Laurance et al. 2010, Condit et al. 2013). Yet, the corresponding effects on the floristic composition of tropical (semi-)arid vegetation communities remain scarcely investigated (Peña-Claros et al. 2012, Muenchow et al. 2013b, Ulrich et al. 2014. This is especially true in terms of how soil properties influence vegetation communities as soon as water is no longer the limiting factor. For instance, nitrogen may increase primary production even in the presence of droughts and grazing pressure (Whitford andSteinberger 2011, Kinugasa et al. 2012). Hence, this study investigates how much the impact of edaphic and topographic variables varies in explaining the floristic composition between dry and humid ENSO years. Additionally, an irrigation-fertilization experiment should help to quantify in more detail the beneficial effects of the water-nutrient interaction on vegetation diversity and cover.
In summary, in this manuscript we study the effects of different ENSO episodes on biodiversity patterns and floristic composition as well as the effect of environmental variables on the floristic composition under varying rainfall patterns. Specifically, we tested the following three hypotheses: 1. We investigated the change in the different levels of biodiversity along a pronounced climatic gradient across different ENSO episodes. Overall, we expected a boost in biodiversity as soon as water was no longer the predominant limiting factor. We were particularly interested in how different ENSO episodes modify the spatial distribution of the floristic composition. Overall, we expected a more pronounced development of vegetation formations along the climatic gradient with wetter conditions (Engelbrecht et al. 2007). 2. Water availability limits plant cover and species richness along the entire climatic gradient in the (semi-)arid study area. Therefore, we were interested in evaluating the probably changing effects of soil characteristics and topography on floristic composition and vegetation cover (biomass production) with varying rainfall patterns as caused by different ENSO episodes. We expected that the influence of topographic and edaphic variables might gain in influence in wetter years, i.e. as soon as water is no longer the limiting factor. 3. To quantify the beneficial effects of the water-nutrient interaction on vegetation cover development and species richness, we executed an irrigation-fertilization experiment simulating different ENSO rainfall scenarios. We expected that the water input of the rainfall scenario representing a Super-Niño episode in combination with the fertilizer would be most conducive to plant cover and species richness development.

Study area
The study area stretches for 120 km along a pronounced bioclimatic gradient in NW Peru and features a desert along the  Pacific coast in the west and TDFs near the Andean foothills in the east (Fig. 1). The climatic gradient is at the same time a gradient of increasing human influence especially in terms of agricultural activity towards the Andean foothills (Muenchow et al. 2013a, c). Agricultural activity includes cultivation of annual crops and fruit trees and extensive grazing caused by free-ranging livestock. The ENSO phenomenon heavily impacts the study area (Richter and Ise 2005). In general, El Niño episodes are associated with more rain than usual while La Niña episodes amplify the already dry conditions. However, even this pattern is highly variable and can occasionally be reversed (see e.g. the years 2005, 2008 and 2012 in Fig. 1d). Additionally, so-called 'local' or 'coastal El Niño episodes' can take place as observed in 2017. Coastal El Niño episodes refer to the regional warming of coastal waters in the Gulf of Guayaquil in December as opposed to the Pacific basin-wide El Niño episode, which occurs when warm North Australian waters reach the South American coast.
Therefore, to put the ecological results into context, it is important to keep in mind the varying rainfall patterns as produced by different ENSO episodes across the studied years. The observed years differed largely in terms of precipitation. 2011 corresponded to a dry La Niña episode and the whole study area experienced drier conditions than normal. By contrast, 2012 was a humid La Niña episode, and 1.5 times and 3 times more rain fell in Piura and Chulucanas, i.e. in the middle to eastern part of the study area (Fig. 1a). 2016 was a moderate El Niño episode and rainfall was double the median input in Paita and Piura, whereas in Chulucanas the median value was almost reached. 2017 coincided with a coastal El Niño episode and was the only studied year during which the whole study area experienced exceptionally wet conditions: Paita received almost 15 times (379 mm), Piura almost 11 times (780 mm) and Chulucanas 6 times (1706 mm) more rain than on median (Fig. 1d). Moreover, it is noteworthy that vegetation development changes drastically in the study area during wet episodes but remains more or less uniform during neutral years (Richter and Ise 2005, Holmgren et al. 2006, Muenchow et al. 2013c).

Landscape and small-scale experiment vegetation sampling
We sampled 50 vegetation plots of 30 × 30 m 2 size along the climatic gradient between Paita and Chulucanas ( Fig. 1). To ensure that all parts of the study area were equally likely to be sampled, we stratified the random sampling by distance to the Pacific Ocean in ten classes (Supplementary material Appendix 1). Prior to the sampling, we excluded riparian, agricultural and urban areas. Each plot was located with a handheld Garmin GPS device in the field. Field sampling took place towards the end of the vegetation period (March-May) in 2011, 2012, 2016 and 2017. In each plot, we recorded all vascular plant species and determined their cover as portion of the plot area in percentage points. Specifically, plant cover refers to the relative projected area covered by a species (Damgaard 2014), and is a proxy for biomass. Prior to the analysis, we transformed the plant cover values in accordance with the decimal scale by Londo (Londo 1976). The nomenclature follows the conventions of the Tropicos online database (Missouri Botanical Garden 2018).
In addition to the plot survey at the landscape scale, we conducted a small-scale irrigation-fertilization experiment to quantify the beneficial effects of the water-nutrient interaction on species richness, plant cover and biomass production. We established the experiment in the middle of the study area on on a protected dryland forest inside the campus of the University of Piura (5°10′S, 80°38′W, see Piura in Fig. 1a). Hence, the experiment was located in the same ecosystem as the landscape study (see previous paragraph) and more specifically in the middle of the observed precipitation gradient, ranging from 26 mm (Paita; annual median precipitation) to 272 mm (Chulucanas; Fig. 1a). The environmental conditions are typical for the region in terms of climate, topography and soil properties (Muenchow et al. 2013a, b). Inside a sandy area devoid of woody vegetation, we placed 12 experimental plots of 3 × 3 m 2 , which were regularly spaced along a rectangular grid of 1 ha. Plots were assigned to one of three irrigation treatments: The first treatment represented a Super El Niño episode (example of 1997/1998: 1780 mm), the second represented a moderate El Niño episode (example of 1991/1992: 258 mm). The last treatment represented the baseline with no additional water input and corresponded to the year in which the irrigation experiment was conducted (2013) which was a neutral episode (54 mm; see also Fig. 1d). The irrigation water stemmed from rain-collecting tanks. In addition to the water treatment, two randomly chosen plots of each irrigation treatment received a fertilization treatment of 200 kg ha −1 granular nitrogen fertilizer (NH 4 NO 3 ). The experiment was carried out between December 2012 and May 2013. Monitoring was executed two times per month and included the identification of plant species (Missouri Botanical Garden 2018) and corresponding coverage (%). Please refer to Supplementary material Appendix 2 for a detailed description of the irrigation-fertilization experiment.

Environmental variables
We randomly took and mixed three soil samples from 15 cm and 30 cm depth in each plot in 2011. Soil composition should not vary too much between the years in which vegetation sampling took place since the study area largely features arenosols, i.e. a very weekly developed soil type. On the other hand, soil chemistry can vary largely between two years. For instance, El Niño episodes can increase soil respiration, i.e. soil CO 2 emissions, by a factor of 100, especially in the vicinity of trees (Salazar et al. 2020). Subsequently, we measured in the laboratory variables known to be important for vegetation development: pH, electrical conductivity (µS), carbonto-nitrogen ratio, soil texture (%), the concentrations of the cations Ca, Mg, K and Na (cmol kg −1 ), and skeletal content (%; measured as gravimetric proportion of stones > 2 mm; see Muenchow et al. (2013a) for a more detailed description).
We computed the topographic variables aspect, altitude, catchment slope, catchment area, soil wetness index, plan and vertical curvature from a digital elevation model with a 30 × 30 m 2 resolution (LP DAAC 2018). We downloaded the normalized difference vegetation index (NDVI) raster imagery from the moderate resolution imaging spectrometer (MODIS; LP DAAC 2018) available for the vegetation period (March, April) for each studied year (2011,2012,2016,2017) at a 250 × 250 m 2 resolution for our study area. Since three NDVI scenes were available for each vegetation period, we computed the mean NDVI for each vegetation period.
Local precipitation values were collected from automated climate stations (Fig. 1d) operated by the Univ. of Piura near Paita ), in Piura (1991) and in Chulucanas (1997.

Biodiversity and plant cover analysis
We calculated three measures of biodiversity for each observed year, namely alpha, beta and gamma diversity (Whittaker 1967). Alpha diversity refers to the mean number of species per plot. Gamma diversity is the total number of observed species across all plots. Beta diversity, in the sense of species turnover, is expressed as the range of the first detrended correspondence analysis (DCA) axis.
The unit of a DCA axis is given in units of standard deviation (SD). 4 SDs correspond to a complete species turnover while a difference between 1 and 1.4 SD units already correspond to half a change in species composition (Legendre and Legendre 2012). Moreover, the floristic composition was assessed with the help of a DCA (Hill and Gauch 1980) to which the decimal scale transformed cover values of the species-plot matrix were subjected (Londo 1976). Floristic composition implicitly represents the named biodiversity measures, and additionally considers plant cover. DCA tries to condense as much as possible the observed variance of the input matrix in a low-dimensional space (usually 2-3 axes), thereby extracting the main apparent gradients. To test if beta diversity and floristic composition, respectively, varied significantly between years, we analyzed the multivariate homogeneity of group dispersions (Anderson et al. 2006).

Modeling and predictive mapping
One objective of this study was the spatio-temporal mapping of floristic composition represented by the scores of the first DCA axis along the climatic gradient in the study area. Please note that the resulting prediction map not only shows species composition but also differences in beta diversity on the predicted pixel level. This is because the units of a DCA axis are given in SDs (see also previous section). Muenchow et al. (2013c) have already shown that NDVI is an excellent predictor for vegetation formations in the study area. Therefore, we modeled the scores of the first DCA axis representing the main observed floristic gradient in the study area with the help of a generalized additive model (GAM) as a function of the interaction between NDVI and year: where i refers to the ith observation, α is the intercept, β 1 is the estimated coefficient for year and f is a spline smoother for the NDVI:year interaction. Please refer to Supplementary material Appendix 3 for a visualization of the non-linear relationship between the response variable (ordination scores) and NDVI by year.

Scores
The sole purpose of the model is the spatial prediction of the ordination scores, not statistical inference or the interpretation of the models' coefficients. We will use variation partitioning later on to assess how the influence of environmental predictors has changed over the studied years (see Impact of environmental variables). Consequently, here we are only interested in a good predictive performance of the model to make sure that it is able to generalize, i.e. that it is able to predict reasonably well unseen data (James et al. 2013, Kuhn andJohnson 2013). To make sure this is the case, we used spatial cross-validation. Overall, cross-validation accounts for over-optimistic predictions on the training set by randomly splitting a data set into partitions used for training and testing (Brenning 2012, Lovelace et al. 2019. Spatial cross-validation is an extension that measures a model's ability to spatially predict the response variable in the presence of spatial autocorrelation (Schratz et al. 2019). We ran a 100-repeated 5-fold spatial cross-validation in which the partitioning is based on k-means clustering (k = 5; Ruß and Brenning 2010). We also made sure that all annual observations of one plot were jointly placed in either the test or the training data set (Meyer et al. 2018). We used the normalized root mean-square error (NRMSE) in percent as performance measure for the spatial cross-validation. where n is the number of observations, e corresponds to the residuals and y refers to the observed values.

Impact of environmental variables
Variation partitioning quantifies the variation explained of the response, in our case the first two DCA axes, by one predictor group, in our case edaphic variables, while controlling for the effect of another predictor group (Peres-Neto et al. 2006, Borcard et al. 2018, in our case topographic variables. Internally, variation partitioning computes a redundancy analysis (RDA) separately for each predictor group and for all predictor groups together, and subsequently quantifies how much each predictor group explains alone and jointly with the other predictor group of the variation of the response. This is frequently displayed with the help of a Venn diagramm. We ran the variation partitioning for each studied year (2011,2012,2016,2017) to assess how edaphic and topographic variables influence vegetation composition under increasingly humid conditions. Prior to the variation partitioning, we log-transformed electrical conductivity, the CN ratio and the catchment area. Subsequently, we excluded all collinear variables which exceeded a variance inflation factor > 3 leaving us with seven edaphic and eight topographic variables. We refrained from a forward selection of explanatory variables (Blanchet et al. 2008) because (i) we only used a reasonable number of meaningful, non-collinear variables that (ii) were available for all four years thus guaranteeing inter-annual comparability.

Irrigation-fertilization experiment
For presenting the results of our irrigation-fertilization experiment we used descriptive statistics and plotted the increase in the number of species, plant cover as well as water input against time.

Biodiversity and plant cover
Mean alpha diversity was lowest during the dry La Niña (2011) and highest during the moderate El Niño episode (2016; Table 1). Compared to the humid La Niña (2012) and the moderate El Niño (2016) episodes, species richness was lower in the east of the study area during the very humid neutral year (2017), however, higher within a distance of 20-25 km to the coast (left panel of Fig. 2). Species turnover (beta diversity), expressed as the range of the first DCA axis (see also Fig. 3), was also lowest during the dry La Niña episode (2011), slightly increased in 2012 (humid La Niña episode) and 2016 (moderate El Niño episode) and reached its highest value in 2017 (very humid neutral year). Still, there were no significant interannual differences in beta diversity in accordance with the test for multivariate homogeneity of groups dispersions (ANOVA F-test: 0.05069). Gamma diversity increased from 2011 (dry La Niña episode) to 2012 (humid La Niña episode), remained almost the same in 2016 (moderate El Niño episode), and decreased in 2017 (very humid neutral year; Table 1). By contrast, the latter showed by far the highest plant coverage along the entire gradient (right panel of Fig. 2 and see Spatio-temporal mapping).
The first two DCA axes explained 68% of the observed variance with the first axis alone contributing 48%. The first axis was mainly associated with precipitation while the second axis was related to topographic variables such as vertical curvature and catchment area. Overall, the DCA scatter cloud became more compressed over the years, and the plots became increasingly positioned along the first axis (Fig. 3). Put differently, the importance of the second axis steadily decreased from 3.22 during the dry La Niña episode (2011) to 1.35 units of standard deviation during the very humid neutral year (2017). By contrast, the first axis was more important during all episodes which were wetter than the dry La Niña episode (Table 1).

Spatio-temporal mapping
Our model explained 76% of the observed variance with a spline smoother absorbing 3 df (equation 1 and see Modeling and predictive mapping). The spatially cross-validated NRMSE also indicated a satisfying fit of 27% (SD 8.53% over 100 repetitions).
Since the predicted DCA scores represent the floristic composition of a pixel, clusters of pixels sharing similar values (visible as similar color tones in Fig. 4) can be understood as belonging to the same vegetation formation. A clear trend towards more pronounced differences in vegetation formations became apparent between the studied years (Fig. 4). The dry La Niña episode (2011) showed a floristic composition that was almost uniform across great parts of the study area (Panel 2011 of Fig. 4). This changed during the wet La Niña episode (2012) where three distinct vegetation formations developed (Panel 2012 of Fig. 4). In the west, we observed a desert-like formation with a sparse herb and grass cover, which at times was accompanied by isolated trees and bushes (blue pixels in Panel 2012 of Fig. 4). This formation developed into an open xeric shrubland (green/yellow pixels in Fig. 4) where Acacia macracantha and Encelia canescens were the dominating species. Further east, the shrubland turned into a TDF formation (orange/red pixels in Fig. 4) with Cordia lutea, Prosopis pallida, Cenchrus echinatus and Antephora hermaphrodita being the most common species. The differentiation of these three vegetation formations was even more pronounced during the moderate El Niño episode (Panel 2016 of Fig. 4). This trend was both amplified and disrupted during the very humid neutral year (2017) during which a dominant grass-herb formation stretched from the coast until far behind the city of Piura (dark blue pixels in Panel 2017 of Fig. 4). TDF formations were also less dominant (orange to red pixels). Though species composition remained roughly the same during the very humid neutral year (2017) compared to previous years, plant cover

Influence of environmental variables
The influence of edaphic and topographic variables on floristic composition was relatively small in the driest year (2011) with 25 and 28% of explained variance, respectively (Fig. 5).
The explained variance roughly doubled in all other years,   9 however, did not further increase with an excessive supply of water in the most humid year (2017; Fig. 5).

Irrigation-fertilization experiment
The experiment revealed that edaphic variables are conducive to vegetation development. The fertilized plots representing moderate El Niño episodes exhibited a slightly greater plant coverage than in their unfertilized counterparts (middle panel of Fig. 6). Adding even more water revealed the true beneficial effects of fertilizing. The fertilized plots simulating Super El Niño episodes showed a threefold increase in plant cover compared to their unfertilized counterparts (right panel of Fig. 6). Nevertheless, there was an upper limit regarding the beneficial effects of nutrient addition. Plant cover, and thus biomass production, stopped increasing or even decreased in the fertilized Super-El Niño plots after 8 March 2013, which corresponded to 1487 mm of irrigation and 200 kg ha −1 of NH 4 NO 3 . Species richness followed a similar pattern -highest numbers in species were recorded in the fertilized Super-El Niño plots, though with just seven observed species the overall number remained low (Fig. 7).

Biodiversity and biomass
Beta diversity (expressed as the range of the first DCA axis), increased with wetter conditions (Fig. 3, Table 1). However, another measure of beta diversity, the multivariate dispersion, did not change significantly between years (see Biodiversity and plant cover) though the spatial prediction yielded drastic changes in the spatial distribution of the floristic composition and beta diversity (see next section). Concordantly, the change of the shape of the DCA point scatter between the relatively humid years (panels 2012 and 2016 of Fig. 3) and the very humid neutral year (panel 2017 of Fig. 3) can be most likely attributed to the increase of plant cover (see also next paragraph). Overall, higher species turnover was mainly Figure 6. Results of the irrigation-ferilization experiment. The blue line refers to the water input, note that the left y-axis was log-transformed. The solid salmon and the gray dashed lines represent the observed plant coverage in % (right y-axis). driven by short-lived, less abundant species. Cleland et al. (2013) reports similar findings from temperate grasslands. Alpha and gamma diversity increased with higher seasonal precipitation amounts. The relationship, however, was non-linear, i.e. the increase in diversity reached a plateau with the rainfall values recorded during the humid La Niña (2012) and moderate El Niño (2016) episodes. Then even more water input as observed during the very humid neutral year (2017) resulted in a slightly decreased diversity but also in a massive and sudden increase of overall plant cover due to higher biomass production. Interestingly, only a few grass (e.g. Aristida adscensionis, Eragrostis cilianensis) and forb species (e.g. Tiquilia paronychioides) were in large parts responsible for this increase (see Spatio-temporal mapping). Our irrigation-fertilization experiment confirms the findings of the landscape scale, i.e. the increase in plant cover was much more pronounced than the increase in species richness (Fig. 6, 7). One explanation might be that under favorable environmental conditions a limited amount of species become the dominant actors of an ecosystem, whereas under intermediate levels of stress and disturbance higher species richness is supported (Brown 2014). A few dominant species outcompeted other species in 2017 by following a ruderal strategy, i.e. by growing as fast and producing as much offspring as possible in the short period of favorable conditions (Schmidtlein et al. 2012). This is in contrast to other dryland studies which showed that dormant and annual species increased biodiversity under more humid conditions (Gutiérrez et al. 2000, Holmgren et al. 2006. Finally, freeranging livestock might also have an impact on vegetation dynamics in the study area, especially in its eastern parts with more than 100 mm of annual precipitation. Livestock might alter floristic composition and diversity while supporting especially plants adapted to grazing. Aside from the spectacular germination of annuals, perennial woody plants also largely benefit from wetter periods, e.g. through increased flower and fruit production (Holmgren and Scheffer 2001). However, we only recorded plant cover, where the corresponding interannual difference is only marginal compared to previous years. Still, rainy periods are very conducive to the regeneration of woody vegetation in (semi-) arid areas (Bowers 1997, Brown et al. 1997, Sitters et al. 2012. For example, the growth and recruitment rate of Prosopis pallida is almost twice as high during wet El Niño episodes than during neutral or dry episodes in northern Peru (Lopez et al. 2006).
Overall, water surplus converts the (semi-)arid landscape interspersed with shrubs and trees temporarily into a tree savannah. Shallow underground water near the Pacific coast is largely responsible for the occurrence of perennial species in the study area. In fact, dendrochronological studies have shown that the growth of the most frequent tree species in the study area, Prosopis pallida, is nearly independent of the interannual climatic variability because frequently the water input of the low annual precipitation and nearby rivers suffices to replenish the groundwater (Brown and Archer 1990, Throop et al. 2012, Salazar et al. 2018b. The ordination and modeling results of our study provide ample evidence how current-year precipitation affects biodiversity dynamics. Still, it might be possible that to a much lesser extent previous year rainfall might also have an effect in the study area (Tenhumberg et al. 2018). For instance, Dudney et al. (2017) observed an increase in grass and a decrease in forb abundance when rainfall was high one year earlier. Similarly, Walter et al. (2011) suggest that grasses might have a 'drought-memory'.

Floristic gradient mapping across ENSO episodes
The spatial prediction maps highlight that specific vegetation formations become more pronounced with increasing water input. In the driest year (2011), a rather uniform distribution of the same floristic composition (green pixels in Fig. 4) developed along large parts of the study area. The corresponding DCA scores close to 0 indicates that plant community composition was primarily made up of species that are largely independent of the climatic gradient (Hill andGauch 1980, von Wehrden et al. 2009). This is because in 2011 rainfall barely occurred along large parts of the study area, and therefore especially perennial species were recorded which can and must endure drought years (Salazar et al. 2018b). By contrast, during the humid La Niña (2012) and the moderate El Niño episodes (2016) the same three distinct vegetation formations developed along the gradient. This can be explained by similar rainfall patterns though 2012 was a humid La Niña and 2016 a moderate El Niño episode. A large water surplus during the very humid neutral episode (2017) led to the disruption of the strictly ordered formations along the humidity gradient. Dense grass and herb cover dominated large parts of the study area. The increased cover dominance of annual species is visible in larger negative DCA score values. Interestingly, vegetation patches (yellow, green and orange pixels) can be found near the coast which are otherwise more typical of the eastern part of the gradient (see Spatio-temporal mapping and Fig. 4). TDF formations (orange and red pixels) are less visible in 2017 because the dense grass/herb cover temporarily converted the landscape into a savannah-like vegetation formation (Espinosa et al. 2011, Salazar et al. 2018b. Another interpretation could be also that the NDVI values, our main predictor for the floristic composition, is similar for the dense cover of annual plants and perennial woody plants in 2017 (Salazar et al. 2018a, Peña andBrenning 2015).
Overall, our spatial mapping approach concedes detailed insights into how floristic composition and beta diversity distribution patterns change in response to varying rainfall patterns caused by different ENSO episodes. Similarly, Lovelace et al. (2019) (chapter 13) predicted floristic composition on a Peruvian fog oasis which revealed that the vegetation belts along the mountain slope are partly interrupted. Consequently, predictive mapping can contribute to a more nuanced understanding of ecosystem dynamics. Naturally, such insights cannot be captured by global beta diversity measures such as the analysis of multivariate homogeneity of groups dispersions which in our case did not yield significant results (see Biodiversity and plant cover and Biodiversity and biomass).
Finally, the spatial prediction of species richness (Supplementary material Appendix 4) supports the findings of the floristic gradient mapping since it revealed similar spatial diversity patterns and dynamics.

Changing impact of environmental variables with wetter conditions
Variation partitioning clearly reveals that environmental predictors (topography, soil) gain in importance for explaining plant diversity and productivity dynamics as soon as water is no longer restricting plant growth (see Influence of environmental variables and Fig. 5). Interestingly, the shared explained variance in the floristic dataset is roughly the same for the years experiencing wetter conditions (humid La Niña episode (2012), moderate El Niño epidose (2016) and the very humid neutral episode (2017; coastal El Niño) though 2017 was the most humid by far. This indicates that there is an upper limit to the beneficial effects of the water-soil and water-topography interaction. Topographic variables such as curvature and catchment area play some role in structuring the vegetation composition in the study area especially in dry years (Muenchow et al. 2013a, c). Similar findings were observed in Mexico where topography also partly explained the composition of dry forests (Mendez-Toribio et al. 2016). However, the most important topographic variable in the study area was altitude since it is highly correlated with mean annual precipitation. Of course, edaphic variables such as pH are also related to precipitation. The influence of edaphic variables on vegetation composition could even increase with more developed soils. However, large parts of the study area consist of sandy soils (arenosols). Only in the east of the study area, soils show more signs of humification and a brownish color due to iron oxide release from primary minerals.
Nutrients can play a major role for vegetation development in drylands (Ronnenberg andWesche 2011, Whitford andSteinberger 2011) as confirmed by our irrigation-fertilization experiment. In fact, the water-nitrogen interaction can result in a plant coverage three times as high as achieved when adding the same amount of water but without additional nutrients (Fig. 6). This is of utmost importance for sustainable yet productive agricultural management in the study area, where farming activity increases along the gradient. Hence, instead of making disappear the last remnants of TDF in the study area, it would be desirable to increase the productivity of already existing farmland. The decision on which crops (annual herbal and grass species) to grow and how much fertilizer to use could be further improved by reliable climate predictions on the occurrence of ENSO episodes. Other studies have also confirmed the importance of soil for the composition of the vegetation in drylands. For example, Ulrich et al. (2014) found that soil fertility and variability in rainfall and temperature are the most important predictors for beta diversity. Similarly, soil plays an important role in Colombian dry forests (Gonzalez-M et al. 2018) Finally, Espinosa et al. (2011) found out that soil moisture, soil temperature and nitrogen explained in large parts the floristic composition of Tumbesian dry forests, which is quite close to our study area.

Conclusions
Alpha and gamma diversity increased with wetter conditions during the humid La Niña (2012) and moderate El Niño (2016) episodes, though, a few dominant species largely prevailed under the very humid conditions of the coastal coastal El Niño episode (2017), which even resulted in a slightly decreased alpha diversity. Still, plant cover was by far greatest during the most humid conditions. Beta diversity also increased with wetter conditions but did not change significantly between years in accordance with the test for multivariate homogeneity of groups dispersions. As useful as such global measures are in general, they are in essence a single descriptive figure which naturally lack the capability of capturing spatial diversity patterns. By contrast, our spatial predictions revealed drastic changes in the distribution of species composition, beta diversity and plant cover. In effect, a large water surplus disrupted the strictly ordered vegetation formations which at first became apparent along the climatic gradient with wetter conditions. Hence, spatial predictions might provide more detailed insights into ecosystem dynamics, especially in semi-arid ecosystems where changes in plant cover and distribution patterns are more pronounced than changes in alpha and gamma diversity.
The influence of topographic and edaphic variables on plant diversity and productivity also gained in importance with wetter conditions but this effect was not further amplified with very humid conditions. Our irrigation-fertilization experiment partly confirms these findings but also revealed the true beneficial effect of fertilizing under wet Super Niño conditions when plant cover was three times as high as in the unfertilized plots.
In conclusion, our predictive mapping approach can identify locations suitable for restoration and conservation, and thus could form the basis for an informed conservation management (forest management plans). At the same time, the results of our experiment can equip farmers with information regarding sustainable agrarian management which might increase the productivity on already existing farmland which in turn might help to prevent the destruction of the last remnants of TDF in the study area.