Community composition changes
To analyse the differences in sample composition at MOTU level among the treatments, we performed Non-metric Multidimensional Scaling (NMDS) ordinations of the trees (N = 72) based on the arthropod datasets of each marker. We used Bray–Curtis clustering to measure the dissimilarity between the trees (function metaMDS , package ‘vegan’). Posteriorly, we ran a Permutation Analysis of Variance (PERMANOVA) (functions adonis at package ‘vegan’) to analyse the compositional differences between samples from different treatments. The effects of Treatment and Vegetation structure on sample composition were explored by correlating them with the NMDS ordination (functionenvfit (McCune & Grace, 2002), package ‘vegan’), evaluating the strength of those relationships through the squared correlation coefficient (r2). Beta-diversity summaries can be found in Table S11.
We calculated the zeta diversity of each group of samples, that is, the number of species shared by a given number of trees (Hui & McGeoch, 2014), and also how the values of zeta diversity changed with increasing number of trees (McGeoch et al., 2019). We used the ratio of the zeta diversity decline to visualize the species retention rate (Figure S2): the probability of finding species shared by all trees as the number of trees increases. We first grouped all trees of the same treatment together (Figure S2A), and later we separated them by Treatment and Farm (Figure S2B). In both cases we followed the ALL model (the addition of trees does not follow any directionality, nor is dependent on proximity to one another), performing a Monte-Carlo chain sampling with 5,000 samples (function Zeta.decline.mc , package ‘zetadiv’ (Latombe et al., 2020)). The better fit to a power-law or exponential regression of the zeta diversity decline was decided on the basis of the Akaike information criterion.
We also visualized the most abundant MOTUs (those whose total abundance across all samples was above the third quartile of the abundance distribution) using a semiquantitative scale in which MOTUs with an abundance an order of magnitude greater appeared as double the size in the graph (Turon et al., 2022) (see Supplemental Information: Figures). To identify the MOTUs most responsible for the differences between treatments we calculated the Dufrene-Legendre indicator value d(Dufrene & Legendre, 1997) of the MOTUs for each treatment (functionindval , package ‘labdsv’ (Roberts, 2016)). Summary of the indicator species analysis can be found in Table S12. We visualized a subset of the datasets containing only these MOTUs using the semiquantitative scale mentioned above.