The relative importance of abiotic and biotic environmental conditions for taxonomic, phylogenetic, and functional diversity of spiders across spatial scales

Both abiotic and biotic conditions may be important for biodiversity. However, their relative importance may vary among different diversity dimensions as well as across spatial scales. Spiders (Araneae) offer an ecologically relevant system for evaluating variation in the relative strength abiotic and biotic biodiversity regulation. We quantified the relative importance of abiotic and biotic conditions for three diversity dimensions of spider communities quantified across two spatial scales. Spiders were surveyed along elevation gradients in northern Sweden. We focused our analysis on geomorphological and climatic conditions as well as vegetation characteristics, and quantified the relative importance of these conditions for the taxonomic, phylogenetic, and functional diversity of spider communities sampled across one intermediate (500 m) and one local (25 m) scale. There were stronger relationships among diversity dimensions at the local than the intermediate scale. There were also variation in the relative influence of abiotic and biotic conditions among diversity dimensions, but this variation was not consistent across spatial scales. Across both spatial scales, vegetation was related to all diversity dimensions whereas climate was important for phylogenetic and functional diversity. Our study does not fully support stronger abiotic regulation at coarser scales, and conversely stronger abiotic regulation at more local scales. Instead, our results indicate that community assembly is shaped by interactions between abiotic constrains in species distributions and biotic conditions, and that such interactions may be both scale and context dependent. Supplementary Information The online version contains supplementary material available at 10.1007/s00442-023-05383-0.


Appendix 1, Supplementary figures and tables
. Phylogeny used to quantify phylogenetic diversity of spider assemblages. The phylogeny was based on genetic sequences from the mitochondrial COI region. Figure S2. Functional dendrogram of spiders used to quantify functional diversity of spider assemblages. The dendrogram was based on a trait matrix consisting of six traits related to size, hunting behaviour, diet preferences and dispersal abilities.  Table S1. Characters used for the spider trait matrix, their data type, and the taxonomic accuracy, quantified as the percentage of taxa for which we had trait information at the taxonomic rank of species, genus, subfamily, family, and nearest family.

Sequence alignment and tree construction
We obtained sequences from GenBank (Sayers et al. 2021) in March of 2021 for all but eleven of the collected spider taxa (Agyneta similis, Diplocentria sp., Gnaphosa lapponum, Gnaphosa leporina, Mecynargus morulus, Mugiphantes whymperi, Oreoneta frigida, Oreoneta sinuosa, Oryphantes angulatus, Porrhomma campbelli, Zornella cultrigera), which lacked available sequences. Observe that the taxonomy of the genus Oreoneta is poorly investigated. For these eleven taxa, we generated sequences by extracting DNA from our own morphologically identified specimens. The DNA extraction is described in the section below. If we obtained multiple sequences, we first aligned them for each taxon separately to generate a consensus sequence to be used in the full alignment. Likewise, sequences for specimens identified to a higher taxonomic rank or belonging to one of two closely related species (Walckenaeria karpinskii/clavicornis) were obtained by first aligning all sequences from our identified species within that taxonomic rank (family/genus etc.) to one consensus sequence before the full alignment. All sequences were aligned using the software Geneious Prime (version 2021.1.1., Biomatters Ltd., Auckland, New Zealand) using a global alignment with free end gaps and a 65 % similarity cost matrix, gap open penalty "12", gap extension penalty "3" and two refinement iterations. This alignment method uses an iterative approach of pairwise alignments, described by Feng and Doolittle (1987). A phylogenetic tree was built using a maximum likelihood approach by running the alignment of DNA sequences in the software IQ Tree (Nguyen et al. 2015) with a combination of ultrafast bootstrap (Hoang et al. 2018) and SH-aLRT branch test (Guindon et al. 2010), both with 1000 iterations. Additionally, the IQ Tree software were let to automatically choose the best fitted substitution model for our data based on BIC (Bayesian Information Criterion).

DNA extraction
For taxa lacking publicly available COI sequences, we generated them using the following methods. DNA extraction was done by removing a couple of legs, depending on the size of the spider. If several specimens were available, one spider from each site was used as replicates. The spider legs were crushed in a tube to increase DNA yield prior to the addition of a lysis buffer containing proteinase K (either Thermo Scientific Cell and Tissue Kit lysis buffer or buffer from Marquina (2022)), after which the sample was incubated under rotation at 56ºC for 2.5 hours or overnight. After incubation the lysate was transferred to a plate and extracted using Thermo Scientific Kingfisher Flex extraction robot, following standard protocol with reagents from the Thermo Scientific Cell and Tissue Kit. The DNA extract was amplified using the COI primers BF3 and BR2 (Elbrecht and Leese 2017;Elbrecht et al. 2019) and using a PCR mix consisting of one Illustra Hot Start Mix RTG bead (GE Healthcare Life Sciences, Freiburg, Germany), 21µl lab grade water, 2µl DNA-extract, 1µl forward and 1µl reverse primer. The PCR were run using the following temperature protocol: initial Taq-activation and denaturalization phase, 95ºC for 5 minutes, followed by 30 cycles of 95ºC for 30 seconds (s), annealing phase at 48ºC for 45s followed by an extension phase of 68ºC for 45s and the program was ended by a final extension phase of 72ºC for 10 minutes. The PCR products were then checked on an agarose gel for amplification success and amplicon length. Negative controls were used throughout the lab work from lysis, extraction, and PCR but were discarded after confirming that they lacked sign of amplification. PCR products were cleaned using ExoSAP-IT (Bell 2008) and sent to Macrogen Inc (Amsterdam, the Netherlands) for sanger sequencing of both forward and reverse complimentary sequences.

Calculations of environmental characteristics
We extracted geomorphological and climate characteristics using 25x25 m square polygons centred on the central coordinate of each sample station. For the intermediate scale, we used the three sampling stations pooled whereas we used individual values for each station for the local scale. These characteristics were calculated per pixel and represented as averages across a total area of 1875 m2 for each transect (intermediate scale) and 625 m2 for each station (local scale). The geomorphological characteristics aspect, slope and TWI were calculated from a digital elevation model (DEM) with a pixel size of 2 m from the Swedish land survey (Lantmäteriet 2022). Aspect and slope were expressed in degrees and calculated following Horn's algorithm (1981) using a 3x3 pixel neighbourhood size. For aspect, the values indicate the azimuth of the slope with 0 degrees representing the geographic north, going clockwise. TWI describes the water accumulation potential across a landscape and was calculated as the natural logarithm of the ratio of an upslope catchment area and the tangent of the slope of the target pixel (Beven and Kirkby 1979). The catchment area was calculated using the multiple flow direction algorithm (Quinn et al. 1991), which has been suggested as appropriate for ecological studies (Kopecký and Čížková 2010). Bedrock silica content was based on bedrock information from the geological survey of Sweden (SGU 2021) and quantified as a factor with four levels based on standard classification of the percentage silica (SiO2), felsic (>65 %), intermediate (52-65 %, set to the average 58.5 %), mafic (45-52 %, set to the average 48.5 %) and ultramafic (<45 %).
Annual temperature and monthly temperature variation values were obtained from a temperature model with a 50 m resolution developed by Meineri and Hylander (2017). Values for monthly average precipitation and within year precipitation variation were obtained from WorldClim 2.1 data with a resolution of 30 arc seconds (~1 km2) as averages for the period 1970-2000 (Fick and Hijmans 2017). GIS analysis was done using QGIS Desktop (version 3.10.14, https://www.qgis.org).
We derived vegetation characteristics from our own surveys of 1 m2 sample plots placed at the same location as the pitfall traps ( Figure 1e). Within each of the sampling plots all vascular plant taxa were recorded to species or closest possible taxonomical rank, following the taxonomy by the database Dyntaxa (2021). The cover of alive vascular plants, mosses, lichens, and dead plant tissue was visually estimated to the nearest half percent using an 10x10cm intercept grid in the sampling plots. Lichen cover was defined as percentage of plot covered by non-crustaceous lichens. We estimated vascular plant diversity using the Shannon index (Shannon 1948), based on the relative abundance of each species in each plot. Relative abundances were estimated using an intercept method modified from the International Tundra Experiment (Naud et al. 2019). We inserted one intercept stick at 25 randomly placed points along the 10x10 cm grid in every plot. The same 25 points within the grid were kept for all plots. In addition, we measured the height of each plant touching the stick at the 25 intercept points. These measurements were used to calculate the variation in vascular plant height and providing information on the physical vegetation structure.