Tree species identity, canopy structure and prey availability differentially affect canopy spider diversity and trophic composition

Forest canopies maintain a high proportion of arthropod diversity. The drivers that structure these communities, however, are poorly understood. Therefore, integrative research connecting tree species identity and environmental stand properties with taxonomic and functional community composition of canopy arthropods is required. In this study, we investigated how the taxonomic, functional and trophic composition of arboreal spider communities is affected by tree species composition and associated differences in canopy structure and prey availability in temperate forests. We sampled canopy spiders as well as their potential prey using insecticidal fogging in monospecific and mixed stands of native European beech, native Norway spruce and non-native Douglas fir. Trophic metrics were obtained from stable isotope analysis and structural canopy properties were assessed with mobile laser scanning. Monospecific native spruce stands promoted local canopy spider abundance and diversity, but native beech and beech–conifer mixtures had the highest diversity at landscape scale. Spider community composition differed between monospecific stands, with broadleaf–conifer mixtures mitigating these differences. Irrespective of tree species identity, spider abundance, taxonomic diversity, functional richness and isotopic richness increased in structurally heterogeneous canopies with high prey abundances, but functional evenness and trophic divergence decreased. Our study shows that canopy spiders are differentially affected by tree species identity, canopy structure and prey availability. Broadleaf–conifer mixtures mitigated negative effects of (non-native) conifers, but positive mixture effects were only evident at the landscape scale. Structurally heterogeneous canopies promoted the dominance of only specific trait clusters. This indicates that intermediate heterogeneity might result in high stability of ecological communities. Supplementary Information The online version contains supplementary material available at 10.1007/s00442-023-05447-1.

The estimation of functional richness (FRic) is based on the total functional space occupied.Values between 0 and 1 represent the space filled by the community.Functional divergence (FDiv) measures the distance between species within the functional space.Values close to 0 indicate that the most abundant species have traits that are centered in the functional trait space, and values close to 1 represent dominance of traits which are distinct to each other.Functional evenness (FEve) quantifies the species value distribution in the functional space.If the values are evenly distributed, FEve is close to 1, whilst clustered values push FEve towards 0 (e.g. if most species have the same functional traits).

Stable isotope analysis
If multiple species with the same abundance added up to the top 80%, we selected the species with highest mean body mass according to Nentwig et al. (2021).Prior to stable isotope analysis, we discarded the opisthosoma and dried the spiders and tree leaves at 60° C for 72 h.The removal of the opisthosoma ensures the exclusion of isotopic bias due to recently digested material (Perkins et al. 2013).In case of large species, we bisected the prosoma longitudinally, whilst for very small species, we pooled multiple prosomas to reach the required minimum dry weight of ~ 0.05 mg.For leaf isotopic analyses, we pooled the tips and bases of two leaves per sampled tree.
All calculations of isotopic metrics are based on the convex hull area, covered by the isotopic ratios of the community (Layman et al. 2007).The interpretation is equivalent to the above described interpretation of functional indices.Further, isotopic uniqueness (IUni) reflects whether isotopic values within the community are unique or if they overlap.IUni approaches 1 when most isotopic values are unique, and IUni values close to 0 show high redundancy of isotopic values (e.g., if one trophic position is covered by multiple species).

Point cloud processing
Based on the laser scanning point cloud, we calculated the overall vegetation volume, the mean effective number of vertical canopy layers (ENL), mean horizontal canopy gap area (mean gap area), the coefficient of variation of horizontal canopy gap area per plot (CV gap area), the CV of intra-canopy gap height in the canopy per plot (CV ICG height) and the box-dimension as a measure of structural complexity.To calculate these variables, we cropped the point cloud to radii from 1-12 m, using the R package "LidR" (Roussel et al. 2020).We subsampled each point cloud using the space tool in CloudCompare, to remove points closer than 0.5 cm to each other.Thereafter, we normalized the point cloud with a digital terrain model as implemented in "LidR".After voxelizing the point clouds to voxels of 50 cm side length, we classified voxels as empty if they contained no point, and non-empty if they contained at least one point, with each point being a laser hit.
We considered non-empty voxels as vegetation and calculated the vegetation volume as the sum of all vegetation voxels times their volume.We calculated the ENL0 following the method of Ehbrecht et al. (2016).The ENL0 represents the number of non-empty voxels in the vertical voxel column.We calculated the ENL0 for each voxel column in the cropped point cloud and averaged it at plot level.For the calculation of horizontal gaps-interrupting the canopy-we considered contiguous cells with no vegetation above two meters as gaps (Brokaw 1982).Using the R package "ForestGapR" (Silva et al. 2019) we calculated the mean and CV gap area per plot.As canopies do not only feature horizontal gaps interrupting the canopy, but also gaps within the canopy, we adapted the quantification of horizontal gaps for the empty space within the canopy.For the calculation of these "intra-canopy gaps" (ICGs), we connected all neighboring empty voxels within the canopy (in horizontal and vertical directions) which had no connection to the empty space above or below the canopy to form 3D gaps.We calculated the CV ICG height per plot.As a measure of structural complexity, we calculated the box-dimension (Seidel 2018).To calculate this fractal index, we fitted the canopy point cloud in boxes of successively decreasing size, starting with a box the size of the bounding box of the point-cloud containing all points, up to a box of 10 cm size length (lower cut-off).For each box size, the number of boxes necessary to englobe all points are counted.The box dimension corresponds to the slope of the regression of the logarithm of the number of boxes as a function of the logarithm of the inverse of the box size.Further algorithmic details can be found in Arseniou et al. (2021).

Model selection and fitting
Using the variance inflation factor (VIF), we assured that no multicollinearity between fixed effects was given (VIFs < 5; Fox and Weisberg 2019).In the modelling step (ii), we reduced the models to the essential predictors based on the step function in the "lmerTest" R package (Kuznetsova et al. 2017).Because prey abundance had a strong positive skewedness, we log-transformed it to avoid disproportionally high impact of high values (Fink 2009).Due to non-normal residual distribution we log-transformed spider abundance and biomass values for analyses at the sample-level.

Table S 1
Stand age classes and coordinates of the study plots.In case of mixtures, age class 1 refers to the coniferous species and age class 2 refers to beech.Age class 3 is always referring to admixed age classes of beech.

Table S 2
Guilds and abundances of all identified and analysed spider species per stand type.Web weavers are split into orb web weavers, sheet web weavers and space web weavers.Hunting spiders are split into ambush hunters and other hunters.

Table S 3
Mean values of all analyzed structural properties per stand type (± SE).CV = coefficient of variation; ICG = intra-canopy gaps.

Table S 4
Model summaries of linear mixed-effects models of spider responses at the local scale (per collecting sheet) with stand type and vegetation volume as fixed effects.

Table S 5
Model summaries of linear models of spider responses at the plot-scale with stand type and vegetation volume as covariates.
Appendix S1: Table S 6Model summaries of linear mixed-effects models of spider responses at the local scale (per collecting sheet) with environmental variables (prey abundance, mean gap area, ENL0, CV gap area, CV ICG height, Boxdimension) as fixed effects and stand type and plot as random effects.Model reduction (step function) might have reduced the number of fixed effects depending on their predictive power.

Table S 7
Model summaries of linear mixed-effects models of spider responses at the plot-scale with environmental variables (prey abundance, mean gap area, ENL0, CV gap area, CV ICG height, Box-dimension) as fixed effects and stand type random effect.Model reduction (step function) reduced the number of fixed effects depending on their predictive power.

Table S 8
Pairwise comparison of the species compositions between stand types (pairwiseAdonis).

Table S 9
Post hoc fitted environmental variables on the NMDS axes.

Table S 10
Model summaries of linear models of functional spider responses at the plot-scale with stand type and vegetation volume as covariates.FRic = functional richness; Eve = evenness; Div = divergence.Model summaries of linear mixed-effects models of functional spider responses at the plot-scale with environmental variables (prey abundance, mean gap area, ENL0, CV gap area, CV ICG height, Box-dimension) as fixed effects and stand type random effect.Model reduction (step function) might have reduced the number of fixed effects depending on their predictive power.FRic = functional richness; Eve = evenness; Div = divergence.

Table S 12
Model summaries of linear models of isotopic spider responses at the plot-scale with stand type and vegetation volume as covariates.IRic = isotopic richness; Eve = evenness; Div = divergence; Uni = uniqueness.

Table S 13
Model summaries of linear mixed-effects models of isotopic spider responses at the plot-scale with environmental variables (prey abundance, mean gap area, ENL0, CV gap area, CV ICG height, Box-dimension) as fixed effects and stand type random effect.Model reduction (step function) might have reduced the number of fixed effects depending on their predictive power.IRic = isotopic richness; Eve = evenness; Div = divergence; Uni = uniqueness.