In jannes-m/ENSO: Spatio-temporal change of biodiversity along a pronounced gradient heavily influenced by ENSO.

library("knitr")
library("kableExtra")
library("tidyverse")
library("here")

# needed for dates in experiment figure species richness, otw months are displayed in German
Sys.setlocale("LC_TIME", "C")

# make sure that figures are allowed to float (hold the assigned position and are not always put on top of a page)
knitr::opts_chunk$set(fig.pos = 'H')

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 [@macarthur_geographical_1984; @whittaker_gradient_1967]. Naturally, water availability drives floristic composition and biomass production in (semi-)arid regions such as NW Peru [@muenchow_coupling_2013; @muenchow_woody_2013], 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_distribution_2011;@espinosa_what_2011;@salazar_tree_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_understanding_2015]. ENSO consists of alternating episodes. El Niño represents the warm and La Niña 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_monitoring_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. \@ref(fig:study-area-map)D; @rollenbeck_climatic_2015].

Especially El Niño episodes provoke drastic ecological responses [@Kogan.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 eco-climatic 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 21^st^ century [@kosaka_recent_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_enso_2018] though the specific response to El Niño episodes may vary from region to region [@phillips_changes_2009]. Overall, the effect of extreme events is even more pronounced in (semi)-arid ecosystems [@holmgren_extreme_2006], and is especially true for ENSO's terrestrial core region in NW Peru and SW Ecuador [@muenchow_woody_2013; @rollenbeck_climatic_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_what_2011;@pennington_tropical_2018]. At the same time, they are the most threatened tropical ecosystem [@miles_global_2006], and in general chronically understudied [@escribano-avila_biodiversity_2017;@muenchow_review_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_predictive_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_mapping_2020]. Furthermore, predictive mapping can reveal ecological changes in space and time which otherwise might be hard or impossible to detect [@lovelace_geocomputation_2019a]. 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_biogeographic_2017;@esquivel-muelbert_seasonal_2017]. 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 [@laurance_influence_2010; @soethe_nutrient_2008; @condit_species_2013]. Yet, the corresponding effects on the floristic composition of tropical (semi-)arid vegetation communities remain scarcely investigated [@muenchow_soil_2013; @pena-claros_soil_2012;@ulrich_climate_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 [@kinugasa_increasing_2012; @whitford_effects_2011]. 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_drought_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 effects of soil characteristics and topography changed in shaping 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.

Methods

Study area

# knitr::include_graphics(here("figures/12_map_grid.png"))
suppressPackageStartupMessages(library("raster"))
library("png")
back = as.raster(readPNG(here("figures/12_ternary_background.png")))
front = as.raster(readPNG(here("figures/12_ternary_foreground.png")))
p = readRDS(here("figures/12_map_grid.rds"))
p

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. \@ref(fig:study-area-map)). 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_coupling_2013; @muenchow_woody_2013]. 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_monitoring_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. \@ref(fig:study-area-map)d). 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.

precip = readRDS(here("images/12_precip.rds"))
rain_17 = group_by(precip, station) %>% 
  # compute the median for each station
  mutate(median = median(precip)) %>% 
  # just keep coastal EN year
  filter(year == 2017)

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. \@ref(fig:study-area-map)A). 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 is the only studied year during which the whole study area experienced exceptionally wet conditions: Paita received almost r rain_17[rain_17$station == "Paita", "dev"] %>% pull %>% round times (r rain_17[rain_17$station == "Paita", "precip"] %>% pull %>% round mm), Piura almost r rain_17[rain_17$station == "Piura", "dev"] %>% pull %>% round times (r rain_17[rain_17$station == "Piura", "precip"] %>% pull %>% round mm) and Chulucanas r rain_17[rain_17$station == "Chulucanas", "dev"] %>% pull %>% round times (r rain_17[rain_17$station == "Chulucanas", "precip"] %>% pull %>% round mm) more rain than on median (Fig. \@ref(fig:study-area-map)D). 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 [@holmgren_extreme_2006; @muenchow_coupling_2013; @richter_monitoring_2005].

Data compilation

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. \@ref(fig:study-area-map)). 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 (see 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_estimating_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_decimal_1976]. The nomenclature follows the conventions of the Tropicos online database [@missouribotanicalgarden_tropicos_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. \@ref(fig:study-area-map)A). 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 r rain_17[rain_17$station == "Paita", "median"] %>% pull %>% round mm (Paita; annual median precipitation) to r rain_17[rain_17$station == "Chulucanas", "median"] %>% pull %>% round mm (Chulucanas; Fig. \@ref(fig:study-area-map)A). The environmental conditions are typical for the region in terms of climate, topography and soil properties (see @muenchow_soil_2013; @muenchow_woody_2013). 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/98: 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. \@ref(fig:study-area-map)D). 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 200kg/ha 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 [@missouribotanicalgarden_tropicos_2018] and corresponding coverage (%). Please refer to Appendix 2 for a detailed description of the irrigation-fertilization experiment.

Environmental variables

We randomly took and mixed three soil samples from 15cm and 30cm 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\textsubscript{2} emissions, by a factor of 100, especially in the vicinity of trees [@salazar_tree_2020]. Subsequently, we measured in the laboratory variables known to be important for vegetation development: pH, electrical conductivity (µS), carbon-to-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_woody_2013 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 [@landprocessesdistributedactivearchivecenter_lp_2018]. We downloaded the normalized difference vegetation index (NDVI) raster imagery from the Moderate Resolution Imaging Spectrometer (MODIS; @landprocessesdistributedactivearchivecenter_lp_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. \@ref(fig:study-area-map)d) operated by the University of Piura near Paita (1997-2018), in Piura (1991-2018) and in Chulucanas (1997-2018).

Data analysis

Biodiversity and plant cover analysis

We calculated three measures of biodiversity for each observed year, namely alpha, beta and gamma diversity [@whittaker_gradient_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 to 1.4 SD units already correspond to half a change in species composition [@legendre_numerical_2012]. Moreover, the floristic composition was assessed with the help of a DCA [@hill_detrended_1980] to which the decimal scale transformed cover values of the species-plot matrix were subjected [@londo_decimal_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_multivariate_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_coupling_2013 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:

$$ \begin{aligned} \text{Scores}{i}\ \sim\ \mathrm{N}\left( \mu{i},\ \sigma^{2} \right) \ \mu_{i} = \alpha + f_{\text{year}{i}}\left( \text{NDVI}{i} \right) + \beta_{1}\text{year}_{i} \end{aligned} $$

(Formula 1)

where i refers to the ith observation, $\alpha$ is the intercept, $\beta$~1~ is the estimated coefficient for year and f is a spline smoother for the NDVI:year interaction. Please refer to Appendix 3 for a visualization of the non-linear relationship between the response variable (ordination scores) and NDVI by year.

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 section \@ref(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_introduction_2013; @kuhn_applied_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_spatial_2012;@lovelace_geocomputation_2019a]. 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_hyperparameter_2019]. We ran a 100-repeated 5-fold spatial cross-validation in which the partitioning is based on k-means clustering (k = 5; @russ_data_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_improving_2018]. We used the normalized root mean-square error (NRMSE) in percent as performance measure for the spatial cross-validation.

$$ NRMSE = \frac{\sqrt{\frac{\sum_{i = 1}^{n}{(e_{i} - y_{i})²}}{n}}}{max(y) - \ min(y)}*100 $$ (Formula 2)

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 [@borcard_numerical_2018; @peres-neto_variation_2006], 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_forward_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.

All statistical analyses were conducted in the open source software R [@rcoreteam_language_2019] using its packages mgcv [@wood_generalized_2017], RQGIS ([@muenchow_rqgis_2017a], sperrorest [@brenning_spatial_2012] and vegan [@oksanen_vegan_2019].

Results

Biodiversity and plant cover

div = readRDS(here("images/13_desc_stats_tab.rds"))
cap = paste("Different biodiversity measures across the studied years.",
            "min = minimum observed number of species over all plots.",
            "max = maximum observed number of species over all plots.")
kable(div, "latex", booktabs = TRUE, escape = TRUE, 
      # align = "c",
      # col.names = linebreak(nms, align = "c"),
      caption = cap)

Mean alpha diversity was lowest during the dry La Niña (2011) and highest during the moderate El Niño episode (2016; Table \@ref(tab:desc-stats)). 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. \@ref(fig:spri)). Species turnover (beta diversity), expressed as the range of the first DCA axis (see also Fig. \@ref(fig:dca)), 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 \@ref(tab:desc-stats)). By contrast, the latter showed by far the highest plant coverage along the entire gradient (right panel of Fig. \@ref(fig:spri) & section \@ref(spatio-temporal-mapping)).

# include_graphics(here("figures/13_spri_cover.png"))
library("lattice")
library("latticeExtra")
my.panel.loess = function(x, y, span = 3/5, degree = 0.5, ...) {
  loess.fit = loess.smooth(x, y, span = span, dgree = degree )
  panel.lines(loess.fit$x, loess.fit$y, ...)
}
pal = rev(c("lightblue",  "#FDD39B", "pink", "lightgray"))
a = readRDS(here("figures/13_spri.rds"))
b = readRDS(here("figures/13_cover.rds"))
plot(a + b)

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. \@ref(fig:dca)). 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 \@ref(tab:desc-stats)).

# include_graphics(here("figures/14_londo_dca.png"))
library("lattice")
load(file = here("figures/14_londo_dca.rda"))
print(p_1)

Spatio-temporal mapping

Our model explained 76% of the observed variance with a spline smoother absorbing 3 df (Formula 1 & section \@ref(modeling-and-predictive-mapping)). The spatially cross-validated NRMSE also indicated a satisfying fit of 27% (SD 8.53% over 100 repetitions).

# include_graphics(here("figures/15_pred_map.png"))
library("sp")
library("sf")
load(here("figures/15_pred_map.rda"))
print(p)

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. \@ref(fig:pred-map)) 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. \@ref(fig:pred-map)). 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. \@ref(fig:pred-map)). This changed during the wet La Niña episode (2012) where three distinct vegetation formations developed (Panel 2012 of Fig. \@ref(fig:pred-map)). 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. \@ref(fig:pred-map)) . This formation developed into an open xeric shrubland (green/yellow pixels in Fig. \@ref(fig:pred-map)) where Acacia macracantha and Encelia canescens were the dominating species. Further east, the shrubland turned into a TDF formation (orange/red pixels in Fig. \@ref(fig:pred-map)) 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. \@ref(fig:pred-map)). 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. \@ref(fig:pred-map)). 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 increased drastically, and particularly in the western part of the study area where highest coverage values were observed for the grasses Aristida adscensionis and Eragrostis cilianensis and the seedlings of the shrub Encelia canescens. Among herbal species, Crotalaria incana, Exodeconus maritimus and Tiquilia paronychioides showed the highest cover values.

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. \@ref(fig:varpart)). The explained variance roughly doubled in all other years, however, did not further increase with an excessive supply of water in the most humid year (2017; Fig. \@ref(fig:varpart)).

# include_graphics(here("figures/16_varpart.png"))
# Variation partitioning for the four londo cts
# attach packages
library("vegan")
library("dplyr")
# attach own functions
devtools::load_all()

load(here("images/16_input_varpart_plot.rda"))
par(mfrow = c(2, 2))
par(mar = c(0, 0, 0, 0))
#par(tcl = -0.25)
#par(mgp = c(2, 0.6, 0))
for (i in seq_along(inp)) {
  mod = ordi_varpart(inp[[i]], choices = 1:2, expl_1 = soil, expl_2 = topo,
                     dca = TRUE)
  # showvarparts(2, bg = c("hotpink","skyblue"))
  # plot(mod$vp, bg = c("hotpink","skyblue"))
  labs = mod$vp$part$indfract[, 3]
  labs = round(labs, 2)
  labs[c(1, 3)] = paste(labs[c(1, 3)],
                        c(as.character(mod$a), as.character(mod$c)))
  plot(mod$vp, labels = labs, bg = c("hotpink","skyblue"),
       Xnames = "", cex = 1)  
  text(x = 0.5, y = 0.8, labels = years[i], font = 2, cex = 1)
  text(x = c(-0.4, 1.1), y = c(0.4, 0.4), labels = c("soil", "topography"),
       cex = 1)
}

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. \@ref(fig:irrigation)). 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. \@ref(fig:irrigation)). 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. \@ref(fig:spri-exp)).

# include_graphics(here("figures/17_nutrient_irrigation_experiment.png"))
load(here("figures/17_nutrient_irrigation_experiment.rda"))

# to plot manually axes labels outside of the panel, we have to disable the
# clipping beforehand
my.setting = trellis.par.get("clip")
my.setting$panel = "off"
trellis.par.set("clip", my.setting)
# trellis.par.get("clip")

# plot
print(p_1)

spri_exp = readRDS(here("images/17_irrigation_spri.rds"))
print(spri_exp)

Discussion

Biodiversity and biomass

Beta diversity (expressed as the range of the first DCA axis), increased with wetter conditions (Fig. \@ref(fig:dca) & Table \@ref(tab:desc-stats)). However, another measure of beta diversity, the multivariate dispersion, did not change significantly between years (section \@ref(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 & 2016 of Fig. \@ref(fig:dca)) and the very humid neutral year (panel 2017 of Fig. \@ref(fig:dca)) can be most likely attributed to the increase of plant cover (see also next paragraph). Overall, higher species turnover was mainly driven by short-lived, less abundant species. @cleland_sensitivity_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 section \@ref(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 (see Figs. \@ref(fig:irrigation) & \@ref(fig:spri-exp)). 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_why_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_mapping_2012]. This is in contrast to other dryland studies which showed that dormant and annual species increased biodiversity under more humid conditions [@gutierrez_variation_2000; @holmgren_extreme_2006]. Finally, free-ranging livestock might also have an impact on vegetation dynamics in the study area, especially in its eastern parts with more than 100mm 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_nino_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_demographic_1997; @brown_reorganization_1997; @sitters_rainfalltuned_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_climatic_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_water_1990; @salazar_effect_2018; @throop_response_2012].

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_timelagged_2018]. For instance, @dudney_lagging_2017 observed an increase in grass and a decrease in forb abundance when rainfall was high one year earlier. Similarly, @walter_plants_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. \@ref(fig:pred-map)) 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_detrended_1980; @vonwehrden_pluralism_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_effect_2018]. 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 (section \@ref(spatio-temporal-mapping) & Fig. \@ref(fig:pred-map)). 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_what_2011; @salazar_effect_2018]. 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 [@pena_assessing_2015; @salazar_intraspecific_2018].

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_geocomputation_2019a (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 sections \@ref(biodiversity-and-plant-cover) & \@ref(biodiversity-and-biomass)).

Finally, the spatial prediction of species richness (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 (section \@ref(influence-of-environmental-variables) & Fig. \@ref(fig:varpart)). 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_coupling_2013; @muenchow_woody_2013]. Similar findings were observed in Mexico where topography also partly explained the composition of dry forests [@mendez-toribio_effects_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_effects_2011; @whitford_effects_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. \@ref(fig:irrigation)). 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_climate_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_disentangling_2018] Finally, @espinosa_what_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.

References {#references .ListParagraph}

jannes-m/ENSO documentation built on Sept. 21, 2020, 4:37 p.m.