Introduction

As climate change and anthropogenic activity alter ecosystems at unprecedented rates, it has become critical to understand the consequences of biodiversity loss on ecosystem processes and the maintenance of ecosystem processes through species interactions. A complex mix of anthropogenic forces are eroding multiple dimensions of global biological diversity, including plant intraspecific, interspecific, and functional diversity (1, 2). Plant diversity affects herbivore abundance and diversity, thereby influencing biomass allocation and energy fluxes between trophic levels (3). Because losses of plant diversity can destabilize the flow of resources to higher trophic levels (3, 4), understanding the connection between biodiversity and trophic interactions is necessary to predict the consequences of the loss of primary producer biodiversity on ecosystem traits such as resilience to environmental perturbation (5, 6). While climate change is associated with increased absolute precipitation in some regions and decreased precipitation in others, IPCC models predict mid-century increases in the frequency of extreme precipitation events in Central and South America (79), a phenomenon already observable in many ecosystems (10). Measuring how precipitation change will affect the relationship between biodiversity and ecosystem function is therefore increasingly important, particularly in diverse tropical systems.

Very few multi-site, manipulative diversity experiments have been reported from tropical areas compared to temperate environments (11), limiting our knowledge of the role of biodiversity in ecosystem function in the most species-rich regions of the planet (12). Here we investigate how multiple dimensions of plant diversity affect herbivory and plant mortality at five neotropical sites and explore how diversity modulates ecosystem responses to changes in water availability at three of those sites, in the context of several established hypotheses (Table S1). By altering plant uptake of nutrients and the induction of plant defenses, abnormal levels of precipitation can alter herbivore pressure, acting to alter the movement of resources into higher trophic levels. As such, climate extremes of dry and wet periods are expected to strongly perturb plant-insect interactions and thereby alter ecosystem function (1315). The insurance hypothesis suggests that greater biodiversity can act to stabilize ecosystems and improve their resilience to environmental change (4, 16). Despite traditional views that interspecific richness has a greater impact on ecosystem processes than intraspecific diversity (17), recent research suggests that plant intra- and interspecific richness can have similar effects on ecosystem productivity and consumer abundance (18, 19), and field studies have demonstrated both positive and neutral effects of plant diversity on ecosystem resilience (20, 21).

To test the insurance hypothesis, the effects of water addition on ecosystem function, and the relative strength of intra- and interspecific richness in driving ecosystem function, we conducted common garden experiments at five sites across Central and South America. Using 33 species in the genus Piper (Piperaceae) as a model system (22), we manipulated intra- and interspecific plant richness in Costa Rica, Ecuador, Peru, and two sites in Brazil. We additionally manipulated water availability in Costa Rica, Ecuador, and Peru (Figs. 1 & S1 A). We predicted that reduced Piper diversity would lead to reduced diversity in higher trophic levels, that water addition would lead to altered herbivore pressure, and that lower Piper diversity would be associated with more extreme changes in herbivory and plant mortality in response to water addition. Additionally, we predicted that changes in intraspecific plant richness would affect ecosystem function with magnitudes similar to interspecific richness.

Treatments of intraspecific richness, interspecific richness, and water addition used in each of the five study sites.

White tiles represent treatment combinations which were not tested in a given site.

Methods

Study sites and focal plant genus

We conducted a large-scale transplant experiment replicated across five sites spanning 42° latitude in the Neotropics (Fig. 1 & Table S2 and encompassing a range in annual precipitation from 1271 mm to 4495 mm (Table S2). At each site, we studied herbivory on planted individuals in the genus Piper (Piperaceae) in response to experimental treatments. Study sites included lowland equatorial humid forest at La Selva Biological Station, Costa Rica; high elevation equatorial humid forest in Yanayacu Biological Station, Ecuador; high elevation equatorial humid forest at El Fundo Génova, Peru; lowland seasonally dry gallery forest in cerrado within the phytogeographic domain of the Atlantic Forest in Mogi-Guaçu Biological Reserve, Brazil; and lowland seasonally semideciduous forest in the transition between the Atlantic Forest and the Cerrado phytogeographic domains in Uaimii State Forest, Brazil. Climate classifications follow the Köppen-Geiger climate model (Table S2; 22).

Multiple species of Piper are found at all sites, ranging from 11 species in Mogi-Guaçu to 50 species in Costa Rica (24). Piper is an ideal genus for large-scale comparative studies, as it is found across the Neotropics and subtropics, from ∼10° N to about ∼32° S. In addition to being widespread, Piper is abundant and diverse across its range, encompassing ∼1000 species in the Neotropics (25). Piper has been the subject of detailed studies of herbivory, and its herbivore fauna has been surveyed across its range (2637). Members of the genus host both specialist caterpillars and beetles as well as generalist caterpillars and Orthopterans (22,26–29,31). Leaf damage patterns produced by the different classes of herbivores are well documented (31)

Experimental design

At each study site, a factorial experiment was implemented to test the effects of plant intraspecific and interspecific richness on herbivory, variation in herbivory, and insect richness. Experimental plants were propagated from cuttings of naturally occurring Piper plants (typically with three nodes and zero to one leaves). In Costa Rica, Ecuador, and Peru a water addition treatment was also tested. In all sites except Uaimii, experimental plots measured 4 m in diameter and contained 12 Piper individuals, planted either in pots with locally derived soil (for experiments manipulating water availability) or directly into the ground (for experiments without water additions in Mogi-Guaçu and Uaimii; Table S2). Plots in Uaimii were 2 m in diameter and contained six Piper individuals, limited by the number of cuttings available. Plants were cultivated without fertilizer or irrigation.

The interspecific richness treatment had three levels at all sites except Uaimii, where there were only two levels due to limited numbers of cuttings. At each site a subset of naturally occurring Piper species was selected to act as a species pool for each experimental plot (Table S3). These species were selected to constitute a breadth of genetic and functional diversity representative of the Piper community present at each site. Combinations of species in each plot for each level of the interspecific richness treatment were selected randomly from these species pools. The first level was a monoculture of Piper consisting of individuals of a single species. In the second level, two Piper species were randomly selected for planting (e.g. six individuals of each species). These species combinations were the same when crossed with the low and high levels of the intraspecific richness treatment (see below). In the high interspecific richness plots, the number of species varied between sites and included all Piper species available to produce cuttings (Table S1Table S2). The minimum number of Piper species in a high richness plot was three in Uaimii, and the maximum was 12 in Costa Rica and Ecuador. Due to the low number of species available in Uaimii no two-species plots were planted. The interspecific richness treatment was replicated three times per site. In Mogi-Guaçu and Uaimii, each replicate contained a monoculture plot for each species.

An intraspecific richness treatment with two levels (clones taken from a single mother plant or individual cuttings all taken from unique mother plants) was crossed with the interspecific richness treatment. In Costa Rica and Ecuador, the high interspecific richness plots only included individuals from unique mother plants because 12 different species were planted. In Mogi-Guaçu and Uaimii, the experiments were fully crossed. In the high intraspecific richness treatment, cuttings were from different Piper individuals growing at least 10 m apart, so they were not from the same root system are were most likely genetically unique individuals.

A water addition treatment was crossed with the inter- and intraspecific richness treatments in Costa Rica, Ecuador, and Peru in order to examine how IPCC predicted regional mid-century increases in precipitation in these regions (38) may affect herbivory and trophic interactions. Plants were cultivated in five to six liter pots with drainage holes. Each potted plant under the water addition treatment was watered with an extra 2 L of water twice per month in Costa Rica and Ecuador. In Peru, 2 L of water were added to plants every three months. Water was rapidly applied as a flooding event to completely saturate the soil. All low richness plantings in Peru were watered, but in Costa Rica and Ecuador the low richness plantings were not watered, so only two levels of the species richness treatment were crossed with the water treatment (Figs. 1 & S1 A). Plant cultivation periods lasted in excess of one year at all sites, and as such water application treatments were performed in both the wet season and dry season.

Experimental plots were randomly located in the study sites; in Mogi-Guaçu and Uaimii they were organized in three replicate blocks. There was a minimum of 20 m between plots in blocks and 100 m between blocks. Randomly located plots were separated by a minimum of 50 m. The plots were spread across a maximum of 1,500 m in Peru and a minimum of 275 m in Ecuador. Experimental periods lasted between 1.4 and 2.8 years, depending on site (Table S2), and all measurements of herbivory were conducted on leaves that were initiated during these periods. Cuttings were replaced if they died during the first three months of the experiment. One of the authors was present at all of the sites for multiple visits to ensure as much standardization of treatments as possible.

To determine the effects of site variation in natural levels of precipitation on the outcomes of the water treatments, the absolute level of precipitation and precipitation anomaly relative to climate normals was collected for each month of the experimental periods in Costa Rica, Ecuador, and Peru (Fig. S2). Climate data for Costa Rica were obtained from La Selva Biological Station. Data for Ecuador and Peru were obtained using TerraClimate (39), and an interpolation error in the precipitation for Peru in February 2016, was corrected using data from the National Service of Meteorology and Hydrology of Peru (SENAMHI).

Measures of herbivory and insect richness

Herbivory was recorded by taking photographs of all the experimental leaves at the end of the experimental period. Additional photos were taken every three months in Uaimii and Mogi-Guaçu, and in the first five months of the experiment in Costa Rica. These data were used to measure herbivory and to determine the types of herbivores feeding on leaves based on patterns of damage. When herbivores were observed on plants, they were photographed but were left to continue feeding so as to not interfere with the experiment. Herbivory was quantified by eye for each type of herbivore on each leaf by a single parataxonomist with extensive experience measuring herbivory on Piper following established protocols (31). The amount of leaf area consumed was measured in function to the total leaf area by visually dividing the leaf into equal sized segments to determine the percent area missing. This was measured as a continuous value to the greatest possible accuracy. Insect herbivores were identified to the lowest taxonomic level possible based on their damage patterns (genus for specialist Lepidoptera, family for other Lepidoptera, family for Coleoptera, order for Orthoptera). Direct observations of herbivores were rare, so only damage patterns were used in analyses. The different types of damage recorded were used as proxy for insect insect richness on Piper. Hereafter, the term “insect richness” refers to the richness of herbivore damage patterns on plants.

Data analyses

The percentage of leaf tissue consumed by each insect taxon (based on damage patterns) on individual leaves was summed to determine the total percentage of herbivory on each leaf, and the variance in herbivory was calculated as the variance in leaf damage within each plant in a plot. Due to high mortality rates across sites, Piper interspecific richness in each plot was the final number of Piper species present in each plot at the end of the experiment (rather than the number of species planted) and was analyzed as a continuous covariate rather than a categorical treatment. Because most plant deaths occurred early in the experiments, final Piper interspecific richness more accurately reflects the local plant richness experienced by herbivores. In order to compare the magnitude of effect of intra and interspecific diversity on predictor variables, interspecific diversity for each plot was analyzed as the proportion of the interspecific diversity present in the species pools at each site. The effect of intra and interspecific diversity on each response variable is reported as the change in response as diversity changes from the lowest possible value at each site to the highest.

The effects of intraspecific richness, interspecific richness, and water addition on percent leaf area consumed, the percentage of damaged leaves, variance in herbivory, and insect insect richness were analyzed across and within sites using hierarchical Bayesian models (HBM). This framework acts as the Bayesian equivalent of a random-effects model where site is a random effect, allowing for generalized parameter estimates across sites. Analyses testing Piper intraspecific richness and interspecific richness were conducted using data from all sites; the effects of water addition and its interactions with richness were conducted using data from Ecuador, Costa Rica, and Peru. Due to high mortality in Peru, interactions between the water addition and intraspecific richness treatments could not be modeled for measures of herbivory or insect richness and were only tested for survival.

We acknowledge that treatment level combinations were not the same across the different sites, but that is the nature of random-effects models, where different levels of random factors, such as site or year, consist of treatment levels that are unique to that site or year. In fact, this type of design goes back to the origins of mixed models (40, 41), and are also a strength of the design, since lack of interactions between fixed and random effects increases generality (42). Even at a single site, manipulated variables in ecological experiments do not even consist of the same manipulations across all of the units of replication, as they suffer from problems such as multiple versions of treatments, interference, and noncompliance (43)

Bayesian structural equation models (BSEMs) were constructed for each site using all treatments as exogenous variables and insect richness and herbivory as endogenous variables. Three path models (Table S1) were developed under the assumptions that 1) intraspecific and interspecific richness may influence herbivory both directly and indirectly by modulating insect richness 2) water addition may influence both herbivory and insect richness and 3) insect richness may influence herbivory directly (Table S1). Models I, II, and III were tested in Ecuador, Costa Rica and Peru, while model III was used to analyze data from Mogi-Guaçu and Uaimii. Additional models incorporating interactions between intraspecific richness, interspecific richness, and water availability were tested in Costa Rica, Ecuador, and Peru.

Models were run at the leaf level for herbivory and insect richness, and at the plant level for the percentage of leaves with damage, variance of herbivory, and Piper survival. For all HBMs and BSEMs, model convergence was estimated visually using traceplots and an R^ discriminatory threshold of 1.1 (44). Model fit was determined via posterior predictive checks (PPC) using the χ2 discrepancy of the residuals as the discriminatory function (45). A PPC near 0.5 indicates a high model fit, while values near one or zero indicate poor fit. BSEM models were further compared using the Deviance Information Criterion (DIC).

All Bayesian models were written in JAGS via the jagsUI package in R (46) using Markov chain Monte Carlo (MCMC) samplers with weakly-informative priors. For the majority of models, MCMC runs were conducted for 10,000 iterations using the first 1,000 iterations as a burn-in phase to generate posterior distributions of parameter estimates for each response variable. HBMs modeling interactions required 20,000 iterations with the first 5,000 as burn-ins for all models to consistently converge. Mean parameter estimates and 95% credibility intervals (CIs) were calculated for all responses. 95% CIs which do not cross the y-axis are typically associated with less than a 2.5% type S error rate (47). A posterior probability of direction (PD) was calculated based on the percentage of the posterior distribution responding in the same direction as the median response. A PD of 95%, for example, indicates that the same direction of response (e.g. a positive or negative response) was observed in 95% of iterations, regardless of the magnitude of the response (48). PDs less than 95%, indicate lower confidence that a relationship exists, but can still be interpreted as the probability that an effect exists.

To analyze the effects of the experimental treatments on plant survival, survivorship curves were constructed in all sites. Analyses of Piper survival were based on the initially planted interspecific richness treatment of each plot because mortality occurred early in the experiment. The effect of water addition, intra- and interspecific richness and species identity on Piper survival were analyzed using Cox proportional hazard models. All data were analyzed using R software v4.1.0 (49).

Results

Considerable variation in herbivory and plant mortality was observed among study sites. Percent herbivory was lowest at Mogi; only 9% of leaf tissue was consumed by herbivores compared to 22% of leaf tissue consumed in Uaimii. Piper mortality was highest in Peru (89%), likely due to an El Niño event, and was lowest in Ecuador (27%; Fig. S3 A). Piper mortality resulted in the reduction of species richness in many plots and the loss of some combinations of treatments, notably the loss of all Piper in unwatered, low intraspecific richness plots in Peru (Fig. S1 B). Plots originally planted as high interspecific richness treatments had a final richness of between 1 and 12 species in Costa Rica, 9 and 11 species in Ecuador, 1 and 3 species in Peru, 2 and 4 species in Mogi, and only 2 species in Uaimii. We identified herbivore damage from 13 insect taxa at the five study sites, as well as damage by leaf miners of unknown orders. A total of 10 taxa were observed in Costa Rica, Ecuador, and Peru, 9 taxa were observed in Mogi, and 8 taxa were observed in Uaimii (Fig. S3 B). The majority of taxa in Costa Rica and Ecuador were Piper specialists, but the proportion of leaf tissue consumed by specialists was only greater than generalist damage in Costa Rica.

Our experiments revealed pronounced heterogeneity in ecosystem responses to changes in water and Piper diversity between sites (Figs. S4 & S5). PPCs for all HBMs were within 0.03 of 0.5, indicating models fit well. Model II (Table S1) was selected as the most parsimonious causal model for BSEMs in Costa Rica (PPC = 0.499), Ecuador (PPC = 0.499), and Peru (PPC = 0.498). Fit was high for model III in Mogi (PPC = 0.5) and Uaimii (PPC = 0.497). When all sites were analyzed together, percent herbivory was 4.2±3.6% (mean ± 95% CI) lower in plots that received additional water (PD = 98.7%), but we did not detect a cross site effect of water on the percentage of leaves with damage (Fig. 2 A-B). We did not detect any direct cross-site effects of water addition, insect richness, or plant diversity on the variance in herbivory (Fig. 2 C), or Piper survivorship (Fig. 2 D). Notably, greater Piper interspecific richness was associated with a 15 ± 18.6% increase in the richness of insect herbivores when comparing high diversity plots to low diversity plots across sites (PD = 95.0%; Fig. 2 E) and an indirect increase in herbivory was mediated by insect richness (Fig. 3). Insect richness was associated with an 8.8±2.8% increase in herbivory per insect taxon present (PD = 100%) and a 6.7±6.9% increase in the percentage of leaves with any damage (PP = 97%) (Fig. 2 A, B). Changes in both intraspecific richness and interspecific richness altered these ecosystem processes considerably, although effects varied in strength and direction across sites (Fig. 3, Table S4). Although intraspecific richness was not quantified, the genetic origin of Piper plants (clones vs. independent propagules) affected herbivory, plant mortality, and insect richness with effect sizes of similar or even greater magnitude than interspecific richness when standardized as the proportion of available species richenss at each site. However, we found that intra- and interspecific richness often had contrasting directions of effect on insect richness and multiple measures of herbivore pressure (Figs. 3, S4 B-C & S5 B-C). For example, in Costa Rica insect richness was 6.0 ± 5.9% (PD = 93.7%) lower in plots with high interspecific plant richness compared to low, while high intraspecific richness increased insect richness by 7.2 ± 6.9% (PD = 98.1%). In contrast, in Ecuador high interspecific richness was associated with a 43.6 ± 6.4% (PD = 100%) increase in insect richness, and high intraspecific richness decreased insect richness by 13.6 ± 5.2% (PD = 100%; Fig. S5 B-C).

Hierarchical Bayesian model parameter estimates for the effects of water availability intraspecific and interspecific Piper richness on herbivory, herbivorous insect richness, and Piper survivorship across sites Violins represent the cross-site posterior parameter distribution for each relationship in site level hierarchical Bayesian models.

Black lines represent the median posterior estimate and white bars represent 95% credibility intervals. Percentages above violins indicate the probability of an effect being positive or negative (as indicated by a negative probability) in response to an increase of the independent variable. Distributions for water addition compare watered and control plots; distributions for interspecific richness compare Piper species richness standardized as the proportion of the maximum richness used at a site; distributions for intraspecific richness compare low and high intraspecific richness treatments; distributions for insect richness compare responses per insect taxon present on an individual leaf

Direct and indirect effects of plant diversity and water availability on insect herbivores at five study sites.

Bayesian structural equation models comparing effects of different drivers of herbivorous insect richness and herbivory at five sites. Standardized path coefficients are means of the posterior distribution for the effects estimated at each causal path. Positive relationships are indicated in blue with triangular heads, and negative relationships are indicated in red with circular heads. Black arrows indicate path coefficients of zero magnitude. Dot plots summarize the standardized mean of the posterior distribution for each causal path with 95% and 80% credibility intervals. Asterisks indicate causal paths where the probability of an effect being positive or negative is >95%.

The effects of Piper richness were especially relevant when examining interactions between Piper richness and water addition (Fig. 4). The effects of water addition were altered by Piper intra- and interspecific richness at all sites. While a reduction in herbivory driven by water availability was observed in Costa Rica (4.7±2.5%, PD = 100%) and Peru (5.1±3.7%, PD = 100%), this effect was only present in Ecuador when intraspecific richness was high (Fig. S6). Across sites, water addition had negligible or negative effects on insect richness at low interspecific richness, but this effect was reduced or reversed when interspecific richness was high. Insect richness increased by 20.8±18.9% more in interspecifically diverse plots that received water compared to unwatered plots in Costa Rica (PD = 99%), 15.2±16.8% more in rich, watered plots in Ecuador (PD = 96%), and 29.3±31.5% more in Peru (PD = 98%) Additionally, Water addition had a negligible effect on insect richness in Ecuador when intraspecific richness was low, but increased insect richness by 12.1% when intraspecific richness was high (PD = 99%) (Fig. 4 B). Water addition had a negative effect on Piper survival in Costa Rica when intraspecific richness was low, but improved survival by 12.1 ± 8.5% in high intraspecific richness plots (PD = 99%; Fig. 4 C). Interactions with water in Peru may have been influenced by an El Niño event which resulted in abnormally low rainfall and high Piper mortality, while the typically wetter sites in Costa Rica and Ecuador experienced greater precipitation than average (Fig. S2). As there was 100% mortality in low intraspecific richness, unwatered plots in Peru, we could not explore interactions between intraspecific richness and water for any response except Piper mortality in this site.

Interactions between intraspecific or interspecific richness and water availability on herbivory, herbivorous insect richness, and Piper survival Bars indicate mean response and standard error of the mean.

Percentages indicate the probability that the two slopes are different as calculated using site level hierarchical Bayesian models. Due to high mortality in Peru, interactions between water and intraspecific richness could not be compared for any responses except mortality. Single and two species richness plots have been combined for visualization purposes only.

Cox analysis of plant survival

Water addition did not affect Piper survival in Costa Rica (z = -1.2, P > 0.2), or in Ecuador (z = 0.2, P > 0.8), but survival was reduced by 48% in plots without water in Peru (z = 3.2, P = 0.001). Intraspecific richness increased survival by 38% in Costa Rica (z = -3.3, P < 0.001) and by 32% in Ecuador (z = -4.9, P < 0.001), but had no effect on survival in Peru (z = - 1.55, P > 0.1), Mogi (z = -0.22, P > 0.8), or Uaimii (z = 0.18, P > 0.8). Interspecific richness had no effect on survival in Costa Rica (z = -1.5, P > 0.1), Ecuador (z = -0.56, P > 0.5), Peru (z = - 0.49, P > 0.6), or Uaimii (z = 0.21, P > 0.8). In Mogi, survival was reduced by 13% in plots with higher interspecific richness (z = -2.16, P = 0.031). There was an interaction between intraspecific richness and water addition in Costa Rica and Peru. Survival in Costa Rica increased in response to water in high intraspecific richness plots, and decreased in response to water in low intraspecific richness plots (z = 3.8, P < 0.001), while the opposite pattern was observed in Peru (z = -2.4, P = 0.02; Fig. S7). Although plant die-offs caused plot level species extirpation at several sites, plant species identity was not related to survival rates with the exception of P. peltatum, which had the lowest survival rate of any species planted in Costa Rica. Further statistical results are available in the supplement.

Discussion

Our results demonstrate two key patterns. First, the strength of effects of intraspecific richness on higher trophic levels is comparable to that of interspecific richness, supporting our predictions and corroborating recent studies demonstrating the importance of intraspecific richness (18, 19). However, the effects of intra- and interspecific richness on herbivory, insect richness, and plant survival can differ dramatically, in contrast to our predictions. Second, we found that perturbations in water availability can have complex effects on herbivores and plant survival, and that these effects can be modulated by plant diversity in a context sensitive manner. While our prediction that water availability would influence herbivory across sites was supported, our results suggest that biodiversity loss and climatic perturbations may have dramatically different effects on ecosystem function at local scales, which may diminish our ability to predict how local communities will change as anthropogenic stressors increase.

Although experimental methods varied somewhat between study sites, we do not believe this can fully explain the level of heterogeneity observed in ecosystem response. For example, the methods employed in Costa Rica and Ecuador were nearly identical, and yet the directions of effect of intra and interspecific diversity on insect richness were reversed in these sites (Fig. S5 B-C). The only meaningful differences between these two sites were the climate and the composition of the local Piper and insect communities. There was considerable variation in both biotic (Fig. S3 B) and abiotic factors (Table S2) across sites, which may have contributed to the heterogeneity observed in ecosystem response to Piper diversity.

As we did not directly measure plant stress or nutrition, it is difficult to determine the exact mechanism through which water addition reduced herbivory. The presence of an El Niño weather pattern during the course of the Peru experiment may have led the water addition treatment to relieve plants from drought stress, while water addition in Costa Rica and Ecuador may have added to water stress in treated plants as these sites received an above average level of rainfall during the course of the experiment. Despite this, water addition consistently suppressed herbivory in Costa Rica, Ecuador, and Peru under natural levels of Piper diversity, suggesting that predicted increases in precipitation in the next century (38) will dramatically alter the flow of resources from primary producers. Although we were only able to record the richness of insect herbivory patterns, this measure is indicative of the functional diversity of insect herbivores on Piper and any changes to this value represent changes in interactions between Piper and higher trophic levels (31). As such, the additional reduction in effects of water addition on insect richness when Piper richness was low suggest that biodiversity loss in tropical systems will alter the ability of higher trophic levels to respond to environmental perturbations.

As herbivore taxa can be differently affected by manipulations of diversity (50), variation in the direction of the effect of intra- and interspecific richness may be due in part to variation in the composition of insect communities and herbivore pressure measured across sites. Changes in neighborhood effects when small numbers of plant species dominate a community can lead either to the reduction or increase of herbivore pressure, dependent both on the nature of plant species lost and herbivore species present in a given community (51). As such, local variation in community composition has the potential to greatly alter the effects of both climate change and biodiversity loss on ecosystem function. Regardless of how biodiversity loss effects ecosystem function at large scales, variation in abiotic and biotic factors at locals scales can alter these effects, reducing our ability to predict how anthropogenic activity will alter ecosystem function. This is especially relevant in tropical systems, which have been the subject of far fewer studies of ecosystem function. As such, our knowledge of local scale effects on the relationship between biodiversity and ecosystem function remains limited in these systems.

A long-standing question in ecology has been the extent to which ecosystem function increases with biodiversity and if this relationship plateaus at a level past which ecological redundancy predominates. Recent results from less complex temperate grassland systems suggest that these ecosystems can be described by a mostly linear relationship between richness and function, where even rare species contribute uniquely to ecosystem function (52). In these systems, high contingency can be expected, where ecosystem-level effects will depend on most of the interacting species. In contrast, we might expect that diverse, tropical communities could be characterized by greater ecological redundancy and thus be subject to less contingency (53, 54). Despite these expectations, our results demonstrate heterogeneity in ecosystem response to changes in both intraspecific and interspecific richness in five tropical sites, suggesting that complexity in these systems may not reduce the contingency effects of biodiversity loss. Understanding the impact of biodiversity loss in tropical forests is fundamental to our ability to conserve those systems, and our findings highlight the importance of approaching the study of changes in ecosystem function as context sensitive responses in complex ecosystems.

Acknowledgements

The authors would like to thank Hacienda San Isidro, Ecuador for access to their forest reserve for experiments, the Meteorological Service from Peru (SENAMHI) for access to their weather data, and IEF, Minas Gerais, for access to Uaimii Forest Reserve during experiments. The collecting permit in Peru was issued by Peruvian Ministry of Agriculture (RDG N 288-2015-SERFOR-DGGSPFFS).

Funding sources

Funding was supplied through NSF grants DEB-1442103 to LAD, LAR, MLF, TLP, AMS, and DEB-1442075 to EJT

Author contributions

L.A.R., L.A.D., A.M.S., T.L.P., E.J.T., M.J.K., and T.J.M. developed the original idea for the research and secured funding. L.A.S., T.J.M, L.A.D, E.J.T., and H.G.L. collected Peru data.

H.G.L. and L.A.D. collected Costa Rica data. W.R.S., T.R.W., and L.A.D. collected Ecuador data. T.J.M and A.R.N collected Mogi data. Y.A., L.B., T.J.M, and A.R.N. collected Uaimii data.

W.R.S measured herbivory for all sites. A.J.G. analyzed the data and wrote the first draft of the manuscript. All authors contributed to analysis ideas and to writing revisions.

Competing interests

The authors declare no competing interests

Supplementary figures and tables

Path models and explanatory hypotheses.

Black arrows indicate causal paths between exogenous and endogenous variables.

Study site characteristics and experimental details.

Species of Piper used at each study location

Mean parameter estimates and probability of direction (PD) for the effects of increases in intraspecific diversity, interspecific richness, water availability and insect richness on measures of herbivory, plant mortality, and insect richness.

Treatments and number of plots used across sites

A) Treatments of intraspecific richness, interspecific richness, and water addition at the beginning of the experimental period at each of the five study sites, and B) final treatments at the end of the experimental period in each site. White numerals indicate number of plots used. Changes in Piper species richness and the loss of some treatment combinations was due to Piper mortality over the course of the experimental period.

Precipitation levels at study sites where the water addition treatment was applied

A) Monthly precipitation over the experimental period at Costa Rica (CR), Ecuador (EC), and Peru (PE). B) difference in monthly precipitation from climate normals at the three sites across the experimental periods for those sites. Bars indicate median values, black points indicate mean values. C) Cumulative precipitation over the course of the experimental period. Dotted lines indicate average cumulative precipitation, shaded regions indicate 95%, 80% and 50% quantiles. Climate normals are based on monthly precipitation for a period from 1958 to 1998.

Overall herbivory, plant mortality, and insect richness at five study sites

A) Percent herbivory, percentage of leaves with any damage, herbivorous insect richness, and percent mortality of Piper in Costa Rica (CR), Ecuador (EC), Peru (PE), Uaimii (UI), and Mogi-Guaçu (MG), across treatments. B) Proportion of herbivorous insect taxa observed at each site measured by feeding damage patterns. Colored regions indicate the percentage of damage observations contributed by each taxon. C) Proportion of herbivory attributed to each taxon at five study sites. Colored regions represent the percentage of each leaf consumed by each taxon at each site. White tiles represent sites where no leaf damage by that insect taxon was observed.

HBM parameter estimates of percent herbivory, percentage of leaves with damage, variance in herbivory, and percent Piper survival against levels of A) water addition, B) Piper intraspecific richness, C) Piper interspecific richness, and D) insect richness at each site.

Violins represent the posterior parameter distribution for each relationship across sites and within sites in Costa Rica (CR), Ecuador (EC), Peru (PE), Uaimii (UI), and Mogi-Guaçu (MG). Black lines represent the median posterior estimate while white bars represent 95% credibility intervals. Percentages above violins indicate the probability of an effect being entirely positive or entirely negative in response to an increase of the manipulated variable. Distributions for water addition compare unwatered control and watered plots; distributions for interspecific richness compare Piper species richness and are standardized as a percentage of the high richness treatment level; distributions for intraspecific richness compare low and high intraspecific richness plots; distributions for insect richness compare responses per insect taxon present at the leaf level

HBM posterior parameter estimates of insect richness compared to levels of A) water addition, B) Piper intraspecific richness, and C) Piper interspecific richness

Violins represent the posterior parameter distribution for each relationship across sites and within sites at Costa Rica (CR), Ecuador (EC), Peru (PE), Uaimii (UI), and Mogi-Guaçu (MG). Black lines represent the median posterior estimate while white bars represent 95% credibility intervals. Percentages above violins indicate the probability of an effect being entirely positive or entirely negative in response to an increase of the manipulated variable. Distributions for water addition compare the unwatered control and watered plots; distributions for interspecific richness compare Piper species richness and are standardized as a percentage of the high richness treatment level; distributions for intraspecific richness compare low and high intraspecific richness plots

Bayesian structural equation models for drivers of insect richness, herbivory and Piper survival at three sites, including interactions between intraspecific and interspecific richness, and water addition.

Path coefficients indicate the standardized mean of the posterior distribution for each causal path. Positive relationships are indicated in blue with triangular heads, and negative relationships are indicated in red with circular heads. Dot plots indicate the standardized mean of the posterior distribution for each causal path in A with 95% and 80% credibility intervals. Asterisks indicate causal paths where the probability of an effect being entirely positive or entirely negative is > 95%.

Percent Piper survival over time in five sites, compared to levels of intraspecific richness and water addition.