Richness and diversity
We draw rarefaction curves (function rarecurve , package ‘vegan’ (Oksanen et al., 2013)) to assess the sequencing depth of each sample for both COI and 16S datasets. We removed samples with a number of reads too low for downstream analyses (Figure S1). To assess the differences in MOTU richness we performed Generalised Linear Mixed Models (GLMMs) and Generalized Linear Models (GLMs), separately for the COI and 16S datasets. We could not pool all MOTUs in one single dataset since some of them could not be identified to the species level. All analyses were done with the datasets rarefied to the lowest number of reads per sample (function rrarefy , package ‘vegan’).
First, we fit Generalized Linear Mixed Models (function glmmTMB , package ‘glmmTMB’ (Brooks et al., 2019)) with Farm as random factor to estimate the effect of Treatment on the MOTU richness in all farms (N = 72) and we ran ANOVA tests with those models (functionAnova.glmmTMB , package ‘glmmTMB’). Then, we performed Generalized Linear Models (function glm ) with Farm and Treatment as independent factors, and richness as dependent variable, including only the experimental and control trees of Farms 1-3 (N = 48), followed by ANOVA tests of the model (function Anova ). We did so because it was in those farms in which there were both control trees and trees subjected to short-term livestock exclosure, so that a full model including both factors (Farm and Treatment) could be performed. All models were fit with binomial negative links. This distribution was chosen after comparing it with Poisson, Generalized Poisson and Conway-Maxwell Poisson distributions. Pairwise post-hoc analyses were made with Tukey tests. We repeated the same analyses at the Family and Order level, merging the COI and 16S datasets. In these cases, models were fit with Gaussian distribution.
We also analysed the alpha diversity of each plot by calculating the Shannon index, H’ (functiondiversity , package ‘vegan’). We fit new GLMM and GLM models following the previous scheme, with Gaussian distribution, and with diversity instead of richness as dependent variables. Model summaries, effects tables, and post hoc analyses can be found in Tables S8-10.