## 2 Materials and Methods
### Data Analysis  
To analyze plant production (peak NDVI per unit area), we used hierarchical generalized additive models (HGAMs) fitted separately to each vegetation type, using the \"GAM\" function available in the R mgcv package [49]. The models included production (peak NDVI) as response (dependent) variable; phenology (time of phenological peak, as number of days counted from 1 September), rainfall accumulated during the growth season (i.e., until the phenological peak), the interaction between phenology and rainfall, and the density (number of adults per unit area) of each of the four ungulate species, as independent, continuous covariates; and management areas as random factors. To account for temporal autocorrelation, we also included \"year\" as a continuous covariate. All models used a Gaussian error distribution and an identity link function. These models were subjected to automatic selection of covariates using the model argument \"Select\" [50]. The predictor variable \"phenology\" (i.e., the day counted from the beginning of the hydrological year at which the phenological peak was observed in the NDVI curve) was taken to represent a surrogate of the evolution (i.e., the specific timing) of rainfall and temperature on that specific hydrological year. This variable was complemented with a second predictor variable, \"accumulated rainfall\" (i.e., the amount of rainfall accumulated until the moment at which the phenological peak took place), taken to represent a surrogate of the water resources available during that year's specific growth period. Hence, the combination of these two variables represents the time and water resources available during the period that preceded the peak NDVI value used as the surrogate of plant production (i.e., as a response variable). The different timing of the (wild and domestic) ungulate censuses vs. plant production estimates introduced a temporal uncertainty concerning which ungulate density values best reflect their impact on plant production. This uncertainty arises because plant production values reflect a biomass-accumulation process taking place from September to peak biomass, in late winter-spring (February-June); whereas ungulate density values during this period fall between measured values of the previous year (before the onset of the growth period) and those at the end of the growth period (June-September). Because all births and most deaths primarily take place during this period (during the summer, most animals are confined for animal-health control and artificially fed), we expected end-of-the-season density values to be more representative of ungulate impacts on vegetation than previous-season values. To objectively test for this assumption, however, we fitted separate models with either of the values, and selected the best-performing models based on the AIC score. Residuals' autocorrelation, if not adequately addressed, may result in biases in parameter estimates and tests of significance. Spatial autocorrelation was addressed by introducing the random factor \"management unit\", which accounted for the spatial dependence of observations in both the dependent and independent variables. Temporal resolution was addressed by (1) including the covariate \"year\" in the model, and (2) computing estimates of the autocorrelation function of each model, using the Acf functions of R's forecast package ([PERSON] and [PERSON], 2007; [PERSON] et al., 2020) [51; 52]. These plots indicated that autocorrelation was adequately dealt with by the models, i.e., there was very low (and, in virtually all cases, non-significant) levels of residuals' autocorrelation. Finally, we used the function vis.gam in the mgcv package to produce contour plots of model predictions on primary production, for different combinations of ungulate abundance (only those with significant effects) and climatology (separate plots for rainfall and phenology). Based on these plots, we performed a visual comparison of historical data and climate change scenarios.