Journal of Animal Ecology 2014

doi: 10.1111/1365-2656.12287

Elements of regional beetle faunas: faunal variation and compositional breakpoints along climate, land cover and geographical gradients Jani Heino1* and Janne Alahuhta2 1 2

Finnish Environment Institute, Natural Environment Centre, Biodiversity, PO Box 413, FI-90014 Oulu, Finland; and Department of Geography, University of Oulu, PO Box 3000, FI-90014 Oulu, Finland

Summary 1. Regional faunas are structured by historical, spatial and environmental factors. We studied large-scale variation in four ecologically different beetle groups (Coleoptera: Dytiscidae, Carabidae, Hydrophiloidea, Cerambycidae) along climate, land cover and geographical gradients, examined faunal breakpoints in relation to environmental variables, and investigated the best fit pattern of assemblage variation (i.e. randomness, checkerboards, nestedness, evenly spaced, Gleasonian, Clementsian). We applied statistical methods typically used in the analysis of local ecological communities to provide novel insights into faunal compositional patterns at large spatial grain and geographical extent. 2. We found that spatially structured variation in climate and land cover accounted for most variation in each beetle group in partial redundancy analyses, whereas the individual effect of each explanatory variable group was generally much less important in accounting for variation in provincial species composition. 3. We also found that climate variables were most strongly associated with faunal breakpoints, with temperature-related variables alone accounting for about 20% of variation at the first node of multivariate regression tree for each beetle group. The existence of faunal breakpoints was also shown by the ‘elements of faunal structure’ analyses, which suggested Clementsian gradients across the provinces, that is, that there were two or more clear groups of species responding similarly to the underlying ecological gradients. 4. The four beetle groups showed highly similar biogeographical patterns across our study area. The fact that temperature was related to faunal breakpoints in the species composition of each beetle group suggests that climate sets a strong filter to the distributions of species at this combination of spatial grain and spatial extent. This finding held true despite the ecological differences among the four beetle groups, ranging from fully aquatic to fully terrestrial and from herbivorous to predaceous species. 5. The existence of Clementsian gradients may be a common phenomenon at large scales, and it is likely to be caused by crossing multiple species pools determined by climatic and historical factors on the distributions of species. Key-words: assemblage composition, compositional gradients, environmental gradients, idealized models, multivariate regression trees, partial redundancy analysis, regionalisation

Introduction Regional biotas develop under the influences of evolutionary, historical and environmental factors. These factors can be understood as a spatiotemporal continuum, with

*Correspondence author. E-mail: [email protected]

evolutionary factors largely driving differences in species composition among broad geographic regions over long time periods, historical factors such as the most recent ice age causing variation among smaller regions, and environmental factors being important not only at large spatial scales but also within smaller regions (Brown & Lomolino 1998; Mittelbach 2012). Examining the relative importance of these factors for faunal composition in different

© 2014 The Authors. Journal of Animal Ecology © 2014 British Ecological Society

2 J. Heino & J. Alahuhta regions is a prerequisite for understanding the interactions between regional species pools and local communities (e.g. Cornell 2013), delineating regions for environmental assessment and conservation (e.g. Bailey 2010), and predicting the effects of global change on species distributions at large spatial scales (e.g. Elith & Leathwick 2009). New insights into those biogeographical patterns can be attained by applying analytical methods that have been more typically applied in the analysis of local ecological communities. Regional faunas show multiple patterns of interest to ecologists. These include variations in species richness, endemism, rarity and composition, all of which may show various patterns in space and time (Brown & Lomolino 1998; Gaston 2000). Elements of faunal composition can also be evaluated with regard to which idealized assemblage pattern they fit best. Leibold & Mikkelson (2002) devised an approach to examine such assemblage patterns, whereby random distributions of species are contrasted with five idealized models. These idealized models include checkerboards (e.g. Gotelli & McCabe 2002), nestedness (e.g. Patterson & Atmar 1986), evenly spaced gradients (e.g. Tilman 1982), Gleasonian gradients (e.g. Gleason 1926) and Clementsian gradients (e.g. Clements 1916). Although this approach of ‘elements of metacommunity structure’ (EMS) was originally aimed at separating different patterns across a set of local communities (Leibold & Mikkelson 2002), the same approach can be adapted to examine distribution patterns and faunal structure at large biogeographical scales, across regions or across sets of islands (Presley & Willig 2010; Meynard et al. 2013). The five idealized spatial patterns and randomness should be tested simultaneously with the same data set (Leibold & Mikkelson 2002; Dallas & Presley 2014), although most previous studies have tested different patterns in isolation. For example, checkerboards and nestedness are examples of such patterns typically examined in isolation. Checkerboards refer to a situation where species pairs have mutually exclusive distributions, whereby two species replace each other across a set of sites (Diamond 1975; Gotelli & Graves 1996). While most studies have examined checkerboards across a set of islands or patchy mainland habitats (Gotelli & McCabe 2002; G€ otzenberger et al. 2012), there is no reason not to examine checkerboards across large spatial scales on continents. For example, Gotelli, Graves & Rahbek (2010) examined checkerboard distributions of birds based on large grain sizes at a regional extent and found that some species pairs occurred together less often than expected by chance expectations. Hence, checkerboards may also occur at large grain sizes and large geographical extents, although the mechanisms underlying such patterns should differ between large-scale biogeographical studies and smallerscale studies encompassing local communities. The nested subset pattern has also been studied extensively in various systems, ranging from islands through habitat patches to continental habitats (Wright et al.

1998). A perfectly nested pattern occurs when species poor faunas comprise subsets of those in progressively richer faunas, and when the species with the broadest range include those of species with progressively smaller ranges (Patterson & Atmar 1986). Nestedness has been examined at large biogeographic scales, too, and these studies have often found a significantly nested pattern (Baselga 2008). However, very few studies have tried to compare such biogeographical variation with other idealised patterns, which have been studied more extensively across local communities along habitat gradients. The remaining spatial patterns can be divided between three main idealized models. Evenly spaced gradients describe a situation, where species range boundaries are hyperdispersed along the underlying environmental gradient (Tilman 1982), indicating maximal differences among species in environmental tolerances (Presley & Willig 2010). Based on the idea of Clementsian gradients, communities are tight collections of two or more species groups, responding similarly to environmental gradients and hence showing discrete community boundaries (Clements 1916). In contrast, based on the idea of Gleasonian gradients, each species responds to the environmental gradients individualistically and communities are mere collections of species, whose ranges happen to overlap in space and time (Gleason 1926). While the debate over community gradients has been ongoing for almost a century, most attention has been devoted to terrestrial plants and examinations of local communities (Allen & Hoekstra 1992; McIntosh 1995; Hoagland & Collins 1997). It is thus important to test the fit of these idealised models in various other organismal groups, ecological systems and spatial scales (Henriques-Silva, Lindo & Peres-Neto 2013; Meynard et al. 2013). The EMS approach is mainly a descriptive, pattern-based approach and it does not tell us much about the underlying mechanisms responsible for the patterns detected. There has hence been recent interest in accounting for variation in community composition using spatial and environmental predictors (Cottenie 2005; Meynard et al. 2013). While most studies have taken a continuum perspective, whereby community composition has been examined using direct gradient analysis methods, such as redundancy analysis (RDA; Dray et al. 2012; Meynard et al. 2013), recently developed constrained clustering methods, such as multivariate regression trees, allow not only description and prediction of variation in community composition, but also provide a means to find breakpoints of variation in assemblage composition (De’ath 2002; Borcard, Gillet & Legendre 2011). Such breakpoints may occur when there is a clear change in assemblage composition that is typically associated with a clear change in environmental conditions. Surprisingly, very few studies have assessed both continuous variation and breakpoints of community composition in the same study, and these studies have mainly considered local habitat gradients within a small region (Davidson et al. 2010; Sahuquillo & Miracle 2013; Heino, Ilmonen &

© 2014 The Authors. Journal of Animal Ecology © 2014 British Ecological Society, Journal of Animal Ecology

Elements of regional beetle faunas 3 Paasivirta 2014). However, it is at large spatial grain sizes and geographical extents where the occurrence of breakpoints in faunal compositions should be most easily discernible, because studies spanning large spatial scales are likely to cross multiple regional species pools determined by the effects of historical factors and climatic constraints on the distributions of species (Brown & Lomolino 1998; Heino 2011). Hence, those breakpoints should also be examined using large-scale biotic, climatic and land cover data. Such examinations are also valuable owing to the fact that they provide important information about biogeographic regionalisation for environmental assessment and biodiversity conservation (Bailey 2010; Rueda, Rodriguez & Hawkins 2010). We examined geographical variation and breakpoints in the regional faunal composition of four beetle groups (Coleoptera: Carabidae, Dytiscidae, Hydrophiloidea, Cerambycidae), showing a wide variety of habitats (both aquatic and terrestrial groups) and trophic modes (ranging from predators through omnivores to herbivores). Our study region comprised the whole Fennoscandia and Denmark, where large gradients in climatic conditions (i.e. from temperate to arctic regions), land cover (i.e. from deciduous forest to tundra landscapes) and land use (i.e. heavily anthropogenically modified to near-pristine regions) occur in a relatively small geographical area. We specifically examined the following questions: (i) What are the contributions of climatic, land cover and spatial variables to variation in faunal composition? (ii) What environmental factors best determine the breakpoints in faunal composition if such breakpoints exist? (iii) Which idealized model best fits the empirical data of the four beetle families? (iv) Are there differences in these patterns among aquatic (Dytiscidae, Hydrophiloidea) and terrestrial (Carabidae, Cerambycidae) beetles? Based on evidence from previous studies, we expected that the aquatic and terrestrial groups would respond relatively similarly to climatic gradients (V€ais€anen, Heli€ ovaara & Immonen 1992; Heino 2001), but show different responses to land cover gradients due to different habitat requirements.

biogeographical analyses (V€ais€anen, Heli€ ovaara & Immonen 1992). Thus, although most insect biodiversity studies are affected by variation in detection of species, we are confident that these data adequately characterize provincial assemblages and are representative of true ecological patterns. Furthermore, there are actually certain benefits of using biogeographical province as sampling units rather than using equal area map grids. First, any change in the spatial grain of the grid system has strong impacts on the final outcomes. Secondly, in case of fragmented data, most of such grid cells tend to appear empty, although they are very likely to be occupied by the target species. Use of more natural geographical units may hence be superior to grid systems, especially when distributional data are fragmented, as is the case of most insect groups. We hence used the 79 provinces as sampling units (i.e. grain size) and, hence, the species recorded from each province were pooled to represent a single assemblage.

the study organisms and species data We analysed the literature data for four beetle groups, including ground beetles (Carabidae; Lindroth 1985, 1986), diving beetles (Dytiscidae; Nilsson & Holmen 1995), longhorn beetles (Cerambycidae; Bily & Mehl 1989) and water scavenger beetles (Hydrophiloidea; Holmen 1987). Ground beetles comprise a large family of mainly terrestrial beetles, which are predatory, omnivorous, granivorous or herbivorous species as adults and mostly predaceous as larvae (Lindroth 1985). Diving beetles occur in fresh waters and, in some cases, brackish water environments, and they are mostly predaceous as larvae and predators or scavengers as adults (Nilsson & Holmen 1995). Longhorn beetles are primarily forest insects that live in timber as larvae and feed on various types of plant material as adults (Bily & Mehl 1989). Water scavenger beetles are mainly predaceous, scavenging or herbivorous as adults and predaceous as larvae. They occur chiefly in freshwater or semi-aquatic environments (Holmen 1987). These four groups thus comprise main life styles and habitats amongst beetles, allowing us to analyse variation in ecologically highly different, yet phylogenetically relatively closely related organisms. The diversity of these four beetle groups is considerable across the 79 provinces considered in this study. The most diverse group is Carabidae (total number of species = 388; mean = 159, SD = 56
9), followed by Dytiscidae (total number of species = 155; mean = 78
9, SD = 19
3), Cerambycidae (total number of species = 118; mean = 46
5, SD = 17
6) and Hydrophiloidea (total number of species = 116; mean = 52
4, SD = 22
6).

Materials and methods study area We analysed biological and environmental data from the 101 biogeographic provinces belonging to Denmark, Sweden, Norway and Finland (V€ ais€anen, Heli€ ovaara & Immonen 1992; V€ais€anen & Heli€ ovaara 1994). Prior to the analyses, we merged small coastal provinces in Norway to provide a better and more accurate representation of species ranges. We also omitted a single small Swedish island province (Gotska Sand€ o), as it was a clear outlier in preliminary ordination analyses. After these modifications, the number of provinces to be analysed decreased from 101 to 79. Each province has typical characteristics of climate and land cover, thereby ‘province’ being a relatively homogeneous study unit. The insect faunas in the study area are relatively well known, offering good data for insect groups used infrequently in

climate variables Climate variables included average annual temperature (°C), maximum temperature of the warmest month (°C), minimum temperature of the coldest month (°C), precipitation of the wettest month (mm) and precipitation of the driest month (mm). The climate variables were average values for each biogeographical province and were derived from WORLDCLIM with 0
93 km 9 0
93 km resolution (Hijmans et al. 2005).

land cover variables Land cover and land use variables consisted of percentages of fresh water, forests, open areas, wetlands, agricultural areas and urban development. Land cover and land use variables were

© 2014 The Authors. Journal of Animal Ecology © 2014 British Ecological Society, Journal of Animal Ecology

4 J. Heino & J. Alahuhta obtained from European CORINE with 100 m resolution. Although European CORINE data are temporally overlapping, yet of recent origin, in comparison to the insect data that have accumulated within a longer time period, we considered our data to represent well provincial differences in land use. This is because the southernmost provinces (i.e. which are most heavily burdened by human land use nowadays) have differed from the northernmost provinces (i.e. which are very little modified even nowadays) in land use for a long time. For example, major deforestation had already began in Denmark 3000 years ago (Bradshaw 2005), whereas widespread human presence was not evidenced in the nature of northern Scandinavia until 500 years ago (van der Linden et al. 2008). We also used distance to the sea as a variable representing marine influence on the insect fauna. The starting point for the distance measurement was the midpoint of a province, which is located at the centre of gravity (i.e. centroid) of the polygon (i.e. true centroid of the province). Finally, average elevation and elevation range within the province were also included as land cover variables, as these variables are related to the range of habitat types along elevation gradients. Elevation variables were obtained from 3D Digital Elevation Model over Europe with 25 m resolution.

province area We considered province area (km2) as a separate group of predictor variables, because at this scale, it is not clearly mechanistically associated with the environmental variables we considered. Province area may also be considered as a proxy for sampling effects (i.e. a larger province is likely to sample more species than a smaller province), although sampling effects should be less important for overall faunal structure (i.e. our focus in this study) than for species richness.

spatial variables We used an eigenfunction spatial analysis by means of the family of Moran’s eigenvector maps to model spatial structures among the provinces (Dray, Legendre & Peres-Neto 2006; Griffith & Peres-Neto 2006). To do this, we ran principal coordinates of neighbour matrix analyses (PCNM) based on Euclidean distance among the provinces as input (Borcard & Legendre 2002; Borcard, Gillet & Legendre 2011). We retained the PCNM eigenvectors with positive eigenvalues as explanatory variables in further analyses (see below). Basically, the first PCNMs with large eigenvalues describe broadscale spatial structures, whereas those PCNMs with small eigenvalues describe fine-scale spatial variation (Peres-Neto & Legendre 2010). The PCNMs are mutually orthogonal and linearly unrelated spatial variables and can be used to model spatial variation in provincial species composition in this study. Significant spatial variation as related to PCNMs may be a consequence of environmental autocorrelation, dispersal limitation or historical effects on provincial species composition (Dray et al. 2012; Peres-Neto, Leibold & Dray 2012). PCNM analysis was conducted using the R package PCNM (Legendre et al. 2013).

statistical methods We first examined variation in species composition as opposed to turnover in species composition (for definitions, see Anderson

et al. 2011) across the provinces, and therefore, we chose to use constrained ordination and constrained clustering methods as our main statistical tools (Legendre & Legendre 2012). Our main statistical method was canonical RDA that analyses variation in species composition in relation to explanatory variables (Rao 1964). Prior to the RDAs, species data of each beetle group were Hellinger-transformed to account for numerous zero values and making the data analysable using linear methods (Legendre & Gallagher 2001). We first selected variables in the final RDA models of each set of variables (i.e. land cover, climate, spatial) based on the forward selection method proposed by Blanchet, Legendre & Borcard (2008). We next used partial redundancy analysis (pRDA) to decompose variation in the provincial species composition of each beetle group among climate, land cover, spatial variables and province area following the protocols explained in detail elsewhere (Borcard, Legendre & Drapeau 1992; Anderson & Gribble 1998; Legendre & Legendre 2012). Variation partitioning among four sets of predictor variables results in individual climate, individual land cover, individual spatial and province area fractions, as well as their shared effects and unexplained variance. Adjusted coefficient of determination (adj. R2) was reported in all variation partitioning analyses, because they are unbiased measures of explained variation (Peres-Neto et al. 2006). RDAs and variation partitioning were conducted using the R package vegan (Oksanen et al. 2013). We constructed Mantel correlograms (Sokal & Oden 1978) for each beetle group to examine spatial autocorrelation in faunal composition across the provinces. Mantel correlograms are a multivariate version of univariate Moran I correlograms and have been shown to be effective in portraying spatial patterns (Borcard & Legendre 2012). We used Hellinger-transformed species data as the basis of Mantel correlograms. In these analyses, the number of distance classes was based on Sturge’s rule (Borcard, Gillet & Legendre 2011), and the significance of Mantel correlations at each distance class was based on Holm correction (Holm 1979). Mantel correlograms were done using the R package vegan (Oksanen et al. 2013). We examined breakpoints of faunal variation along climatic and land cover gradients using multivariate regression tree analysis (MRT; De’ath 2002). We used MRT based on Euclidean distance of Hellinger-transformed multivariate species data (response variables) and province area, climate and land cover variables as explanatory variables. We also ran trial analyses with latitude and longitude included among the explanatory variables, but their inclusion did not lead to clearly better MRT models. MRT forms clusters of sites by repeatedly splitting the data, each split is defined by a simple rule based on the values of explanatory variables, and the splits are chosen to minimize the dissimilarity of sites within clusters (De’ath 2002; Borcard, Gillet & Legendre 2011). MRT results in a tree whose terminal site groups or ‘leaves’ are composed of subsets of sites that are chosen to minimize the within-group sums of squares. Each successive partitioning of data is defined by a threshold value of one of the explanatory variables. This method retains a solution with the greatest predictive power and can handle a wide variety of situations, including a situation where assemblage–environment relationships are nonlinear (Borcard, Gillet & Legendre 2011; Legendre & Legendre 2012). In brief, the computation of MRTs consists of (i) constrained partitioning of the data and (ii) cross-validation of the results (De’ath 2002). In this study, we chose the ‘best’ tree with the minimum cross-validated error (CVRE). MRTs were constructed using the R package MVPARTwrap (Ouellette & Legendre 2013).

© 2014 The Authors. Journal of Animal Ecology © 2014 British Ecological Society, Journal of Animal Ecology

Elements of regional beetle faunas 5 We also examined elements of provincial faunal structure based on the original EMS approach and the ‘range perspective’ in our analyses (Leibold & Mikkelson 2002). The analysis of the elements of faunal structure was based on three metrics: (i) coherence, (ii) turnover and (iii) boundary clumping. Prior to calculating these three metrics, provinces-by-species presence– absence matrix is ordinated through reciprocal averaging (i.e. correspondence analysis; Gauch 1982). Correspondence analysis defines a latent gradient, ordering sites and species along that gradient (Leibold & Mikkelson 2002; Presley & Willig 2010). The first metric that we evaluated is coherence. It is based on calculating the number of embedded absences in the ordinated matrix and, subsequently, comparing the observed value to a null distribution of embedded absences from randomizations. An embedded absence means that there is a gap in a species range (Leibold & Mikkelson 2002; Presley & Willig 2010). A small number of embedded absences lead to positive coherence, while a large number of embedded absences mean negative coherence. Significantly negative coherence thus provides support to a checkerboard distribution of species, non-significant coherence to randomness, and significantly positive coherence to nestedness, evenly spaced gradients, Gleasonian gradients or Clementsian gradients (Leibold & Mikkelson 2002). A second metric is called turnover. This metric is evaluated if coherence is positive, and it is used for deciding which gradient model best fits the data. Species turnover is measured as the number of times one species replaces another between two sites in an ordinated matrix (Presley & Willig 2010). Significantly negative turnover is related to nestedness, whereas non-significant (i.e. moderately positive) or significantly positive turnover points to evenly spaced, Gleasonian or Clementsian gradients (Leibold & Mikkelson 2002; Henriques-Silva, Lindo & Peres-Neto 2013). These last three types of gradients can be separated based on evaluation of boundary clumping (Leibold & Mikkelson 2002; Presley, Higgins & Willig 2010). Boundary clumping is analysed using Morisita’s index and a chi-square test comparing observed and expected distributions of range boundary locations. Values of this index that are not different from 1 indicate randomly distributed range boundaries (i.e. Gleasonian gradients), values significantly larger than 1 indicate clumped range boundaries (i.e. Clementsian gradients) and values significantly < 1 indicate hyperdispersed range boundaries (i.e. evenly spaced gradients). The significance of the index values for coherence and turnover was tested using the fixed-equiprobable null model (Gotelli 2000), where the species richness of each site was maintained (i.e. row sums are fixed), but species occurrences were assumed to be equiprobable (i.e. column sums are not fixed). Because species richness typically varies along environmental gradients, this null model makes sense ecologically (Presley et al. 2009). We also used a strict and conservative fixed–fixed null model to find out if the null model used affected the results. In the fixed–fixed null model, the species richness of each site was maintained (i.e. row sums are fixed) and species occurrences were the same as species frequencies of occupancy (i.e. column sums are fixed). We used 999 simulations to provide random matrices for testing coherence and turnover. Elements of faunal structure were evaluated for each beetle group along primary (axis 1) correspondence analysis axis based on each null model. Elements of faunal structure were analysed using the R package metacom (Dallas 2013) in the R environment (version 3.0.1, R Core Team 2013). Although the original EMS approach (Leibold & Mikkelson 2002) has been criticized by some authors on the basis that a

data set may simultaneously show different significant patterns (Ulrich & Gotelli 2013), other researchers have recently shown that it is a good means to indicate ‘the best fit’ of an empirical data set with one of the idealized patterns (Meynard et al. 2013; Dallas & Presley 2014; de la Sancha et al. 2014). Furthermore, one major difference between the EMS approach and the traditional nestedness and co-occurrence analyses is that the former approach examines different patterns along a single major gradient in the data, while the latter two analyses examine a given pattern in the whole sites-by-species matrix simultaneously (Presley, Higgins & Willig 2010; Presley & Willig 2010). We used a combination of k-means partitioning, chi-square tests and correspondence analysis to assess the robustness of the MRT results and to associate the findings from the MRT and EMS analyses regarding faunal breakpoints. We hence used k-means partitioning to cluster the provinces based on the number of groups corresponding those from the MRTs at the first node and based on the number of final MRT leaves. K-means partitioning works by forming groups through identifying high-density regions in the data and, contrary to other clustering methods, it is based on a pre-determined number of final groups (Borcard, Gillet & Legendre 2011). We used either Hellinger-transformed or chi-square-transformed data as the input matrix in k-means partitioning (for the transformations, see Legendre & Gallagher 2001). The Hellinger transformation connected the Hellingerbased RDAs and MRTs and the chi-square transformation connected the findings with correspondence analysis that is based on chi-square distance (Gauch 1982). We compared the matches between the MRT and k-means clusters using chi-square tests (Borcard, Gillet & Legendre 2011). We also compared the matches between the four beetle groups using chi-square tests. K-means clustering was done using the R package Rcmdr (Fox 2005) and chi-square tests with permutation using the R package coin (Hothorn et al. 2008). Finally, we plotted the province and species scores of each beetle group on correspondence analysis ordination plots to show that the MRT breakpoints at the first node were associated with breaks in species distributions and provincial species composition detected in the EMS analysis.

Results Eigenfunction spatial analysis produced 28 eigenvectors with positive eigenvalues. These ranged from those describing large-scale spatial trends (i.e. those with large eigenvalues) to those related to small-scale spatial relationships among sites (i.e. those with small eigenvalues). These spatial variables were subsequently used to describe variation in the provincial species composition of each beetle group. The spatial models in the RDAs included the following numbers of significant spatial variables: diving beetles 16 spatial variables, ground beetles 20 spatial variables, water scavenger beetles 17 spatial variables and longhorn beetles 13 spatial variables (Appendix S1, Supporting information). In general, the spatial eigenvector variables (SEVs) with the largest eigenvalues (SEV1 and SEV2) were selected first for each beetle group, followed by a variable combination of intermediate and small-scale spatial variables depending on the beetle group. Highly similar sets of significant climate and land cover variables were selected in the final RDA models of each

© 2014 The Authors. Journal of Animal Ecology © 2014 British Ecological Society, Journal of Animal Ecology

6 J. Heino & J. Alahuhta beetle group (Appendix S1, Supporting information). Variation in the provincial species composition of diving beetles was significantly accounted for by (in order of importance) mean annual temperature, maximum temperature, precipitation of the wettest month and minimum temperature, as well as by agriculture, forests, wetlands, altitude range and open area. For ground beetles, the selected climate variables were the same as for diving beetles, but the set of land cover variables included agriculture, forests, wetland, altitude range, open area and water systems area. For longhorn beetles, minimum temperature, maximum temperature, mean annual temperature and precipitation of the wettest month were significant climate variables, whereas agriculture, forests, wetland, altitude range, open area, urban area and water systems area were significant land cover variables. The provincial species composition of water scavenger beetles was closely associated with mean annual temperature, maximum temperature, precipitation of the wettest month and minimum

temperature, as well as significantly related to agriculture, forests, wetlands, altitude range, open area and water systems area. Variation partitioning based on pRDAs showed that the three main groups of predictor variables and province area accounted for variation in the species composition of each beetle group highly similarly (Table 1). Thus, the contributions of climate, land cover and spatial variables were very high, accounting for 31–51% of variation in provincial species composition. Of the individual effects of the explanatory variable groups (Table 1), spatial variables were clearly the most important ones (variation accounted for: 7
5–11
6%), followed by the minor individual effects of climate (0
7–2
7%) and land cover variables (0
9–1
8%). There was also joint variation among two of the variable groups, and even these joint effects typically slightly overcame the individual effects of climate and land cover variables (Table 1). Province area had clearly weaker effects on provincial species composition than any

Table 1. Results of variation partitioning of species composition in each beetle group among four predictor variable groups (land cover, climate, spatial and province area). Shown are total, individual and joint fractions of variation. Province area was not considered an environmental variable, but was included to account for the species–area relationship. Significant individual fractions are in bold Diving beetles

Fractions Land cover = L Climate = C Spatial = S Area = A L+C L+S L+A C+S C+A S+A L+C+S L+C+A L+S+A C+S+A L+C+S+A Individual fractions L| C+S+A C| L+S+A S| L+C+A A| L+C+S Joint fractions L∩C C∩S L∩S L∩A C∩A S∩A L∩C∩A L∩C∩S C∩S∩A L∩S∩A C∩L∩S∩A Residuals

Ground beetles

Scavenger beetles

Longhorn beetles

d.f.

Adj. R2

d.f.

Adj. R2

d.f.

Adj. R2

d.f.

Adj. R2

5 4 16 1 9 21 6 20 5 17 25 10 22 21 26

0
376 0
424 0
462 0
188 0
459 0
509 0
403 0
526 0
426 0
465 0
534 0
462 0
518 0
528 0
537

6 4 20 1 10 26 7 24 5 21 30 11 27 25 31

0
364 0
398 0
513 0
192 0
443 0
543 0
388 0
546 0
403 0
52 0
559 0
449 0
549 0
554 0
564

6 4 17 1 10 23 7 21 5 18 27 11 24 22 28

0
332 0
343 0
419 0
139 0
400 0
466 0
351 0
456 0
347 0
425 0
475 0
404 0
475 0
466 0
482

7 4 13 1 11 20 8 17 5 14 24 12 21 18 25

0
314 0
342 0
398 0
157 0
397 0
441 0
344 0
450 0
348 0
407 0
473 0
402 0
451 0
459 0
477

5 4 16 1

0
009 0
019 0
075 0
003

6 4 20 1

0
010 0
016 0
116 0
006

6 4 17 1

0
015 0
007 0
077 0
007

7 4 13 1

0
018 0
027 0
076 0
005

0 0 0 0 0 0 0 0 0 0 0 0

0
045 0
040 0
027 0
000 0
006 0
000 0
006 0
134 0
017 0
001 0
169 0
463

0 0 0 0 0 0 0 0 0 0 0 0

0
018 0
044 0
035 0
002 0
000 0
000 0
000 0
133 0
019 0
000 0
168 0
436

0 0 0 0 0 0 0 0 0 0 0 0

0
034 0
047 0
042 0
004 0
002 0
002 0
007 0
120 0
013 0
003 0
127 0
518

0 0 0 0 0 0 0 0 0 0 0 0

0
025 0
031 0
035 0
004 0
005 0
000 0
005 0
109 0
020 0
002 0
131 0
522

© 2014 The Authors. Journal of Animal Ecology © 2014 British Ecological Society, Journal of Animal Ecology

Elements of regional beetle faunas 7 0

Coefficient of determination (R 2)

Mean annual temp (oC) < 3·035

Mean annual temp (oC) >3·035

(1) 24·33%

24·33 Open area (%) >0·185

Open area (%) < 0·185

(2) 6·74%

31·07 Prec. of driest month (mm) >54·09

Prec. of driest month (mm) < 54·09

(3) 5·37%

36·44 Min temp (oC)< –5·755

(7) 4·99%

Min temp (oC) >–5·755

41·43 Agriculture (%) < 0·015 45·78

Alt range (m) >1783

48·12

Agriculture (%) >0·015

(5) 4·35%

(11) 2·35%

Alt range (m) < 1783

n = 11

n = 11

n=3

n = 13

n=3

n = 19

n = 19

(#3)

(#5)

(#7)

(#8)

(#10)

(#12)

(#13)

R 2 × 100 = 48·1 % Error = 0·519 CV Error = 0·702 SE = 0·0371 Fig. 1. Multivariate regression tree for diving beetles. This is the ‘best’ tree with minimum cross-validated error.

of the three main predictor variable groups (Table 1). The total amounts of variation attributable to the four explanatory variable groups were relatively high, ranging from 47
7% for longhorn beetles to 56
4% for ground beetles. The Mantel correlograms were highly similar among the beetle groups (Appendix S2, Supporting information). There was significant positive autocorrelation at the smallest distance classes, no significant autocorrelation at the middle classes and significant negative autocorrelation at the largest distance classes. The MRT analyses showed that, for each beetle group, temperature variables were most clearly associated with faunal breakpoints. For diving beetles, mean annual temperature alone accounted for 24% of variation in species composition, followed by percentage open areas,

precipitation of the driest month, minimum temperature, percentage agriculture cover and altitude range (Fig. 1). These other variables were much weaker contributors to the total explained variation of 48% than mean annual temperature. For ground beetles, mean annual temperature was again associated with the strongest faunal breakpoint, accounting for 23% of variation in species composition, followed by clearly minor contributions (< 8%) to the total explained variance of 51% by percentage urban areas, mean annual temperature, percentage agriculture cover, distance to the sea, altitude range and minimum temperature (Fig. 2). For water scavenger beetles, maximum temperature accounted for 19% of variation, followed by minimum temperature, maximum temperature, distance to the sea and percentage agricultural

© 2014 The Authors. Journal of Animal Ecology © 2014 British Ecological Society, Journal of Animal Ecology

8 J. Heino & J. Alahuhta

Coefficient of determination (R 2)

0

23·61

Mean annual temp (oC)< 3·395

Urban (%) < 0·005

Urban (%) >0·005

(2) 7·79%

31·4

37·71

41·25 44·39

Mean annual temp (oC) >3·395

(1) 23·61%

Mean annual temp (oC) < 7·015

Agriculture (%) < 0·035 Sea (km) < 334·5

(8)

Mean annual temp (oC) >7·015

Agriculture (%) >0·035

(4) 3·54%

3·14%

(3) 6·31%

Sea (km) >334·5 Altrange (m) < 938·2

(5)

Altrange (m) >938·2

2·90%

47·29

Altrange (m) >1099

49·71 51·61

Altrange (m) < 1099

(6) 2·42% Min temp (C)< –8·668

(13) 1·90%

Mintemp (C) >–8·668

n=8

n = 10

n=3

n = 12

n=7

n=3

n=8

n = 12

n = 16

(#5)

(#6)

(#7)

(#9)

(#10)

(#13)

(#15)

(#16)

(#17)

R 2 × 100 = 51·6 % Error = 0·484 CV Error = 0·716 SE = 0·0338 Fig. 2. Multivariate regression tree for ground beetles. This is the ‘best’ tree with minimum cross-validated error.

cover (Fig. 3). These variables explained 43% of variation in species composition. For longhorn beetles, minimum temperature alone accounted for 21% of variation in species composition, followed by percentage forest cover, percentage urban land cover, maximum temperature and distance to the sea (Fig. 4). These variables accounted for a total of 41% of variation in species composition. The results of the MRT analyses remained largely the same when latitude and longitude were included in the set of the explanatory variables. Their inclusion did not significantly increase the explanatory power of the MRTs, and they appeared only in the lower nodes of the trees (results not shown). The EMS analysis indicated that each beetle group showed (i) positive coherence (i.e. the number of embedded absences was lower than expected by chance), (ii) more turnover than expected by chance (i.e. the number

of replacements was higher than expected by chance) and (iii) significantly higher boundary clumping than 1 (i.e. based on Morisita’s I) (Table 2). The same result held true for both fixed-equiprobable and fixed–fixed null models. These findings indicated that the beetle faunas fitted best with Clementsian gradients along the primary axis of variation in provincial species composition. The groupings based on the MRTs and k-means partitioning were very similar, and significant chi-square tests strengthened this conclusion (Table 3). This result held with both distance transformations of species data that were used as the basis of k-means partitioning. Chi-square tests also indicated that the breakpoints at the first MRT node were highly similar between the beetle groups (Table 4; see also Appendices S3 and S4, Supporting information), and the same was true for k-means clusters (Table 4). Finally, plotting the two groups of the first

© 2014 The Authors. Journal of Animal Ecology © 2014 British Ecological Society, Journal of Animal Ecology

Elements of regional beetle faunas 9

Coefficient of determination (R 2)

0

Maxtemp (oC) < 18·8

19·42

27·86

Min temp (oC) < –9·762

Min temp (oC) < –13·67

35·11 37·94

Max temp (oC) >18·8

(1) 19·42%

40·75

8·44%

Min temp (oC) >–13·67

(2) 7·25%

Max temp (oC) < 14·57 Sea (km) < 360

Min temp (oC) >–9·762

(3)

Max temp (oC) >14·57

(5) 2·83%

Sea (km) >360

(4) 2·81%

Agriculture (%) < 0·43

43·50

(7) 2·75%

Agriculture (%) >0·43

n=7

n=6

n=5

n = 11

n = 20

n = 18

n = 12

(#4)

(#5)

(#7)

(#8)

(#10)

(#12)

(#13)

R 2 × 100 = 43·5 % Error : 0·565 CV Error : 0·78 SE : 0·0465 Fig. 3. Multivariate regression tree for water scavenger beetles. This is the ‘best’ tree with minimum cross-validated error.

MRT node on CA ordination plots clearly indicated that the breaks were strongly associated between those two methods (Fig. 5), also indicating that MRT and EMS found the same breaks in the beetle data sets.

Discussion While a large number of studies have examined variation in species richness across large spatial grains and extents (see review by Field et al. 2009), fewer studies have examined species compositional variation at such broad scales (V€ais€anen, Heli€ ovaara & Immonen 1992; Baselga 2008; Rueda, Rodriguez & Hawkins 2010). We found that approximately half of variation in the provincial species composition of the four beetle groups was accounted for by the four groups of explanatory variables. In particular,

the shared effects of the variable groups were considerable in partial RDAs, implying that both climate and land cover vary geographically across the study area and affect provincial species composition simultaneously. In fact, there are clearly parallel gradients in land cover and climate across the study area. Hence, it is difficult to show decisively which explanatory variable set is mostly responsible for variation in species composition. However, given that both climate and history often are related to variation in species richness (Gaston 2000; Fattorini & Ulrich 2012a; Cornell 2013) and species composition (V€ ais€ anen, Heli€ ovaara & Immonen 1992; Heino 2001; Fattorini & Ulrich 2012b) along large geographical gradients, we assume that they are also responsible for variation in beetle species composition across the provinces in our study area. These findings concur with previous studies in the

© 2014 The Authors. Journal of Animal Ecology © 2014 British Ecological Society, Journal of Animal Ecology

10 J. Heino & J. Alahuhta

Coefficient of determination (R 2)

0

Min temp (oC) < –8·59

21·06

27·73

Min temp (oC) >–8·59

(1) 21·06%

Forests (%) < 0·12

Urban (%) < 0·005

Urban (%) >0·005

(2) 6·39%

34·12

38·38

Forests (%) >0·12

(3) 6·66%

Max temp (oC) < 15·64

Sea (km) < 334·5

(7) 4·25%

Max temp (oC) >15·64

Sea (km) >334·5

(4) 3·11%

41·49 n=9

n = 10

n = 26

n = 11

n=5

n = 18

(#4)

(#5)

(#6)

(#8)

(#10)

(#11)

R 2 × 100 = 41·5 % Error : 0·585 CV Error : 0·762 SE : 0·0379 Fig. 4. Multivariate regression tree for longhorn beetles. This is the ‘best’ tree with minimum cross-validated error.

same area, which have suggested the influences of climate on the regionalisation of various animal and plant groups (Lahti, Kurtto & V€ais€anen 1988; V€ais€anen, Heli€ ovaara & Immonen 1992; Heino 2001). Our findings are also similar to those at a broader spatial scale across Europe, where climatic and historical factors were important in determining species richness, endemism and composition of regional faunas of longhorn beetles (Baselga 2008) and darkling beetles (Fattorini & Baselga 2012). In contrast, the degree to which land cover affected species composition is more difficult to infer. We could expect that if land cover was important, we would have seen differences in the particular land cover variables between aquatic and terrestrial beetles. This was not the case, however, and

generally the same land cover variables appeared in the final constrained ordination models. There was significant positive autocorrelation at the smallest distance classes and significant negative autocorrelation at the largest distance classes in all beetle groups. This finding suggests (i) that contiguous provinces have similar species composition, as a consequence of the fact that that neighbouring provinces share similar environmental conditions and (ii) that provinces far from each other have increasingly dissimilar faunas, suggesting that climatic differences between provinces increase faunal differences. Thus, it can be assumed that due to their climatic similarities or differences, significant spatial autocorrelation results from climate forcing on species

© 2014 The Authors. Journal of Animal Ecology © 2014 British Ecological Society, Journal of Animal Ecology

Elements of regional beetle faunas 11 Table 2. Results of the ‘elements of metacommunity structure’ analysis of patterns in provincial species composition. All analyses on the first correspondence analysis axis based on each null model point strongly towards Clementsian gradients Coherence

Null model Fixed-Equiprobable Diving beetles Ground beetles Scavenger beetles Longhorn beetles Fixed–Fixed Diving beetles Ground beetles Scavenger beetles Longhorn beetles

Turnover

P-value

Sim mean

Sim SD

Rep

26
71 31
21 21
85 15
55

0
001 0
001 0
001 0
001

5475 16 406 4405 4833

83 255 99 125

746 4 282 343 420

18
43 22
31 8
89 8
20

0
001 0
001 0
001 0
001

5102 14 113 3262 3989

101 254 116 135

746 4 282 343 420

Abs

Z-value

3234 8438 2224 2875 3234 8438 2224 2875

Boundary clumping

Z-value

P-value

Sim mean

Sim SD

Index

P-value

d.f.

882 211 100 386

69
67 65
50 24
87 19
34

0
001 0
001 0
001 0
001

26 168 30 49

10 62 12 19

346 794 560 159

2
45 2
85 8
28 4
43