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.