We set up experimental plots in August 2007 to exclude medium and large bodied mammals in two contiguous areas within the same site in the Brazilian Atlantic Forest: Itamambuca (45°5′16″W; 23°19′29″S; hereafter non-defaunated forest), and Vargem Grande (45°14′39″W; 23°26′16″S; hereafter defaunated forest, Fig. 1). Both sites are located in the northern portion of the Serra do Mar State Park, a continuous area of 3154 km2 of Atlantic forests (Ribeiro et al., 2009). The study sites are located only 15 km apart but are separated by a heavy-traffic highway that minimizes movements of mammals between the sites. The sites thus share geology, elevation and precipitation regimes, providing an excellent natural experiment to test local and landscape-wide effects of defaunation. Both forests have similar mammalian communities composition (including agoutis (Dasyprocta spp.), pacas (Cuniculus paca), brocket deer (Mazama spp.) and collared peccaries (Pecari tajacu), amongst other species), but the non-defaunated forest has larger biomass of terrestrial mammals −147.37 kg/km2− than the defaunated forest −46 kg/km2− (Galetti et al., 2017; Rocha-Mendes et al., 2015). In addition the white-lipped peccary (hereafter WLP) represents 93% of the crude mammalian biomass in the non-defaunated forest but goes undetected at the defaunated forest (Rocha-Mendes et al., 2015).

Fig. 1.

(a) Location of non-defaunated (1) and defaunated (2) sites in the continuous Atlantic rainforest. Both sites are located at Parque Estadual Serra do Mar, São Paulo, Brazil. The areas are 15 km apart and split by a highway (BR383). (b) Experimental design. 15 paired 5 m × 3 m open control and closed plots were set up at every site. Within each plot, damage of artificial seedlings, seedling recruitment, plant biomass and species richness were measured in three different subplots of 1 m2 every 6 months. Forest structure was measured in a 10 m radius from the center of the experimental plots.

(0.29MB).

At each site, we set up 15 pairs of exclosure and control plots (intra-pair distance <5 m, between pair distance >200 m), under closed forest canopy away from any type of forest edge (Fig. 1). Pairs of plots were spaced at least 200 m and located under closed forest canopy away from any type of forest edge including forest gaps, edges or dirt roads. Rectangular exclosure plots of 3m × 5m were constructed with 1 m high wire mesh attached to wooden stakes at the periphery of the plot. The mesh had an opening of 5 cm and prevented large and medium-sized mammals (agoutis, pacas, deer, peccaries, tapirs) but not small mammals (<200 g, rodents and marsupials) or arboreal species from entering the plots. Control plots of the same dimensions were marked with 0.5 m high wooden stakes at the corners without mesh as to allow free access. Each plot was divided into eight 1m × 1m subplots with a 0.5 m rim for us to work all around the exclosure without trampling emerging seedlings (Fig. 1). In each plot, three subplots were randomly assigned to measure trampling of artificial seedlings, while the rest of measurements were followed in other three randomly assigned subplots. All plant cover from exclosure and control plots was removed at the beginning of the experiment to accurately estimate ANPP and to overcome differences in regeneration histories between sites.

We characterized structural variables known to drive seedling recruitment (Rother et al., 2015): (1) canopy cover (using a densiometer), (2) depth of litter layer (e.g. number of leaves perforated by a fine metal bar with a sharp end buried vertically into the ground with a single firm movement), (3) abundance of bamboo (Guadua, Merostachys or Chusquea spp.), (4) the density of trees with diameter at breast height (DBH) ≥ 30 cm, (5) density of adult palms of Euterpe edulis (the dominant palm in the area), and (6) the diversity of tree species. Canopy cover and depth of the litter layer were estimated at the center of each plot, all other variables were recorded within a circular plot with 10 m of radius centered at each plot. The gradient in palm density discussed throughout the manuscript was not established on purpose to test our hypothesis, but rather emerged from random sampling of the forest at locations where the experimental plots were placed, thus representing the underlying natural gradient and frequency distribution of palms density found on our forest sites.

In order to estimate herbivores activity in our plots, we systematically distributed 100 artificial seedlings (Clark and Clark, 1991) within the three subplots assigned for this experiment. An artificial seedling consisted of a 40 cm straight piece of galvanized wire (gauge 14) with a plastic flag attached to one end. Each artificial seedling was buried 10 cm into the ground leaving 30 cm exposed. Every six months, from February 2008 and until August 2009, we counted the number of artificial seedlings trampled by mammals at one different randomly assigned subplot with artificial seedlings. We considered an artificial seedling as damaged when the flag was less than 15 cm above the ground either because the wired was bend, buried or the combination thereof (Roldán and Simonetti, 2001).

WLP and tapir presence at experimental plots was also assessed using camera-traps deployed on study sites between February and July 2008. Fifteen cameras were located close to experimental plots (∼ 10 m) and operated with an equivalent sampling effort across sites (416 camera-days per site).

Seedling recruitment was also assessed three times, one subplot in each plot every six months from February 2008 until August 2009. From the randomly assigned subplot we harvested all newly grown seedlings and saplings. Every harvested plant was morphotyped (by a voucher specimen collection) and dried. We counted the number of new individuals and species richness, calculated diversity, and seedling productivity measured as the overall plant dry mass (g/m2). Experiments in controlled environments show that tthough seedlings from tropical trees can develop from seed reserves on early stages of development, these reserves can be exhausted as early as 40 days after planting up to 90 days for 90% of species (Green and Juniper, 2004; Ichie et al., 2001), which is half of the period between successive samplings in our experiment. Indices of alpha diversity such as Shannon’s H and Simpson’s D are non-linearly related to the number of species or abundance of individuals, and hence preclude deriving conclusions about the magnitude of the impact of herbivores based on statistical tests (Chao et al., 2014). Instead we calculated Dstar, Hstar and Estar, as indicators of effective number of species (also called Hill numbers or true diversity), Shannon diversity, and evenness, respectively (Chao et al., 2014; Mendes et al., 2008). Such indices provide unbiased equivalents with standardized magnitudes amenable to statistical comparisons using linear models (Mendes et al., 2008), and represent a gradient where the relative abundance of every species in the community is increasingly factored in, so as to emphasize the dominance and evenness patterns of the community.

Forest structure was summarized with a principal component analysis based on the variance and covariance matrix and scores of the first and second principal components (hereafter PCAF1 and PCAF2, respectively) were used as surrogates of these drivers in the analytical models.

To investigate the response of large herbivores to the density of adult palms and to other known drivers of seedling recruitment we modeled the difference in trampling scores between control and exclosure pair plots as a function of micro-environment descriptors using generalized mixed-effects models (GLMMs (Pinheiro and Bates, 2002)). We fitted two sets of alternative explanatory variables in the fixed component of the GLMMs. A first set included the interaction between the first descriptor of forest structure (PCAF1) and site (defaunated vs. non-defaunated), the second descriptor of forest structure (PCAF2) and site, and time as an additive covariate (factorial: 6, 12 or 18 months). A second set included palms instead of PCAF2. The logic behind these models is that herbivore activity could respond to PCAF1, PCAF2 and palm density simultaneously. However preliminary analyses showed PCAF2 strongly correlated with palm density (R2 = 0.75) but without completely accounting for effects of other structural variables on PCAF2 (such as canopy cover, see results section). Hence by specifying alternative models with either PCAF2 or palm density we avoided co-linearity issues without losing insight into the mechanisms behind the patterns observed. Simplified versions of these alternative models to the data were combined into the same model set so as to assess model fit based on AICc scores, and model-averaged estimates are reported accordingly (Burnham and Anderson, 2002). Models included random intercepts for every experimental plot pair location, and normally distributed errors.

In addition, we tested if WLPs and tapirs responded to forest structure and palm density using camera trap records. We modeled WLP or tapir presence as a response and PCAF1 plus PCAF2 or PCAF1 plus palm density (as explained above) as explanatory variables. In this case, we used Generalized Linear Models (GLMs) and time was excluded as a covariate to avoid overpamaterizing the models and convergence problems that were identified during data exploration. For these models, we specified a binomial error structure, as analyses showed no evidence of overdispersion. Again, AICc was scores were used so as to assess the fit of simplified versions of these models and to derive model averaged estimates.

For every response variable (seedling recruitment, species richness, Hstar, Dstar, Estar and seedling productivity) and experimental pair plot location, we calculated weighted estimates of the variable being analyzed as the difference between the control and exclosure plots divided by their average. This allowed us to test simultaneously the relative magnitude (%) and direction of herbivore impacts whilst accounting for inherent differences between and within defaunated and non-defaunated sites. Positive departures from zero indicate higher values of the metric in control plots and negative departures indicate higher values in the exclosure plots. We followed an analogous statistical modeling strategy as for trampling. Weighted differences in every response variable were modeled as a function of two sets of alternative explanatory variables in the fixed component of the GLMMs. A first set included the interaction between PCAF1, trampling and site, plus the interaction between PCAF2, trampling and site, and time as an additive factorial covariate. A second set included palms instead of PCAF2. Both sets were combined to assess model fit and derived model averaged estimates based on AICc scores.

To test the influence of WLP and tapir on every response variable we fitted again several sets of GLMMs. Because camera trap sessions took place between 6 and 12 months from the beginning of the experiment, we excluded all records for all variables from time = 18 months from this last analysis. Because tapir camera records were scarce, we compared the effects of WLP presence against the additive effects of WLP plus tapir presence (a dummy variable hereafter referred as “TWLP”), so that the net additive effects of tapirs could be inferred from comparing the fit of models including one or the other variable. This prevented spurious correlations arising from small sample size detected during data exploration. A first model set included the interactions between WLP presence/absence and PCAF1, WLP and PCAF2, and site as an additive covariate. A second set included palms instead of PCAF2. The third and fourth sets included TWLP instead of WLP. Time was excluded as a covariate to avoid overpamaterizing the models and convergence problems that were identified during data exploration. All sets were combined to assess model fit and derived model averaged estimates based on AICc scores.

We performed additional analyses to test if impacts of large herbivores on seedling productivity changed with the level of productivity, e.g., if herbivores might lead to changes in productivity in relation to baseline “background” levels in the absence of herbivores. Using GLMMs, we modeled the absolute difference between open and exclosure plots (on a logarithmic scale) as a function of the interaction of trampling, site, and seedling productivity at the exclosure plot (hereafter referred as “background productivity”, also on a logarithmic scale). Additionally, an equivalent set of models were fitted including WLP or TWLP presence/absence as predictors instead of trampling to test how their interaction with background productivity influenced seedling productivity. In these models site was included as an additive covariate to avoid overparameterisation and convergence problems. We used AICc to assess model fit and derive model averaged estimates based on scores.

To test for landscape-scale modulation of seedling productivity, we tested whether mean differences and variances in absolute seedling productivity between open and exclosure plots were different for the defaunated and non-defaunated site. We also estimated the degree of kurtosis in the distribution of these differences. Smaller variance and higher degree of kurtosis in the non-defaunated site would support the hypothesis of landscape-scale modulation by large herbivores.

All GLMMs described above included random intercepts for every experimental plot pair location, and normally/binomial distributed errors as required. No collinearity issues were detected between explanatory variables included in the same model. Previous to performing model selection, all umbrella models were validated using exploratory model fit tools (Zuur et al., 2009). All statistical analyses were done in the Program R using the packages nlme, lme4, MuMIn and vegan (Barton, 2015; Bates et al., 2015; Oksanen et al., 2017; Pinheiro et al., 2019; R Core Team, 2020).

Principal component analysis of environmental variables around experimental plots accounted for 83% of the variability in the first two principal components. Structural variables with the highest correlation coefficients with the scores in the PCAF1 were canopy cover (0.71) and density of trees (0.69), both increasing with increasing values of PCAF1; for the PCAF2 these were the density of Euterpe edulis palms (0.64), canopy cover (−0.60) and density of trees (-0.48), these two decreasing with increasing values of PCAF2. All other correlation coefficients were an order of magnitude lower than those above (see Table S1, Appendix S1), and there was no substantial variance in tree diversity between plot locations. Density of palms and PCAF2 were highly collinear (R2 = 0.75).

Large herbivores showed a strong response to the density of E. edulis palms. Best models included an interaction between site and palm density (Table 1), and suggested that the activity of large herbivores strongly increased with increasing palm density on the non-defaunated site (model averaged estimate [standard error] = 0.804 [0.340], Fig. 2). The site*palm density interaction estimate indicated that this response almost disappeared on the defaunated site (−0.735 [0.309], Fig. 2). The influence of PCAF1 on trampling was not consistent across the set of models considered (estimate [standard error] = 1.32670 [1.420]) and models accounting for effects of PCAF2 performed much worse. Tests also supported an additive effect of time. Like the full community of herbivores, both WLP and tapir presence showed strong positive responses to palms and PCAF2, and WLPs also to PCAF1 (Table 1 and Fig. 2).

Fig. 2.

Response of large herbivores to the gradient in adult palm density: (a) all herbivores, non-defaunated and defaunated sites; (b) white-lipped peccary presence; (c) tapir presence. For statistics see Table 1 and results section.

(0.16MB).

Differences in seedling recruitment between controls and exclosures changed in response to defaunation and the forest structure gradient described by PCAF2, with some minor influence of PCAF1 and trampling, and no effect of time whatsoever (Table 2). The best model, including only site effects, suggested a positive effect of landscape-level defaunation on seedling recruitment: on average, on the non-defaunated site exclosures had 61% more seedlings than controls, but on the defaunated site this difference dropped to 6%. However, a closely competing model included an interaction between PCAF2 and defaunation (Table 2). Model averaged estimates accounting for this uncertainty suggested that on non-defaunated locations with large negative PCAF2 micro-environmental conditions (few palms and closed canopy), seedlings at exclosure plots almost doubled those at controls (93% more abundant) but as palm density increased and canopies became more open this difference decreased, being only 28% more abundant on exclosures (non-defaunated site, intercept = −0.519 [0.233]; slope of PCAF2 effect = 0.342 [0.354], Fig. 3). For defaunated locations, estimates of PCAF2 effects suggested the opposite trend: areas with few palms and large cover showing larger recruitment in controls, and areas with abundant palms and low cover showing larger recruitment in exclosures (difference between defaunated and non-defaunated sites, intercept = 0.566 [0.276]; slope PCAF2 = −0.628 [0.315]). In addition, WLPs affected recruitment though model selection again showed some level of uncertainty (Table 2). Models with palms as predictor fitted better the data than models with PCAF2, and models with additive tapir effects failed to improve substantially the fit of models with WLP alone. Estimates suggested that WLP had large negative effects on recruitment in areas with few palms (100% decline in controls respect to exclosure plots) whilst this negative effect declined in areas with abundant palms (slope palm density effect = −0.020 [0.022]; WLP presence effect = −0.062 [0.670]; slope interaction = 0.1017 [0.056], Fig. 3). In contrast, in areas where WLP were absent the effect of palms on recruitment was negligible.

Fig. 3.

Changes in the response of plant communities to environmental gradients of PCAF2 and adult palm density as a result of experimental defaunation. Responses are weighted differences between control and exclosure plots, where positive values denote larger magnitude in controls and negative values larger magnitude on exclosures (see methods section for more details). The figure shows: % changes in recruitment (a), species richness (c) and Hstar (e) along the PCAF2 gradient for non-defaunated and defaunated sites; % changes in recruitment (b), species richness (d) and Hstar (f) along the palm density gradient for sites where white-lipped peccaries and tapirs where present.

(0.4MB).

Differences in species richness between controls and exclosures changed in response to defaunation, forest structure and trampling, and some specific interactions between these factors, but not time. On average, on the non-defaunated site exclosures had 42% more species than controls, but on the defaunated site controls only had 7% more species than exlosures. The set of best models suggested and interaction between PCAF2 and site, and between PCAF1 and trampling (Table 2). With similar trends as those described for recruitment but less uncertainty, the first interaction indicated that on the non-defaunated site richness was on average 87% larger on exclosure than in control plots in areas with few palms and closed canopy, but as palm density increased and canopies became more open differences between controls and exclosures decreased down to 16% (non-defaunated site, intercept = −0.417 [0.163]; slope of PCAF2 effect = 0.375 [0.250], Fig. 3). In contrast on the defaunated site richness was on average 24% larger on control than in exclosure plots in areas with few palms and closed canopy, but as palm density increased and canopies became more open exclosures held 20% more species than controls (difference between defaunated and non-defaunated sites, intercept = 0.482 [0.190]; slope PCAF2 = −0.506 [0.219]). The second strong interaction (PCAF1 and trampling, Table 2) indicated that in areas with high density of trees and large canopy cover trampling had a substantial positive effect on richness, but in low density and low cover areas, trampling had a large negative effect (PCAF1 effect estimate = −0.021 [0.182]; trampling effect estimate = −0.001 [0.011]; interaction estimate = 0.025 [0.011]. In addition, both tapirs and WLPs affected species richness (Table 2). Results supported a strong TWLP*palm interaction, with e.g. estimates for the palm interaction suggesting an average decline of more than 115% in controls respect to exclosure plots on areas with few palms but an increase up to 25% in controls respect to exclosure plots on areas with abundant palms where tapirs and WLPs were present (e.g. slope palm density effect = −0.009 [0.016]; TWLP presence effect = −0.898 [0.523]; slope TWLP*palm interaction = 0.1113 [0.037], Table 2 and Fig. 3). As for recruitment, the differential effects of tapirs and WLPs was linked to palm density rather than PCAF2, and in areas where these were absent the effect of palms on species richness was negligible.

Defaunation at the landscape scale showed the largest influence on Hstar (equivalent to Shannon’s diversity) and Dstar (equivalent to Effective Number of Species), whilst the influence of forest structure, palm density, trampling or time were not consistent across models (Table 2). Model averaged estimates suggested that on the non-defaunated site, Hstar and Dstar were 53% and 31% larger on exclosures than controls, respectively, whilst for defaunated site differences between exclosures and controls were negligible (7% larger on controls for both metrics).

When considering the effect of WLPs and tapirs, there was a strong effect of the interaction between TWLP (but not WLP) and palms on Hstar, suggesting that tapir more than WLPs had large negative effects on Hstar on areas with few palms (114 % more on exclosures) but positive effects on areas with many palms (Fig. 3). The impact of tapirs and peccaries also changed with PCAF1, with negative effects on Hstar decreasing as canopy cover and density of trees decreased. In the absence of tapirs and peccaries Hstar showed similar values in controls and exclosures regardless of palms, canopy cover or density of trees. For Dstar, there was evidence of a TWLP*palms interaction fitting better than WLP*palms interaction, though the strength of the interaction was lower than for Hstar (Table 2). Model averaged estimates suggested that negatives effects of both tapir and WLP at locations with low density of palms (55% reduction in Dstar) changed to positive at high palm density (20% increase). In contrast to other variables, E-star didn’t show a consistent response to any of the predictors considered, including site.

Seedling productivity responded strongly to defaunation, forest structure and trampling. Model averaged estimates suggested that on the non-defaunated site seedling productivity was on average 36% larger on exclosures than controls, whilst for defaunated site differences between exclosures and controls were negligible (4% larger on exclosures). However, best models supported a strong interaction between site and trampling, and between trampling and PCAF2 (Table 2), and marginal or no support to other predictors. The interactions suggested that at non-defaunated areas trampling had a marginal positive effect on seedling productivity (slope of trampling effects on non-defaunated site = 0.005 [0.020]), but on defaunated areas the effect of trampling was strongly negative (difference in slopes of trampling effects between defaunated and non-defaunated sites = −0.129 [0.042]). The second interaction suggested that, regardless of site, in areas with high canopy cover and low palm density (negative PCAF2), trampling had a negative effect on seedling productivity, whilst on areas with high palm and low cover, trampling had an unexpected strong positive effect (PCAF2 effect estimate = −0.607 [0.465]; interaction estimate = 0.061 [0.027]). The combined effect of these two interactions on seedling productivity is shown in Fig. 4: at non-defaunated plots with large negative PCAF2 values, trampling had a strong negative impact, so that seedling productivity in exclosure plots was on average almost 50% larger than on open controls if trampling was high. Surprisingly, as PCAF2 values became more positive (high palm density and open canopy) the effects of trampling reversed from negative to positive, so that in controls open to herbivores seedling productivity was almost twice that of exclosures if trampling was high. At the defaunated site however, the strong negative impact of trampling dominated, so that increasing trampling had always a negative impact on seedling recruitment regardless of forest structure (Fig. 4).

Fig. 4.

Changes in seedling productivity along the PCAF2 gradient as a result of: (a) experimental defaunation in the non-defaunated site; (b) experimental defaunation on the defaunated site; (c) presence or absence of white-lipped peccaries (left and right, respectively). Responses are weighted differences between control and exclosure plots, where positive values denote larger magnitude in controls and negative values larger magnitude on exclosures (see methods section for more details). For figures (a) and (b) subpanels denote a gradient in PCAF2 (e.g., leftmost subpanel showing largest negative PCAF2 value range, rightmost subpanel showing largest positive PCAF2 value range); within subpanels, the x-axis indicates the effect of the intensity of trampling for every particular range of PCAF2 values.

(0.43MB).

In addition, WLPs but no tapirs had strong differential spatially structured effects on seedling biomass (Table 2). Results supported a WLP*PCAF2 and WLP*palm interaction, with e.g. estimates for the palm interaction suggesting an average decline of more than 80% in controls respect to exclosure plots on areas with WLPs and few palms but an increase up to 25% in controls respect to exclosure plots on areas with WLPs and abundant palms (Fig. 4). In contrast, on areas with no WLPs, increasing palm density decreased seedling biomass in controls (e.g. slope palm density effect = −0.045 [0.026]; WLP presence effect = −1.141 [0.803]; slope WLP*palm interaction = 0.124 [0.070]).

When considering differences in seedling productivity between experimental treatments in the context of a gradient in background productivity, best models supported the influence of all two-way interactions between defaunation, trampling and productivity (Table 2). First, background productivity had smaller negative effects on seedling productivity in plots of the defaunated site (productivity effect estimate on non-defaunated site = −0.617 [0.288]; productivity*site interaction estimate = 0.362 [0.308]). In addition, trampling had larger negative effects on seedling productivity at higher background productivity (trampling*productivity interaction estimate = −0.022 [0.016]), but also more negative effects on seedling productivity at plots on the defaunated site (trampling effect estimate on non-defaunated site = −0.018 [0.027]; trampling*site interaction estimate = −0.098 [0.070]). The combined effect of these three interactions on seedling productivity is shown in Fig. 5: at non-defaunated sites with low productivity, trampling had a strong positive impact, e.g. so that at when trampling was high seedling productivity in control plots open to herbivores almost doubled that of exclosures; however, in non-defaunated sites with high productivity increasing trampling had always net negative impacts on productivity so that seedling productivity in exclosures almost doubled that of controls. At the defaunated site however, increasing trampling always had a negative impact on seedling productivity regardless of background levels of productivity (Fig. 5).

Fig. 5.

The modulating impact of large herbivores. Figures (a) and (b) show changes in log seedling productivity along the background seedling productivity gradient (sprod) as a result of experimental defaunation in the non-defaunated and defaunated sites, respectively. Figure (c) shows the density distribution of these changes (left figure) and the density distribution of changes in log seedling productivity at plots where white-lipped peccaries (WLP) where present or absent (right figure). Responses are weighted differences between control and exclosure plots, where positive values denote larger magnitude in controls and negative values larger magnitude on exclosures (see methods section for more details). For figures (a) and (b) subpanels denote a gradient in background log seedling productivity (e.g., leftmost subpanel showing lowest productivity value range, rightmost subpanel showing largest productivity value range); within subpanels, the x-axis indicates the effect of the intensity of trampling for every particular range of productivity values.

(0.4MB).

Mean landscape-scale differences in absolute seedling productivity between open and exclosure plots were not different for both sites (non-defaunated mean [s.e.] = −0.337 g m−1 y−1 [4.557], defaunated additive estimate = −2.570 g m−1 y−1 [6.288], t = −0.409, P = 0.686). However differences in variances were significant (larger on defaunated, F = 0.014, P < 0.001) and excess kurtosis (“peakedness” in the distribution) was large in the defaunated forest but extreme in the non-defaunated forest (2.49 and 8.23, respectively, Fig. 5).

WLP but no tapirs had a strong role in modulating landscape-level productivity (Table 2). Where WLPs were absent, controls had larger seedling productivity than exclosures at low background productivity but exclosures had larger productivity than controls at high background productivity (intercept: −0.501 [0.5055]; slope background productivity effect = −0.599 [0.216]). WLPs tended to reduce this gradient (WLP presence effect = −0.603 [0.664]; slope WLP*palm interaction = 0.460 [0.289]). As a result, whilst mean landscape-scale differences in seedling productivity between open and exclosure plots were similar in areas with and without WLPs (WLPs absent mean [s.e.] = −0.585 g m−1 y−1 [3.874], WLP present additive estimate = −3.201 g m−1 y−1 [6.599], t = −0.485, P = 0.631), differences in variances were substantial (larger where WLPs were absent, F = 0.412, P = 0.034). Excess kurtosis (“peakedness” in the distribution) was large in locations without WLPs but extreme in locations where WLPs were present (4.65 and 11.17, respectively, Fig. 5).

We tested for spatially structured differential effects of large mammalian herbivores on tropical forest communities along a gradient in density of palms, whose fruits are a key resource for large herbivores and frugivores on Neotropical forests. Our results match those from previous studies (Akkawi et al., 2020; Bodmer, 1990), demonstrating that the density of palms strongly drives large herbivore activity, both at the whole community level (measured through trampling) but also for WLPs and tapirs in particular (through camera traps). WLPs in addition selected areas that had simultaneously larger density of palms and lower canopy cover. It is possible that per-capita palm fruit production in areas with open canopies might increase in response to larger availability of light, which is a limiting resource for plants productivity in many tropical forests and elsewhere (Borer et al., 2014; Tang and Dubayah, 2017; Wright and Carel, 1994). WLPs live in large herds, have large home ranges, move across large tracts of the landscape and might be able to locate, move between and memorize hotspots of productivity (Beck, 2005; Reyna-Hurtado et al., 2009). Thus, it is possible that WLPs might track palm fruit availability across large spatial scales more efficiently than other medium/large terrestrial herbivores with more restricted home ranges and/or solitary habit on these forests. Tapirs, for example, showed a strong response to palm density but only when palm density was very large (clusters of over 12 individuals in a 10 m. radius), suggesting that tapir preference for E. edulis fruits might be high only when these fruits are highly abundant. Likewise, Bodmer (1990) showed that tapirs consumed a large amount of Mauritia flexuosa palms occurring in an extremely aggregated pattern in the Amazon, indicating that tapirs preference for extremely productive palm patches might be conserved across biomes. WLPs, in contrast showed a strong saturating functional response to palms typical of consumers with high searching efficiency but large handling time able to closely track their main resource regardless of its density (Holling, 1965).

By affecting large herbivore activity, the effects of E. edulis palm indirectly influenced the seedling community, which showed concomitant treatment, site and spatially structured tendencies in responses in recruitment, species richness and productivity. The responses detected are consistent with a process of landscape-scale modulating role of large herbivores dominated by their interaction with palms. Ignoring this critical “within-site” interaction and focusing on landscape-average effects would have led to the wrong inference that the effects of large herbivore consumers on plant communities are strongly negative across the whole forest. For example, on the non-defaunated site, landscape-scale average seedling recruitment, species richness, diversity [Hstar/Dstar] and productivity were larger on exclosures than controls (61%, 42%, 53%, and 36%, respectively). On the defaunated site, on the contrary, differences between open plots and exclosures were negligible. However, by accounting for the within-landscape natural gradient in palm density and forest structure, we gained insight into strong underlying spatial herbivore-plant dynamics in tropical forest systems.

The patterns indicate that on areas where large herbivores were more abundant their impact was less negative or null. This suggests that the intensity of top-down limitation of these herbivores on plant communities at the ground level varies spatially, decreasing where their key resources (e.g. palm fruits) are more abundant. For grassland ecosystems, many studies have shown that large herbivores increase productivity and diversity with increasing productivity of resources (Augustine and McNaughton, 1998; Bakker et al., 2006; Frank et al., 1998; Pringle et al., 2007; Hillary S. Young et al., 2013). However for forest ecosystems such patterns haven’t been explored in detail, and the assumption is that large herbivores exert negative effects regardless of resource productivity (Augustine and McNaughton, 1998; Crête, 1999; McLaren and Peterson, 1994; Pastor et al., 1993; Ritchie et al., 1998). Our results reject such tacit assumption. Several non-mutually exclusive mechanisms based on previous empirical studies might help to explain the patterns observed.

First, rather than direct consumption, competitive exclusion between large and small herbivores might have contributed towards these patterns. For example, previous studies in the Atlantic Forest suggest that large herbivores like peccaries reduce the abundance and biomass of small mammals that feed on seed and seedlings (Bovendorp et al., 2019; Galetti et al., 2015). Indeed, small mammal declines in response to ungulate densities is a common pattern in most ungulate-small mammal studies reported worldwide (Daskin and Pringle, 2016; Villar et al., 2014; Young et al., 2015). If the impact of small mammal predation on seeds and seedlings is strong (as suggested by (Paine et al., 2016)), then their exclusion by large tropical herbivores might lead to relative increases in plant recruitment, and as a result, in plant species richness. Since large herbivores preferentially select areas where palms are abundant, exclusion of small mammals and relative increases in seedling recruitment would be expected on such areas, a prediction matched by our results. Indeed, previous work from grassland ecosystems suggest that the negative impact of large herbivores on small mammal densities increases with resource productivity (Daskin and Pringle, 2016; Young et al., 2013), though evidence from forests ecosystems is lacking. On the other hand, some studies from the Atlantic forest strongly suggest that average plant recruitment does not decrease in response to defaunation or herbivore exclusion (Brocardo et al., 2013; Villar et al., 2020), contradicting a critical prediction from this hypothesis. Regardless, further work is needed so as to investigate how large herbivores and small mammals interact across gradients in resource availability in forest ecosystems.

A second prominent mechanism that might explain patterns observed in our study are the effects of large herbivores on nutrient cycling. Recently, Villar et al. ( 2021) demonstrated experimentally that large frugivores like WLPs increase soil nitrogen cycling and nitrogen availability to plants, and that this pattern was also spatially structured across gradients of density of palms. Thus, despite large herbivore activity, trampling and top-down consumption of seeds might increase on areas with abundant palms, by increasing the soil nitrogen pool available to seedlings, large herbivores might stimulate seedling recruitment and productivity on those areas, effectively compensating for their antagonistic top-down effects, compensated by indirect stimulation of bottom-up forcing. On some grassland ecosystems large herbivores promote nutrient availability to plants, improve soil moisture and stimulate plant re-growth, leading to a positive feedback process that can compensate for large herbivore consumption rates (Schrama et al., 2013). Though fundamental differences exist between grassland and woodland systems, it is possible that some of these mechanisms operating in grasslands might also operate beneath the canopy of tropical forests where frugivory is ubiquitous, and hence contribute to increase seedling productivity in hotspots of large herbivore activity under open canopies. WLPs can forage in large herds (from 50 to hundreds of animals), plowing roots, trampling the soil in their search for seeds and defecating so they might be able to influence nutrient cycling and stimulate the growth of some plant species adapted to large herbivore disturbance.

A third mechanism that might explain the positive effects of large herbivores is an indirect facilitation process towards seedlings from species with higher wood density in detriment of palm seedlings. Light wooded species have higher competitive advantage in areas of open canopies due to faster photosynthetic rates (Santiago et al., 2004; Svenning, 2002), so that in the absence of substantial herbivory, palm seedlings may have a competitive advantage respect to harder-wooded seedlings. However, it is possible that increased large herbivore trampling and predation of palm seeds and seedlings under parent trees predicted under Janzen-Connell density dependence (Comita et al., 2014; Connell, 1971; Janzen, 1970; Terborgh, 2012) might lead to proportionally decreased abundance of palm recruits relative to harder-wooded competitors, causing increases in net seedling productivity. Since WLPs exert a strong negative impact on palm seeds and seedlings (Beck, 2006; Keuroghlian and Eaton, 2009; Silman et al., 2008), it is possible that they might indirectly contribute towards increasing recruitment of other species.

Lastly, this mechanism might interact with herbivore-assisted seed dispersal. Our analyses suggest that tapirs had a strong positive influence on diversity on areas with abundant palms, and contributed towards the positive effects of palms on species richness. Thus, by dispersing seeds from other species into palm hotspots tapirs might contribute towards the recruitment of harder-wooded species under palm trees and open canopies. As noted previously, tapirs consume a variety of fruits and therefore are bound to be effective long-distance seed dispersers. Indeed, this hypothesis implicitly reinforces the complementary functional role of WLPs and tapirs, though by highlighting the leading role of WLPs our results might imply that the role of tapirs is subordinated to the presence of WLP activity, as has been previously suggested (Villar et al., 2020).

Despite their effects on recruitment, species richness and productivity, structuring effects on diversity were negligible. Results suggested that tapirs appeared to have some structuring effects on seedling diversity across palm density gradients, but neither WLPs nor the whole guild of large herbivores did, rather reduced landscape-scale seedling diversity average by circa 50%. A decade-long experiment in the Atlantic Forest also suggested that, in the long-term, large herbivores do not have a substantial effect on sapling alpha diversity, but rather increased beta diversity (Villar et al., 2020). Other additional studies in Neotropical forests have shown contradictory results (Kurten and Carson, 2015), and so the debate on whether large herbivores increase or decrease plant diversity in tropical forests remains yet unresolved. It has been argued that such contradictory results might be due to possible differences in the community composition of the medium, small-sized and arboreal herbivore guilds between defaunation and exclosure studies, though it is also possible that the duration of the studies and the ontogenic phase of the plant community inspected might affect the outcome.

We found that WLPs impacted more strongly on seedling recruitment and productivity, contributed towards the positive effect of large herbivores on species richness and had no effect on diversity. WLPs consume large quantities of palm seeds, damage palm seedlings and change soil conditions below parent palm trees (Beck, 2006; Villar et al., 2021). Thus, it is possible that WLPs effects on seedling communities are linked to their strong association to palm trees and might be mediated, at least in part, by their impacts on palm demography. In contrast, we found that that tapir impacts were restricted to diversity, and also contributed towards plant richness. Given their well-documented seed dispersal potential, it is possible that tapirs might have increased seedling species diversity at focal locations by dispersing seeds from their ingested species pool. In contrast tapirs’ solitary habit and low population densities might limit their impact on recruitment and productivity of seedlings. Previous studies in other systems suggest complementary functions of different herbivores (Brodie et al., 2009; Burkepile and Hay, 2008; García and Martínez, 2012), and, on the Atlantic Forest, a long-term positive effect of large herbivores on plant beta diversity has been attributed to the joint impact of WLPs and tapirs (Villar et al., 2020). Because tapirs and WLPs affected different seedling community processes in our experiment, it appears that they might play a complementary function in the regulation of regenerating seedling communities in the Atlantic Forest.

One of the most striking results from our experiment is the discovery of positive productivity feedbacks between large herbivore activity and seedling productivity on areas with abundant palms and low background productivity at the forest where WLPs and tapirs were abundant. Regardless of the mechanism(s) operating (discussed above), these spatially structured feedbacks are far from trivial and comparable to other documented on well-studied grassland systems. For example, minimum and maximum weighted differences in seedling productivity control minus exclosure were −194% and +185%, respectively, whilst herbivores increased landscape-level primary productivity by an average 94% (range = −8%, +344%) in tropical grasslands and 33% (12%, 85%) in temperate grasslands (Frank et al., 2002, 1998; McNaughton, 1985). Furthermore, though absolute estimates of ground-level productivity in control plots at our non-defaunated study site were obviously orders of magnitude lower than for these grassland systems (e.g. 0.62 g/m−1/y−1 [CI = 0.380 1.00], compared to 692.57 g m−1 y-1 and 180.68 g m−1 y−1 in Serengeti and Yellowstone grazed controls), spatially structured impacts of large herbivores on regenerating tropical seedlings will most likely multiply exponentially at later ontogenic stages of development of these seedlings into adult trees, so as to substantially affect the spatial distribution of carbon stocks in mature forests (Bello et al., 2015).

Our results also suggest that spatially structured effects of large herbivores on seedling productivity across gradients of palm density led to the emergence of a process of herbivore modulation of seedling productivity at the landscape-scale. By increasing productivity at the lower end of the productivity spectrum and decreasing productivity at the higher end, large herbivore activity affected the density distribution of seedling productivity, increasing the frequency of intermediate productivity patches and decreasing the landscape-scale variance in productivity at the ground level on the non-defaunated site. Villar et al. (2021) found a very similar pattern of landscape-scale modulation of nitrogen availability by large herbivores in the Atlantic Forest. Like in that study, our results suggest that this process was primarily driven by WLPs, whose herding behavior, large home ranges, frequent large distance movements and large memory might contribute to forest-wide scale regulation of plant communities, nutrient cycling and productivity.

Our study demonstrates that palms are a key driver of the activity of large tropical herbivores, and that large herbivore-palm consumer-resource interactions spatially structure plant communities at the understory of tropical forests. Such findings support the view that palms might act as foundation species in tropical forests, and that their interaction with fruit-consuming large herbivores might be critical in this foundational role. In the absence of such herbivores, the spatial structure of plant communities along palm density gradients collapsed. There is a growing concern for the long-term consequences of illegal hunting, the subsequent loss of seed dispersal function by large vertebrates and their implications for plant recruitment, community composition and carbon storage in tropical systems (Bello et al., 2015; Galetti and Dirzo, 2013; Harrison et al., 2013; Kurten, 2013; Stoner et al., 2007). Our work indicates that increasing pressures placed by simultaneous illegal hunting on large herbivores and illegal harvesting of palms (Galetti and Fernandez, 1998) might have major effects on tropical forests dynamics. It is indeed possible that in the long term these changes might lead to an alternative ecosystem state where plant communities are regulated by alternative processes, as has been evidenced in other ecosystems (James A. Estes et al., 2011). The patterns and mechanisms discussed in this study deserve closer attention, and future studies should inspect whether they are applicable to other forest ecosystems so as to infer generality. Our study also suggests that it is paramount to characterize and integrate how defaunation affects alpha and beta diversity, species turnover and functional traits, nutrient cycling to understand the regulation of spatiotemporal dynamics in tropical forests.

On a wider ecological context our study demonstrates that the sign and magnitude of large mammalian tropical herbivore effects on forest plant communities are strongly spatially structured and change dramatically and predictably across resource gradients. This highlights the importance of considering spatially-structured consumer-resource dynamics to understand trophic cascades, ecosystem processes and overall dynamics in forest ecosystems.