Small to large-scale patterns of ground-dwelling spider (Araneae) diversity across northern Canada

This study was performed as part of the Northern Biodiversity Program, a collaborative research initiative created to document spatial and temporal arthropod diversity patterns in northern Canada (Ernst and Buddle 2013; Timms et al. 2013). This study was carried out in northern Canada at 12 sites encompassing three broad ecoclimatic regions: north boreal, subarctic, and arctic (see Ernst and Buddle 2015 for a map of study sites). We selected sites to maximize longitudinal and latitudinal coverage in each region, characterized by unique climate, soil type, and vegetation composition (Strong and Zoltai 1989). The hierarchical nested sampling design consisted of three spatial scales of diversity: sample (i.e., trap and grid), site, and regional (i.e., ecoclimatic region). In each ecoclimatic region, four sites were selected for sampling. At each site, we sampled ground-dwelling arthropods in mesic and wet natural open habitat (i.e., no canopy cover by trees) to capture a greater diversity of arthropods (Ernst et al. 2016). Mesic open sites were characterized by elevated topography and well-drained soils. The vegetation on mesic sites was a discontinuous cover of dwarf shrubs (e.g., willows (Salix spp.) and Labrador tea (Ledum spp.)) and perennial forbs like Dryas, saxifrages (Saxifraga spp.), or wood rush (Luzula spp.), with a layer of fructose lichen (Cladonia and Cladinas spp.). Open wet habitats were water saturated, with a continuous cover of moss (including sphagnum), saxifrages, and sedge-like cotton grass (Eriophorum spp.).

Each habitat (mesic and open wet) was sampled in three locations (grids) separated by at least 500 m for a total of six grids per site. Each grid consisted of 18 traps, with nine pitfall traps and nine yellow pan traps placed randomly, 15 m apart. We used pitfall traps because they collect ground-dwelling arthropods in a reliable and repeatable manner even at large spatial scales (Gotelli and Ellison 2002; Bowden and Buddle 2010b). The pitfall trap was a plastic circular container 10 cm in diameter, with an inner sampling cup (7 cm depth), covered by a square plastic corrugate to provide protection from rain. Each pan trap was bright yellow, 20 cm wide and 3 cm deep. All traps were sunk into the ground with the rim at the surface level. Propylene glycol diluted with water (2:1) was used as preservative with a small amount of surfactant to break surface tension. We emptied the traps approximately every four days. A total of 1296 traps (108 traps per site) were used over 30 degrees of latitude and 80 degrees of longitude across northern Canada. Each site was sampled for two weeks in 2010 or in 2011 during the peak of spider abundance (Table S1).

We identified all mature spiders, both males and females, to species or, in some cases where a species name could not be applied, to morphospecies (Dondale and Redner 1978; Dondale and Redner 1982; Dondale and Redner 1990; Platnick and Dondale 1992; Paquin and Dupérré 2003; Dondale et al. 2004). Nomenclature followed the Word Spider Catalog (2018), and voucher specimens were deposited at McGill’s Lyman Entomological Museum (Sainte-Anne-De-Bellevue, Québec, Canada). Data are available from the Dryad digital repository (Loboda and Buddle 2018).

Mature specimens were identified and pooled by trap for the two-week sampling period for each site. To calculate the completeness of the inventories, we estimated the absolute species richness in each ecoclimatic region using the nonparametric estimators Jackknife 1 and Chao 1. Completeness values were derived from the ratio between the observed and Jackknife 1 estimated richness. Individual-based rarefaction curves were used to compare the observed and estimated richness for each ecoclimatic region (Gotelli and Ellison 2002; Buddle et al. 2005).

We used ordination methods to analyze the assemblage composition. Prior to ordination analysis, singletons were removed to better detect overall composition similarities and to remove the potential effect of differential sampling effort between sites. Moreover, species abundance was standardized as the number of individuals caught per day per sampling period. We performed a non-metric multidimensional scaling ordination (NMDS) with two dimensions based on Bray–Curtis distance to determine species composition similarities where each point is a composition within a grid to visualize patterns at site and regional scales. NMDS was conducted using the metaMDS function in the vegan software package, which uses several random starting configurations to select a solution with low stress (Oksanen et al. 2016).

We used the Partition 3.0 software (Veech and Crist 2009) to additively decompose the total diversity of spiders into local diversity (the average diversity within (α1) and among (β1) traps and the diversity among grids (β2)), regional diversity (average diversity among sites (β3)) and northern diversity (average diversity among ecoclimatic regions (β4)). Partitioning diversity can be summarized as:

For each level, observed α and β diversities were compared with predicted α and β diversities obtained with a random distribution generated using 10 000 randomizations (Crist et al. 2003). We evaluated differences in diversity partitioning between regions by doing the analyses separately for each region.

We used variance partitioning to determine the relative influence of multiple complementary sets of environmental factors on spider composition (Peres-Neto et al. 2006). Three explanatory matrices with vegetation, climate, and spatial variables were tested. The effect of one matrix independent of the influence of others (marginal effect) and the conditional variation explained by the combined effect of several matrices were calculated. The species matrix was then transformed to the Hellinger distance, the most suitable transformation because of the high environmental heterogeneity between sites, reflected in the species matrix by a high proportion of zeros (Legendre and Gallagher 2001). A direct gradient analysis (redundancy analysis (RDA)) was performed on all data as well as separately for each region to determine which variables best explain the variation in spider composition. We used the permutation forward selection procedure to reduce the number of explanatory variables for the spider species matrix. We retained significant variables with a cumulative R2 adjusted below the R2 adjusted obtained with all variables. Monte Carlo tests with 10 000 permutations were used to determine the significance of axes of the RDA. Statistical analyses, other than additive partitioning, were done using the R software (R Core Team 2018).

At similar latitudes, species richness was generally higher in western Canada (west of the 96° meridian) than eastern Canada (east of the 96° meridian). For example, the Norman Wells site (NW) in western Canada had a species richness of 132 at 65° latitude, whereas Iqaluit (IQ) at 62° latitude in eastern Canada had a richness of only 22 species (Fig. 2). With 169 species (55% of the total species richness), the Linyphiidae family was the most diverse family of spiders collected (Fig. 3). The most abundant family was the Lycosidae, with 15 778 individuals and 32 species, representing 67% of the total abundance (Fig. 2). The proportion of the total diversity (pie charts) and total abundance (bar charts) represented by Linyphiidae species increased from north boreal sites to arctic sites (Fig. 2).

Ecoclimatic regions showed clear differences in species composition (Fig. 3). Variation in composition between sites was higher in the arctic than in the subarctic or north boreal ecoclimatic regions (Fig. 3).

Diversity partitioning results showed that spider diversity was structured (i.e., non-random) at the regional, site, and sample scales with lower diversity than expected by chance at lower scales and higher observed diversity at larger spatial scales (Table 2). The observed values of β4 and β3 were significantly higher than expected with a random distribution (Table 2). The ecoclimatic region was the most important spatial scale in structuring diversity; β4 contributed to 50% of the total diversity of γ diversity meaning that half of the species are restricted to one ecoclimatic region. The site scale was the second most important spatial scale with β3 contributing to 30% of the total richness. Additive partitioning of diversity within each ecoclimatic region did not show differences between regions. β3 contributed to almost half of the diversity within each region (Table 2).

The three explanatory matrices (vegetation, climate, and spatial variables) explained 59.6% of the total variation in spider composition at the regional scale. The vegetation matrix based on local measures explained 38.1% of variation in ground-dwelling spider diversity compared with 35.2% for the climate matrix. Latitude was not an important driver and the two spatial variables explained only 17.1% of the variation in diversity at a large spatial scale. Variance partitioning within each region showed similar results in different proportions. The effects of vegetation explained 61% and 46% of the variance in the arctic and north boreal ecoclimatic regions, respectively.

The constrained ordination was performed using all variables, as no variable was discarded by the permutation selection (Fig. 4a). The first canonical axis explained 16% of the total variance of spider composition, and the first two canonical axes combined explained 26.7% of the variance. The significant variables selected to find parsimonious models were different within each ecoclimatic region (Figs. 4b, 4c, and 4d). Variables selected in each region were a mix of several local variables and one or two climatic variables, with the exception of the arctic ecoclimatic region where longitude was also significant. For the arctic and the north boreal ecoclimatic regions, maximum temperature of the warmest month seemed to influence differences between the eastern and western sites (Figs. 4b and 4d). The first two canonical axes explained 62% of the variance in the arctic (Fig. 4b), 52% in the subarctic (Fig. 4c) and 50% in the north boreal ecoclimatic regions (Fig. 4d).