Defaunation model
Using population matrices from the 64 tree species, we simulated the
effects of defaunation by altering particular matrix elements.
Specifically, we used the following parameters to estimate defaunation
impacts: (1) the change in seed predation rate after defaunation
(p ), (2) the proportion of seeds dispersed before
(dbase ) and after
(ddefaun ) defaunation, (3) the difference in
survival between dispersed and undispersed seeds (a ), (4) the
difference in survival between dispersed and undispersed seedlings
(b ), (5) the strength of conspecific negative density dependence
(β ), and (6) the change in trampling-induced seedling mortality
after defaunation (t ). For p , we compiled a list of values
based on published meta-analyses (Kurten 2013; Gardner et al.2019) and our own literature search (Table S4). For other parameters
(dbase , ddefaun , a ,b , and β ), we compiled lists of values for based on
published meta-analyses (Kurten 2013; Comita et al. 2014; Lamannaet al. 2017). Because research that quantified the effects of
reduced seedling trampling was limited, we compiled a list of values fort based on published studies (Roldán & Simonetti 2001; Rosinet al. 2017). Table S5 provides a full list of model parameters
and the sources of the data that we used to estimate them.
We ran 10,000 iterations of the model, randomly drawing values in each
iteration for p , dbase ,ddefaun , a , b , β , andt from lists of parameter values (see Table S5). Within each
iteration, we assembled a ‘baseline’ population matrix for each tree
species from the chosen parameter values. To make the baseline
population matrix, we used the original matrices for each species, but
we divided both seeds and seedlings into dispersed versus undispersed
stages (Figure 1). For each species, we used original vital rates and
the parameters dbase , a , and b , to
calculate the numbers of seeds dispersed or undispersed, the baseline
survival rates for dispersed and undispersed seeds, and the maximum
survival rates for dispersed and undispersed seedlings (Equations
S1-S8). Seedling survival can depend on the density of adult trees
(A ) due to Janzen-Connell effects (Zhu et al. 2015); so we
incorporated density dependence using a Ricker function (Ricker 1954):
sl = αeβA (1)
where sl is the survival of seedlings (dispersed
or undispersed), α is the maximum survival of seedlings
(dispersed or undispersed) in the absence of adult trees, and βis the strength of conspecific negative density-dependence.
For each species in each iteration, we also created a ‘defaunation’
population matrix (Figure 1). The defaunation matrices had three key
differences from the baseline matrices. First, the probability of seed
dispersal for ‘large vertebrate dispersed’ tree species was changed toddefaun . For other species, the probability of
dispersal remained set to dbase . Second, we usedp , the change in seed predator-induced mortality between baseline
and defaunation scenarios, to alter the survival of dispersed and
undispersed seeds (Equations S11-S12). Third, we used t , the
change in trampling-induced seedling mortality between the baseline and
defaunation scenarios, to alter the maximum survival (before accounting
for density-dependence) of dispersed and undispersed seedlings
(Equations S13-S15). With these new maximum survival values, actual
seedling survival was calculated using the same density-dependent
function (Equation 1).
Once the baseline and defaunation matrices of each species had been
built, we simulated 120 years of population growth under baseline and
defaunation scenarios. Each simulation began with an initial density of
25 adults ha-1, based on typical values assembled by
LaManna et al. (Lamanna et al. 2017), while the rest of the
initial population vector was calculated based on the proportions of
adults to other life stages at stable stage distribution. We then
calculated the ‘defaunation response’ as the final number of adults in
the defaunation population divided by the final number of adults in the
baseline population.
For each of the 10,000 interactions, we drew new values for our
defaunation parameters (p , dbase ,ddefaun , a , b , β , andt ) to set a new defaunation scenario, simulated baseline and
defaunation population dynamics for each species, and calculated
defaunation response of each species for that scenario.