Episodic herbivory, plant density dependence, and stimulation of aboveground plant production

Abstract Herbivory is a major energy transfer within ecosystems; an open question is under what circumstances it can stimulate aboveground seasonal primary production. Despite multiple field demonstrations, past theory considered herbivory as a continuous process and found stimulation of seasonal production to be unlikely. Here, we report a new theoretical model that explores the consequences of discrete herbivory events, or episodes, separated in time. We discovered that negative density (biomass) dependence of plant growth, such as might be expected from resource limitation of plant growth, favors stimulation of seasonal production by infrequent herbivory events under a wide range of herbivory intensities and maximum plant relative growth rates. Results converge to those of previous models under repeated, short‐interval herbivory, which generally reduces seasonal production. Model parameters were estimated with new and previous data from the Serengeti ecosystem. Patterns of observed frequent and large magnitude stimulated production in these data agreed generally with those predicted by the episodic herbivory model. The model thus may provide a new framework for evaluating the sustainability and impact of herbivory.

Here, we consider herbivory to be episodic that is an event in which some fraction of biomass is removed (herbivory intensity), and production is the biomass accumulation (growth) over some time interval both before and following that event. We argue that herbivory can be inherently episodic, particularly in terrestrial environments where herbivores consume a significant fraction of production.
Individual plants may not experience daily or even weekly loss of tissue; rather, plants often lose a significant fraction of their biomass over short time periods. For example, plants lose tissue as bites, sometimes several in succession (Belovsky, 1986;Hobbs, Gross, Shipley, Spalinger, & Wunder, 2003) to mammalian herbivores, and specialist insect herbivores may impose most of the consumption of their host plant in their final instar of development (Agrawal, 1998;Brown, Gange, Evans, & Storr, 1987). Episodes of herbivory can occur for both herbaceous and woody plants in terrestrial environments (Broadbent, Bork, Cooke, & Willms, 2019;Broadbent, Bork, & Willms, 2018;Call & St Clair, 2018;Herrero-Jauregui, Schmitz, & Pineda, 2016;Huttunen et al., 2013;Mudongo, Fynn, & Bonyongo, 2016;Teague, Dowhower, & Baker, 2016) as well as plants in aquatic systems (van Tussenbroek & Morales, 2017). Fisheries models suggest that pulsed harvests cannot stimulate yields (Braverman & Mamdani, 2008). However, this model includes no explicit feedback between harvest and population growth rate, as might occur where plant growth is resource-limited and herbivory affects access to resources. Such feedback appears to be important but no theoretical framework for episodic herbivory exists from which to formulate and test hypotheses about episodic herbivory effects on plants and the ensuing plant responses.
To explore the consequences of episodic herbivory, we modify a previous model of plant response to herbivory (Hilbert et al., 1981), based on an assumption of exponential plant growth, to consider effects of herbivory on density (biomass)-dependent growth by plants. We then connect this more general approach to explicit resource limitation of plants to explore how resource supply and availability might influence plant density dependence and the impacts of herbivores. We derive first a simple description of plant growth as resource-limited and therefore density-dependent, and then explore the conditions of herbivory intensity, frequency, plant maximum relative growth rate, biomass at the time of herbivory, and steady-state plant biomass in the absence of herbivory that influence whether herbivores stimulate plant production by herbivores.
We confront the episodic herbivory model's predictions with both new and published field data from the Serengeti ecosystem (McNaughton, 1985;Ritchie, 2014). Measurements of community (all species') biomass at monthly intervals in fenced plots were used to estimate relationships between relative growth rate, RGR, and biomass and to estimate maximum RGR across seven sites that differ in rainfall and soil nutrients. These data were used to parameterize the model of episodic herbivory, which was then tested with reported herbivory (grazing) intensity and grazed and ungrazed production at twenty sites across the Serengeti from the classic study by McNaughton (McNaughton, 1985).

| Herbivory and biomass dynamics
We begin with a classic logistic description of the density-dependent rate of change in plant biomass where S is biomass (mass/area), r is a maximum relative growth rate under ideal conditions (mass . mass −1. time −1 ), t is time, and S K is the steady-state biomass in the absence of herbivory. Mass-specific, or relative, growth rate, RGR is given as Note that all these parameters are measurable in the field and have ecological meaning, including r, which represents a plant trait.
(1) dS dt = S r 1−S S K (2) RGR = 1 S dS dt = r 1−S S K The value of r can be attributed to a single species in the case of a specialist herbivore on its host plant or to an agricultural monoculture grazed by livestock, or it can reflect a community-weighted mean trait from a multi-species assemblage (Muscarella & Uriarte, 2016). RGR is assumed to decline linearly with increasing biomass, presumably due to greater resource limitation at higher biomass ( Figure 1a). The linear assumption is an approximation that allows time-dependent trajectories of biomass to be solved analytically and thus generally. Field results reported here support this approximation.

| Episodic herbivory
From some initial biomass S 0 at the time of herbivory, biomass S(t) at a future time t can be found by solving the separable differential equation in Equation (1), yielding An herbivory event removes a proportion G (herbivory intensity) of biomass at initial value S 0 , yielding biomass following herbivory of S 0 (1-G) ( Figure 1a). Thus, we can determine biomass at a future time t following the herbivory event as Inspection of Equation 4 reveals that increasing herbivory intensity, G, has both positive and negative effects on the factors affecting productivity. A positive effect is that herbivory increases RGR Note that production is low if biomass is near S K due to the balance of growth with tissue loss due to resource limitation. This is commonly observed in forests where biomass remains approximately constant over the season following initial leaf-out in the Graphical representation of the theoretical episodic herbivory model. (a) Under density-dependent plant dynamics, relative growth rate RGR (solid line) declines from a maximum ate r at S = 1, to zero at steady-state biomass without herbivory, S K . Biomass removal as a proportion G (herbivory intensity) by herbivores from an initial biomass S 0 at the time of the herbivory event to yield a biomass immediately following herbivory S 0 (1-G). The expected increase in RGR from that at initial biomass S 0 (RGR 0 ) to that at biomass following herbivory S 0 (1-G), RGR g , reflects the reduction in density-dependent effect of biomass on RGR. (b) Biomass accumulation over a time t without herbivory (solid curve) and following an herbivory event that removes a proportion G of biomass (intensity) (dashed curve). Seasonal production is estimated from the difference between biomass at time t without herbivory [P U = S(t) -S 0 ] and following herbivory [P g = S g (t) -S 0 (1-G)]. (c) Calculation of the difference in seasonal production over time t between plants following the herbivory event (dashed curve) and plants without herbivory (solid horizontal line) absence of herbivory (Xiao et al., 2004). Production under herbivory, P g , is therefore Production under herbivory (Equation 6) is a nonlinear function of herbivory intensity G (Figure 1c). Light herbivory intensity (small G) does not reduce biomass enough to achieve significant release from density-dependent biomass inhibition and increase RGR. If r is high enough, and t long enough, intermediate G can lead to higher productivity than in the absence of herbivory. However, under intense herbivory (large G) the reduction in biomass can be large enough that the absolute increase in biomass is reduced, despite a higher RGR, resulting in herbivore-induced decline in production.

| Field Measurements
As an example, empirical data measurements to estimate model parameters, such as r, t, G, and S 0 /S K , were obtained from two different studies in the Serengeti ecosystem, including Serengeti National Park in Tanzania

| Community biomass relative growth rate
Relative growth rate of community biomass, or that of all species, as a function of biomass was determined from biomass measurements made once, at least 27 days apart, in each of consecutive months April, May, and June in 2000, 2001 in three 4 × 4 m fenced plots at each of seven sites from the long-term grazing exclosure (LTGE) experiment (Ritchie, 2014;Veldhuis et al., 2019). Four 25 × 25 cm quadrats were clipped to the ground surface in each 4 × 4 m plot, with green material separated from litter. Different quadrats were clipped from different locations within plots in each month. Green material was air-dried at 45°C for a week and then weighed. Sites varied in mean annual rainfall from 490 to 890 mm/ yr, soil N from 0.05% to 0.22% and soil P from 0.005% to 0.15% (Ritchie, 2014).
RGR is operationally defined as ln (S m /S m-1 )/t where S is biomass, m is month, and t is the number of days between biomass samples.
Maximum community biomass relative growth rate, r, was estimated as the intercept of a regression of RGR versus S m-1 calculated for all twelve paired monthly samples (April to May, May to June, three plots, and two years). Points are not independent because the independent variable S m-1 is also present in the calculation of the dependent variable RGR. However, the regressions, conducted separately for each of the seven sites in the LTGE experiment, establish whether plant growth was density-dependent (negative slope) (Rees, Condit, Crawley, Pacala, & Tilman, 2001;Umana et al., 2018).
Herbivore effects on production as a function of herbivory intensity, G, as driven by (a) initial relative biomass S 0 /S K , (b) time interval following herbivory t, and (c) maximum relative growth rate or r. The horizontal dashed line represents equal productivity

| Herbivory and community biomass production
The impact of herbivory on production in the Serengeti was measured during 1974-1976 in the classic study by Sam McNaughton (McNaughton, 1985). Production estimates were made at twenty different sites, many of which are within 2 km of the LTGE experimental sites used to estimate maximum RGR, or r. Production was estimated using the moving cage method, which uses temporary, moveable exclosures to estimate biomass production over monthly increments in grazed grassland and permanent exclosures to estimate production in ungrazed control plots (McNaughton et al., 1996). Initial biomass was interpreted as S 0 (1-G) in grazed plots and as S 0 in permanently fenced plots. Subsequent biomass measurements at time t = 20-40 days later were interpreted as S g (t) in temporary fences erected at time t = 0, and as S(t) in the permanent exclosures at time t.
Production with and without grazing was estimated using Equations 9 and 10 and ΔP from Equation 11. Herbivory intensity G was measured as 1 -[(biomass in presence of herbivore)/(biomass in permanent exclosures)] at time t = 0, per McNaughton (1985). The original biomass and production data were not available, so values were read from graphs in the published paper and are shown in Table 1.
Analysis of patterns in the difference in aboveground monthly production in the presence and absence of herbivory, ΔP, was conducted with generalized linear models in SPSS 24. We used a model to explain ΔP that included hypothesized independent variables S 0 /S K , G, and G 2 (to capture nonlinear response of ΔP to herbivory intensity).

| Thresholds for stimulated production
Solving for the conditions under which ΔP > 0 in Equation 7 yields a threshold herbivory intensity, G', below which herbivory stimulates production: The relative biomass S 0 /S K is a proportion that reflects how close initial biomass S 0 is to steady-state biomass in the absence of herbivory S K . At higher S 0 /S K , the threshold herbivory intensity G increases and the magnitude of increase in ΔP also increases ( Figure 3) because herbivory more strongly weakens density dependence and increases RGR. Conversely, herbivory events that occur at low S 0 are more likely to leave remaining biomass at too low a biomass to yield large absolute increases in biomass during growth following herbivory and thus are unlikely to stimulate production.
G' also increases with time following the herbivory event, t ( Figure 2b), implying that longer rest following herbivory yields a greater chance of overproduction. Defoliated plants that have more time for their higher RGR to be in effect effectively "catch up" to the biomass plants would grow to in the absence of herbivory. Conversely, if herbivory events are frequent, with shorter rest times, there is insufficient time for higher RGR to produce absolutely more biomass.
Herbivory frequency, a largely underappreciated aspect of herbivory, therefore should strongly affect whether overproduction occurs.
Maximum RGR or r strongly influences overproduction.
Intuitively, higher r should make overproduction more likely because RGR is increased more strongly for a given reduction in plant biomass by herbivores if r is higher. Figure 3b shows that higher r yields a greater range of grazing intensities and initial biomasses under which overproduction can occur, and the magnitude of overproduction also increases with r ( Figure 2c). Equation 8 can be rearranged to yield a threshold value for r, r', as a function of S 0 /S K , G, and t that must be exceeded for overproduction to occur (Figure 4).  Figure 2 highlights the phenomenon of "herbivory (grazing) optimization," that is, an herbivory intensity exists that maximizes production above that expected in the absence of herbivory. Factors that increase the likelihood of overproduction also increase the herbivory intensity that yields maximum production. For longer times (t > 20)

| Optimized production under herbivory
between herbivory events, optimal herbivory intensity is more than 50% and the increase in production relative to that in the absence of herbivory approaches 50%. This magnitude of difference in production is also obtained when initial biomass approaches 70%-90% of the steady-state biomass without herbivory S K (S 0 /S K = 0.7-0.9). Despite the potential for strong overcompensation under such conditions, the model also predicts precipitous drops in productivity with increased herbivory intensities just beyond optimal values (Figures 1c and 2).

| Repeated herbivory events
The framework presented thus far allows consideration of the steady state of plant biomass S g (t) Equation 4 under repeated herbivory events of intensity G separated by a time interval t for plants at a given (10)

F I G U R E 3
Thresholds of herbivory intensity G' below which a grazing event will increase productivity over a time t. (a) G' in response to increasing initial relative biomass S 0 /S K , where overproduction occurs at all G below the curve. Effects on G' of (b) increasing time interval t following the herbivory, and (c) increasing maximum relative growth rate, or r Analysis reveals that above a certain G, G max = 1-e −rt , repeated herbivory impacts will reduce biomass to zero, and thus such intensity is not sustainable (Figure 5a). More interestingly, repeated herbivory intensities that can produce sustained overproduction occur under a more restricted set of conditions (Figure 5b). The critical threshold of G* below which overproduction occurs when at steady state during repeated herbivory events is obtained by substituting Equation 6 for S 0 /S K in Equation 8 and solving for G Analysis reveals that G* is an increasing function of both r and t, the time interval between herbivory events (Figure 5c).

| Testing with Field Data
Regressions of RGR versus S m-1 were significantly negative at all seven LTGE sites sampled in 2000-2002 ( Figure 6). All intercepts were significantly greater than zero and slopes significantly negative, as expected under moderate to strong density dependence in biomass production (Table 2). Intercepts represent estimates of r, and the mean intercept (± SE) across the seven sites was 0.049 ± 0.005 (N = 7). In

| D ISCUSS I ON
We present here a framework for considering both herbivory impact and time intervals between impacts. This departs from most previous models of plant-herbivore dynamics (Loreau, 2010;de Mazancourt et al., 1998;Ritchie, 2014;Schwinning & Parsons, 1999) where herbivory is assumed to be a continuous process, subtracting biomass on small time scales (dt). Continuous herbivory may not capture some more realistic dynamics of plants with their herbivores.

| Herbivory and plant production
Our analysis highlights several key predictions that differ from previous evaluations of herbivore effects on plant production. First, shortterm, highly intense herbivory can lead to overproduction as long as herbivory intervals are sufficiently long (Figures 2-4). Conversely, frequent but less intense herbivory is unlikely to stimulate productivity ( Figure 2). Plant traits that contribute to higher RGR at low (11) G * = 1 − e −rt 1 + e −rt F I G U R E 4 Critical thresholds of r for overproduction (Equation 13) as a function of herbivory intensity. (a) Plants with r above the threshold will overproduce following herbivory. (b) Biomass relative to that at steady state without herbivory, S 0 /S K and (c) time interval following herbivory, t biomass and thus higher maximum RGR makes overproduction more likely. Finally, overproduction is more likely when herbivory events occur at plant biomass closer to steady-state biomass in the absence of herbivory. Consequently, the timing of herbivory during the season may be an important factor in determining whether overproduction occurs.
Alteration by herbivory of plant traits that affect r (Adler, Milchunas, Sala, Burke, & Lauenroth, 2005;Diaz et al., 2007;Koricheva & Nykanen, 2004;de Mazancourt et al., 1998), or enhancement of resource availability by herbivores (de Mazancourt et al., 1998;McNaughton et al., 1997;Pastor & Naiman, 1992;Ritchie et al., 1998) are hypothesized mechanisms by which herbivory might stimulate seasonal production. However, a key result of our analysis is that herbivore stimulation of production does not require either of these mechanisms. Indeed r in the model is assumed to be constant with or without herbivory, so any increase in production with herbivory does not depend on an effect of herbivory on the inherent capacity for plant growth. Likewise, there is no assumption that plants use stored resources to regrow following herbivory, as biomass accumulation following the herbivory event depends only on growth associated with residual aboveground biomass. Similarly, there is no requirement that herbivory alters carbon allocation in plants, such as might reduce root biomass to compensate for lost aboveground tissue (Milchunas & Lauenroth, 1993). This prediction is supported by multiple studies that find either no or weak impacts of intense aboveground herbivory on root biomass and production (McNaughton, Banyikwa, & McNaughton, 1998;Milchunas & Lauenroth, 1993;Ritchie, 2014).
The episodic herbivory model presented here (Equation 7) suggests that not only is overproduction possible but also it may be expected in situations where herbivory is reasonably intense and initial plant biomass is close to S K . What is required for significant stimulation of productivity is strong density dependence, or reduction in RGR with increasing biomass, herbivory occurring when biomass is near S K , and long time intervals between herbivory events. This result provides an important alternative prediction to the resource limitation hypothesis (Wise & Abrahamson, 2005 and may explain why compensation is sometimes observed to be stronger in resource-poor conditions (Luo et al., 2012;de Mazancourt et al., 1998;Schwinning & Parsons, 1999;Yamauchi & Yamamura, 2004). Recent simulated herbivory studies of both temperate and tropical grass species support these hypotheses (Broadbent et al., 2018(Broadbent et al., , 2019Mudongo et al., 2016); higher intensity clipping of multiple grass species followed by longer rest periods increased production, while less intense and more frequent defoliation generally either had no effect on or reduced production.
The results of the 1972-1975 Serengeti field study of the effects of migratory grazing mammals on production were consistent with the predictions of the episodic herbivory model. The higher the initial biomass relative to peak or steady-state biomass without herbivory, the larger the magnitude of overproduction due to herbivory (Figure 6a). Likewise, stimulation of productivity was stronger at intermediate herbivory intensities (Figure 6b). The episodic herbivory model also predicted that overproduction may occur across a wide range of herbivory intensity and initial biomass conditions if r > .05 ( Figure 4); estimates of daily mean

| Time intervals between herbivory episodes
Long time intervals between intense herbivory events can occur in a variety of contexts. Herbivores can be highly aggregated in space through a variety of mechanisms. Specialist insects can exhibit a high degree of spatial aggregation (Berenbaum & Isman, 1989;Dennis, Young, & Gordon, 1998;Hassell, Comins, & May, 1991;Hunter & Price, 1998) within a population of hosts, and a given set of hosts can experience high variability in herbivore densities over time (Berryman et al., 1987). Even for small generalist herbivores that show population cycles, such as small mammals (Laine & Henttonen, 1983;Pimm, 1991), intense herbivory can be separated by years of little herbivory. Any of these scenarios can result in long "rest" from herbivory during times when specialist populations are low or concentrated in space on other patches of plants. Likewise herbivores that travel in herds typically present locally high densities that may impose intense herbivory over short periods followed by much longer periods of recovery.
In contrast to these scenarios, frequent, intense herbivory, which may maintain biomass at levels far below S K , is unlikely to However, insufficient time is available before the next herbivory event for biomass to accumulate at a rate higher than that of plants without herbivory. This is because seasonal production depends on both available biomass and RGR, and under frequent herbivory total plant biomass following herbivory events may be insufficient to yield high biomass accumulation over a short recovery time t. Frequent herbivory would be characteristic of systems of high densities of stationary herbivores, such as grazing lawns or shrubs grazed or browsed by resident wild mammalian herbivores (Archibald, Bond, Stock, & Fairbanks, 2005;McNaughton, 1985;Veldhuis, Hulshof, Fokkema, Berg, & Olff, 2016), grazing by diverse fish assemblages in coral reefs (Diaz-Pulido & McCook, 2003;Hixon, 1996), or continuously grazed livestock systems (Milchunas & Lauenroth, 1993).

| Implications for conservation and management
The theoretical results presented here may be highly relevant to understanding and managing grazing systems (Briske, 1993;Noy-Meir, 1993;Zegler et al., 2018). They suggest that systems with migratory grazers with herding behavior, such as occur in the Serengeti (Figure 6), Yellowstone (Frank et al., 2002(Frank et al., , 2016, and other natural ecosystems, may be more likely to enhance productivity than systems with "continuous" herbivory by sedentary consumers. Herding increases the short-term intensity but shortens F I G U R E 7 Patterns of herbivoryinduced difference between biomass production in grazed and biomass production in ungrazed (ΔProduction > 0; Equation 11 for N = 20 sites in the greater Serengeti ecosystem. (a) ΔProduction exhibits a significant (p = .03) association with initial relative biomass S 0 /S K . (b) ΔProduction is associated with measured grazing intensity. Statistical justification of a nonlinear function is provided in Table 3.
(c) Association of ΔProduction predicted by the episodic herbivory model, using parameters r = 0.06, t = 30 days, and initial relative biomass (S 0 /S K ) and grazing intensity (G) from each of the 20 sites, and observed ΔProductivity for the same sites. The dashed line is the 1:1 relationship  [1972][1973][1974][1975] the length of the herbivory event, and migration may increase the time interval between events. Such conclusions are consistent with studies of mammalian grazing by migratory ungulates in natural ecosystems where migrations are still present (Frank et al., 2002(Frank et al., , 2016de Mazancourt, Loreau, & Abbadie, 1999;McNaughton, 1985;Stewart, Bowyer, Ruess, Dick, & Kie, 2006), but see Knapp et al. (2012).
Experimental studies of "rotational" grazing, which feature variable durations of both herbivory events and rest in fenced paddocks, have yielded mixed results (Briske et al., 2008;Teague, Provenza, Kreuter, Steffens, & Barnes, 2013;Teague et al., 2011 (Broadbent et al., 2018(Broadbent et al., , 2019Mudongo et al., 2016;Teague et al., 2013). Previous studies generally lack the necessary data to determine the intensity and time intervals during grazing episodes, so it is not possible to test whether these published rest-rotation results contradict the predictions of our model. Regardless, the theoretical framework presented here provides a basis for measuring biomass and defining grazing events that could easily be tested in the field (Frank et al., 2016).
Our key result is that herbivory can lead to stimulated plant production when plant growth is resource-limited and thus density-dependent, even without herbivore enhancement of plant growth potential (r) or nutrient availability. This simple framework of herbivory events punctuated by periods of "rest" for plants that are otherwise resource-limited can be expanded in future work to consider different herbivory contexts in more detail, such as insectplant host dynamics, or grazing and browsing by large mammals, or herbivory on woody versus herbaceous plants. Likewise, the model could be modified to evaluate the impact of potential herbivore-induced modifications of plant traits and plant resources. This episodic herbivory framework provides new hypotheses and the basis for field measurements that could help understand the consequences to both agro-ecosystems and natural ecosystems of modifications to herbivory, including loss of migration corridors, introduction of exotic herbivores, management of domestic herbivores, and loss of large herbivores.

ACK N OWLED G EM ENTS
We thank Doug Frank and Jason Fridley for comments. The work was supported by US NSF grant DEB1557085 and National Park Service Cooperative Agreement P15AC01545.

CO N FLI C T O F I NTE R E S T
The authors have no competing interests.

This article has been awarded Open Data and Open Materials
Badges. All materials and data are publicly accessible via the Open Science Framework at https://doi.org/10.5061/dryad.qrfj6 q5bj.