**Towards a Better Understanding of Beta Diversity: Deconstructing Composition Patterns of Saproxylic Beetles Breeding in Recently Burnt Boreal Forest**

Ermias T. Azeria1, Jacques Ibarzabal2, Jonathan Boucher2 and Christian Hébert1 *1Natural Resources Canada, Canadian Forest Service, Laurentian Forestry Centre,Québec,*  <sup>2</sup>*Université du Québec à Chicoutimi, Département des sciences fondamentales, Québec, Canada* 

#### **1. Introduction**

74 Research in Biodiversity – Models and Applications

United Nations. (2008). Population Division of the Department of Economic and Social

United Nations Environment Programme-World Conservation Monitoring Centre. (2008).

Affairs of the United Nations Secretariat, world population prospects: the 2008

UNEP-WCMC Species Database. 20.2.2010, Available from: http://www.unep-

United Nations (2006). UN data. 20.2.2010, Available from:

wcmc.org/isdb/Taxonomy/

http://data.un.org/Data.aspx?d=CDB&f=srID%3A29922

revision. 20.2.2010, Available from: http://esa.un.org/unpp

Biodiversity patterns vary across space and time, and this variation is thought to be driven by several ecological/biogeographical processes that act upon species distributions and leave their imprint at various spatial and temporal scales (Huston, 1994; Ricklefs, 2004). The regional or geographical (gamma diversity) diversity patterns, in addition to local diversity (alpha diversity), are consequent upon the extent of "spatial" variation in species composition among sites (beta diversity) (Whittaker, 1960; 1972). Arrays of ecological hypotheses have been proposed as determinants of beta diversity patterns. It can arise from variation in environmental/habitat heterogeneity among sites, spatial constraints, disturbance regimes, and corresponding species-level variations in life history traits (e.g., dispersal, niche-breadth) (Harrison *et al.*, 1992; Nekola & White, 1999; Chase, 2007; Veech & Crist, 2007; Baselga, 2008) as well as neutral and/or stochastic processes (Hubbell, 2001; Chase, 2010). Nevertheless, beta diversity patterns and the mechanisms by which ensembles of local communities maintain their variations and, consequently, influence diversity at regional scale remains a central challenge in ecological and conservation studies (Wilson & Shmida, 1984; Veech & Crist, 2007; Baselga, 2010; Chase, 2010; Tuomisto, 2010; Anderson *et al.*, 2011).

Despite its apparent conceptual simplicity, the empirical assessment of beta diversity has been mired with extensive debates over a perplexing array of indices, analytical approaches, or scales of analysis (Wilson & Shmida, 1984; Lande, 1996; Loreau, 2000; Koleff *et al.*, 2003; Tuomisto & Ruokolainen, 2006; Legendre *et al.*, 2008; Veech & Crist, 2010; for a detailed outline of the issues, see recent review by Anderson *et al.*, 2011). One of the most fundamental and recurrent challenges in beta diversity studies has been distinguishing and quantifying the pure "spatial" turnover component of beta diversity from that caused by variation in species richness between local communities (Wilson & Shmida, 1984; Koleff *et al.*, 2003; Veech & Crist,

Towards a Better Understanding of Beta Diversity: Deconstructing

caused by a harsh environment (also see Chase *et al.*, 2011).

diversity expected under null distribution given richness gradients.

defined primarily by tree-species, burn-severity, and tree-size classes.

an emerging ecological issue in forest management (Lindenmayer *et al.*, 2008).

The study was conducted in 72 sites within four forest burns (burned in 2005) within the western spruce–moss bioclimatic subdomain of northwestern Quebec, Canada (49o15'-50 o

**2. Material and methods** 

**2.1 Study area and bole sampling** 

Composition Patterns of Saproxylic Beetles Breeding in Recently Burnt Boreal Forest 77

similar than expected by random chance, probably due to a deterministic filtering effect

Clearly, a common ground between the partitioning framework and null model analysis of beta diversity needs to be established in order to promote the state of the art and practical understanding of beta diversity patterns. The idea behind null model analysis is to disentangle beta diversity patterns beyond richness gradients, which should be reflected in the turnover component of beta diversity. The nestedness component of beta diversity, on the other hand, is primarily driven by richness differences, which may be related to beta

The present study examines the overall beta diversity components as well as the turnover and nestedness-driven components of beta diversity of saproxylic beetles emerging from tree boles following forest fire. In the present study, we have two major goals. First, we will establish the relationship of the overall beta diversity as well as its two components, turnover and nestedness, with beta diversity expected under and beyond random assembly of communities. We partition or deconstruct the overall beta diversity pattern following the partitioning framework of Baselga (2010). We perform null model analyses for beta diversity with a similar methodological basis as in Chase *et al.* (2011) to estimate the extent of beta diversity expected under and beyond random assembly, given variations in species richness among sites and the regional frequency of species. Second, we will demonstrate if the two components of beta diversity have different underlying causes, here habitat attributes

Our study system is located in the western spruce–moss bioclimatic subdomain (northwestern Quebec, Canada) where fire is an important natural disturbance generating mosaics of stands with high structural and compositional heterogeneity in the forest landscape (Bergeron *et al.*, 2004). The post-fire environment is characterised by huge amounts of freshly killed and stressed trees, which are very important habitat attributes that sustain saproxylic beetles that feed directly on the bark/wood of dead and dying trees, saprophagous, mycophagous and their predators (McCullough *et al.*, 1998; Grove, 2002; Saint-Germain *et al.*, 2004; Boulanger *et al.*, 2010). The early post-fire environments are characterised by a diversity/abundance of saproxylic beetles that rapidly attack or colonize burned forests. Yet, the influence of the more-local, post-fire habitat legacies such as host-tree species, treesize, and burn-severity gradients on beta diversity patterns of saproxylic beetles is poorly known. In this study, we examine their effect on the overall, turnover and nestedness components of beta diversity by using distribution data of saproxylic beetles that actually breed in the burned forest, i.e., beetles that develop in, and emerge from, naturally fire-killed tree boles. Understanding the importance of local factors is undoubtedly crucial in conserving these essential groups, whose feeding activity may facilitate the decay of dead trees and, consequently, could have a substantial role in nutrient cycling of disturbed forests (Grove, 2002; Cobb *et al.*, 2010). On the other hand, post-fire salvage logging is currently increasing to maintain timber supply and the negative ecological consequences of this practice has become

2010; Chase *et al.*, 2011). Various indices and frameworks have been proposed to estimate beta diversity patterns independent of richness gradients (Wilson & Shmida, 1984; Koleff *et al.*, 2003). In this context, the recently proposed framework by Baselga (2010) appears to provide an improved means to systematically distinguish between the two independent components of beta diversity - the pure "spatial" turnover of species (due to species replacement) from the nestedness-driven pattern due to "species loss or gain" associated to variation in species richness. These two distinct patterns might also have different underlying causes.

The nestedness component of beta diversity essentially reflects differences in the composition of nested assemblages, i.e., when communities in species-poor sites are subsets of those found in species-rich sites (Patterson & Atmar, 1986). In other words, the difference in composition between species-poor and species-rich sites arises because the latter contains species not present in species-poor sites and not the converse. Nestedness patterns can result from mechanisms/factors that sort species hierarchically such as differential extinction or colonization rates, nested distributions or suitability of habitat gradients (Atmar & Patterson, 1993; Hylander *et al.*, 2005; Azeria *et al.*, 2006; Ulrich & Gotelli, 2007; Azeria & Kolasa, 2008). It should be emphasized that the nestedness component of beta diversity is the dissimilarity due to a nestedness effect or richness difference and not a measure of nestedness itself (Baselga, 2010).

The "spatial" turnover component of beta diversity, on the other hand, reflects distinct community ensembles and their corresponding underlying causal factors or processes , which depart notably from hierarchical effects underlying nestedness (idiosyncartic communities; Atmar & Patterson, 1993; Azeria *et al.*, 2006; 2009b). Potential causes include distinct habitats sustaining distinct communities or spatial segregations due to interspecific interactions (within same habitat template) (Loreau, 2000; Baselga, 2010; Chase *et al.*, 2011). It is clear that distinguishing between the two components is important for our understanding of biodiversity patterns. For instance, a recent study by Baselga (2010) has demonstrated that beta diversity of longhorn beetles in southern Europe was primarily caused by spatial turnover (and associated endemics due to historical effects) while in the northern Europe it was caused by spatial turnover but also by nestedness (ordered loss of species towards the north) (Baselga, 2008; Baselga, 2010). This crucial information was not evident in the un-partitioned beta diversity, which indicated a similar pattern in northern and southern Europe.

Another promising approach to decipher beta diversity patterns beyond richness gradient effects is by using null model analysis (Chase, 2007; Vellend *et al.*, 2007; Anderson *et al.*, 2011; Chase *et al.*, 2011). Null model tests have, however, been remarkably neglected in beta diversity studies (but see Chase, 2007; Chase *et al.*, 2011) despite their wide applicability in many areas of ecology and biogeography (Connor & Simberloff, 1983; Gotelli & Graves, 1996; Gotelli, 2001). Null models quantify and assess whether observed patterns depart from random expectations by comparing them against patterns emerging by randomization of observation data. This approach has been helpful in establishing non-random species cooccurrence patterns beyond random expectations (Connor & Simberloff, 1983; Gotelli & Graves, 1996; Azeria *et al.*, 2009a) and has recently been extended in "biodiversity deconstruction" for identifying species groups that are distinct in some ecological sense (Azeria *et al.*, 2009a; 2011). Similarly, null models can be useful in deciphering beta diversity patterns that are beyond (higher or lower) random chance expectations (given variation in species richness). Such beyond chance patterns may imply a deterministic causal processes; for example, Chase (2007) demonstrated that communities in drought ponds were more

2010; Chase *et al.*, 2011). Various indices and frameworks have been proposed to estimate beta diversity patterns independent of richness gradients (Wilson & Shmida, 1984; Koleff *et al.*, 2003). In this context, the recently proposed framework by Baselga (2010) appears to provide an improved means to systematically distinguish between the two independent components of beta diversity - the pure "spatial" turnover of species (due to species replacement) from the nestedness-driven pattern due to "species loss or gain" associated to variation in species

The nestedness component of beta diversity essentially reflects differences in the composition of nested assemblages, i.e., when communities in species-poor sites are subsets of those found in species-rich sites (Patterson & Atmar, 1986). In other words, the difference in composition between species-poor and species-rich sites arises because the latter contains species not present in species-poor sites and not the converse. Nestedness patterns can result from mechanisms/factors that sort species hierarchically such as differential extinction or colonization rates, nested distributions or suitability of habitat gradients (Atmar & Patterson, 1993; Hylander *et al.*, 2005; Azeria *et al.*, 2006; Ulrich & Gotelli, 2007; Azeria & Kolasa, 2008). It should be emphasized that the nestedness component of beta diversity is the dissimilarity due to a nestedness effect or richness difference and not a measure of

The "spatial" turnover component of beta diversity, on the other hand, reflects distinct community ensembles and their corresponding underlying causal factors or processes , which depart notably from hierarchical effects underlying nestedness (idiosyncartic communities; Atmar & Patterson, 1993; Azeria *et al.*, 2006; 2009b). Potential causes include distinct habitats sustaining distinct communities or spatial segregations due to interspecific interactions (within same habitat template) (Loreau, 2000; Baselga, 2010; Chase *et al.*, 2011). It is clear that distinguishing between the two components is important for our understanding of biodiversity patterns. For instance, a recent study by Baselga (2010) has demonstrated that beta diversity of longhorn beetles in southern Europe was primarily caused by spatial turnover (and associated endemics due to historical effects) while in the northern Europe it was caused by spatial turnover but also by nestedness (ordered loss of species towards the north) (Baselga, 2008; Baselga, 2010). This crucial information was not evident in the un-partitioned beta

Another promising approach to decipher beta diversity patterns beyond richness gradient effects is by using null model analysis (Chase, 2007; Vellend *et al.*, 2007; Anderson *et al.*, 2011; Chase *et al.*, 2011). Null model tests have, however, been remarkably neglected in beta diversity studies (but see Chase, 2007; Chase *et al.*, 2011) despite their wide applicability in many areas of ecology and biogeography (Connor & Simberloff, 1983; Gotelli & Graves, 1996; Gotelli, 2001). Null models quantify and assess whether observed patterns depart from random expectations by comparing them against patterns emerging by randomization of observation data. This approach has been helpful in establishing non-random species cooccurrence patterns beyond random expectations (Connor & Simberloff, 1983; Gotelli & Graves, 1996; Azeria *et al.*, 2009a) and has recently been extended in "biodiversity deconstruction" for identifying species groups that are distinct in some ecological sense (Azeria *et al.*, 2009a; 2011). Similarly, null models can be useful in deciphering beta diversity patterns that are beyond (higher or lower) random chance expectations (given variation in species richness). Such beyond chance patterns may imply a deterministic causal processes; for example, Chase (2007) demonstrated that communities in drought ponds were more

diversity, which indicated a similar pattern in northern and southern Europe.

richness. These two distinct patterns might also have different underlying causes.

nestedness itself (Baselga, 2010).

similar than expected by random chance, probably due to a deterministic filtering effect caused by a harsh environment (also see Chase *et al.*, 2011).

Clearly, a common ground between the partitioning framework and null model analysis of beta diversity needs to be established in order to promote the state of the art and practical understanding of beta diversity patterns. The idea behind null model analysis is to disentangle beta diversity patterns beyond richness gradients, which should be reflected in the turnover component of beta diversity. The nestedness component of beta diversity, on the other hand, is primarily driven by richness differences, which may be related to beta diversity expected under null distribution given richness gradients.

The present study examines the overall beta diversity components as well as the turnover and nestedness-driven components of beta diversity of saproxylic beetles emerging from tree boles following forest fire. In the present study, we have two major goals. First, we will establish the relationship of the overall beta diversity as well as its two components, turnover and nestedness, with beta diversity expected under and beyond random assembly of communities. We partition or deconstruct the overall beta diversity pattern following the partitioning framework of Baselga (2010). We perform null model analyses for beta diversity with a similar methodological basis as in Chase *et al.* (2011) to estimate the extent of beta diversity expected under and beyond random assembly, given variations in species richness among sites and the regional frequency of species. Second, we will demonstrate if the two components of beta diversity have different underlying causes, here habitat attributes defined primarily by tree-species, burn-severity, and tree-size classes.

Our study system is located in the western spruce–moss bioclimatic subdomain (northwestern Quebec, Canada) where fire is an important natural disturbance generating mosaics of stands with high structural and compositional heterogeneity in the forest landscape (Bergeron *et al.*, 2004). The post-fire environment is characterised by huge amounts of freshly killed and stressed trees, which are very important habitat attributes that sustain saproxylic beetles that feed directly on the bark/wood of dead and dying trees, saprophagous, mycophagous and their predators (McCullough *et al.*, 1998; Grove, 2002; Saint-Germain *et al.*, 2004; Boulanger *et al.*, 2010). The early post-fire environments are characterised by a diversity/abundance of saproxylic beetles that rapidly attack or colonize burned forests. Yet, the influence of the more-local, post-fire habitat legacies such as host-tree species, treesize, and burn-severity gradients on beta diversity patterns of saproxylic beetles is poorly known. In this study, we examine their effect on the overall, turnover and nestedness components of beta diversity by using distribution data of saproxylic beetles that actually breed in the burned forest, i.e., beetles that develop in, and emerge from, naturally fire-killed tree boles. Understanding the importance of local factors is undoubtedly crucial in conserving these essential groups, whose feeding activity may facilitate the decay of dead trees and, consequently, could have a substantial role in nutrient cycling of disturbed forests (Grove, 2002; Cobb *et al.*, 2010). On the other hand, post-fire salvage logging is currently increasing to maintain timber supply and the negative ecological consequences of this practice has become an emerging ecological issue in forest management (Lindenmayer *et al.*, 2008).

## **2. Material and methods**

#### **2.1 Study area and bole sampling**

The study was conducted in 72 sites within four forest burns (burned in 2005) within the western spruce–moss bioclimatic subdomain of northwestern Quebec, Canada (49o15'-50 o

Towards a Better Understanding of Beta Diversity: Deconstructing

**3.1 Deconstruction or partitioning of beta diversity** 

null model.

**3. Statistical analysis** 

Development-Team, 2010).

**3.2 Null model analysis of beta diversity** 

Composition Patterns of Saproxylic Beetles Breeding in Recently Burnt Boreal Forest 79

omitted them as such species often tend to have an unduly large influence in multivariate analysis and, consequently, distort the overall pattern (e.g. in ordinations). In addition, it is advisable that species that are uncharacteristic of the species pool (here the species might be unlikely to colonize fire-killed tress in burned forests) are omitted to minimize their effect in null model analysis of beta diversity (see Chase et al 2011 for discussion related to the latter issue). We thus restrict our analysis to 32 species that occurred in more than two sites. These species accounted for 94% of the presence–absence matrix and for 97% of the total abundance (1639 individuals). We also omitted one site that only had a single species, which was recorded in nearly all sites (71 of 72 sites, thus its composition dissimilarity was always zero with respect to this site) as it was not possible to generate a "null community" and consequently "null" beta diversity values for the site given the constraints imposed in the

Our first goal was to partition the overall beta diversity pattern of the saproxylic beetles into two components: the "spatial" turnover and nestedness component. We also wanted to examine their relationship to beta diversity values expected *under* and *beyond* random community assembly given a null model. We emphasize that throughout the paper, beta diversity refers to site-to-site composition dissimilarity. We partitioned beta diversity of saproxylic beetles following the framework of Baselga (2010) as: βsor = βsim + βnes; where βsor (Sorensen dissimilarity) represents the total difference in species composition between two sites, and βsim (Simpson dissimilarity) and βnes (nestedness-driven dissimilarity) are its "turnover" and 'nestedness" components, respectively. We computed these components using the function provided by Baselga (2010), as implemented in R version 2.11 (R-

We used null models to disentangle beta diversity values expected under null distributions and beyond random assembly given constraints set by null models ("delta" beta diversity; also see Chase et al 2011). In principle, the null model analysis for site-to-site beta diversity is methodologically the same as that used for species-to-species co-occurrence patterns in

First, we computed the beta diversity for the observed data using the Sorensen dissimilarity index (βsor), consistent with beta diversity partitioning framework. Second, null models are applied to randomize the observation data to generate "null" communities (n=1000) for which beta diversity will be calculated. We used two null models: the Fixed-Fixed (FF) and Fixed-Equiprobable (FE) null models. Both null models maintain species richness of sites from the observation matrix. The FF null model also maintains species frequency as in the observation data, while FE null model sample species from the regional species pool equiprobably. We used the function *permatfull*, a wrapper for *commsimulator*, in the R-package Vegan (Oksanen *et al.*, 2010) to generate 1000 null matrices according each null model. For the FF null model, we used the quasi swap algorithm (Miklós & Podani, 2004), which generates matrices that are

We computed beta diversity (using the Sorensen dissimilarity index) for each of the 1000 null matrices (βsor-null.mat), from which the beta diversity expected by "random" chance was

the context of a "biodiversity deconstruction" framework (Azeria *et al.*, 2009a).

independent of each other and different from the original matrix.

40'N and 75 o 00'W-73 o 45'W). This forest subdomain is typically dominated by black spruce (*Picea mariana*) and, to a lesser extent, jack pine (*Pinus banksiana*)*.* The forest landscape also contains various combinations of Balsam fir (*Abies balsamea*), trembling aspen (*Populus tremuloides*) and paper birch (*Betula papyrifera*). In this forest subdomain, forest fire is an important disturbance agent and occurs in a relatively short cycle (120-180 years) and, as consequence, the landscape is dominated by even-aged forest stands (Bergeron *et al.*, 2004).

We sampled sites using a systematic factorial design: 2 tree species X 3 burn levels X 4 tree size categories X 3 replicates (Total = 72 forest sites). Our focus was to examine saproxylic beetles that actually develop in, and emerge from, fire-killed trees following forest fire. Therefore, from each sampling site, we retrieved a 50-cm bole section from five fire-killed trees that we felled (Total 360 trees). This was carried in early June 2006, one year after fire, and when initial colonization or attack by saproxylic beetles was achieved naturally (Boulanger and Sirois 2007). The retrieved bole sections were then enclosed in rearing (emergence) cages placed in a field insectarium, where they were suspended while respecting their natural vertical orientation. The offspring were bred out, and the emerging adults/larvae were collected monthly from June to November in 2006 and 2007 in a vial (with preservative) placed under the bole section. Vouchers are conserved in the insect collection of the Laurentian Forestry Centre. All collected specimens were identified to the lowest taxonomic level (species or genus level) whenever possible, otherwise they were identified only to the family level. The "species" lists of boles were pooled per-site in subsequent beta diversity analysis.

Our sampling protocol thus enabled us to investigate simultaneously the effects of variation in tree species (black spruce and jack pine), burn severity (low, moderate and high) and tree size (db1=8-12 cm; db2=12-16cm; db3=16-20 cm; db4=20-24 cm) for beta diversity of these saproxylic beetles. The tree size classes were based on diameter at breast height (dbh); and burn severity classes were visually defined following criteria adopted by *Ministère des Ressources naturelles de la Faune du Québec* (MRNF) but with modifications of the classes. These habitat variables have been shown to differentially influence the distribution of saproxylic beetles in burned forests (Saint-Germain *et al.*, 2004; Boulanger *et al.*, 2010); and we expect them to also differentially contribute to the turnover and nestedness component of beta diversity. For example, different tree species might attract or host different beetle species and, consequently, might contribute more to the turnover than nestedness component of beta diversity. The effects of burn severity and tree size class is to generally create a gradient of habitat suitability, the lowest suitability being in small trees and/or severe-burns and highest suitability being in large trees of low severity classes (Boulanger *et al.*, 2010). Such a hierarchical suitability gradient would potentially contribute to the nestedness component of beta diversity. On the other hand, the within habitat class contribution for beta diversity components might differ; for example, there could be more dispersion or turnover withinsmall than within-large tree size classes due to intensive competition in smaller trees that cause segregated species co-occurrences, or some species being biased towards species-poor sites, termed as idiosyncratic species (Azeria *et al.*, 2006; 2009b).

#### **2.2 Saproxylic beetle (emerging) distribution data**

We compiled a presence–absence matrix (sites x species) of all saproxylic beetles in each sampling site by pooling the respective data from five fire-killed tree boles. The total number of "species" recorded was 62 species. Nearly half of these species (30 species) were very rare recorded from only one or two sites, and many were primarily common in dead trees in unburnt forest such as those generated by gap dynamics (Boucher, 2011). We omitted them as such species often tend to have an unduly large influence in multivariate analysis and, consequently, distort the overall pattern (e.g. in ordinations). In addition, it is advisable that species that are uncharacteristic of the species pool (here the species might be unlikely to colonize fire-killed tress in burned forests) are omitted to minimize their effect in null model analysis of beta diversity (see Chase et al 2011 for discussion related to the latter issue). We thus restrict our analysis to 32 species that occurred in more than two sites. These species accounted for 94% of the presence–absence matrix and for 97% of the total abundance (1639 individuals). We also omitted one site that only had a single species, which was recorded in nearly all sites (71 of 72 sites, thus its composition dissimilarity was always zero with respect to this site) as it was not possible to generate a "null community" and consequently "null" beta diversity values for the site given the constraints imposed in the null model.

## **3. Statistical analysis**

78 Research in Biodiversity – Models and Applications

40'N and 75 o 00'W-73 o 45'W). This forest subdomain is typically dominated by black spruce (*Picea mariana*) and, to a lesser extent, jack pine (*Pinus banksiana*)*.* The forest landscape also contains various combinations of Balsam fir (*Abies balsamea*), trembling aspen (*Populus tremuloides*) and paper birch (*Betula papyrifera*). In this forest subdomain, forest fire is an important disturbance agent and occurs in a relatively short cycle (120-180 years) and, as consequence, the landscape is dominated by even-aged forest stands (Bergeron *et al.*, 2004). We sampled sites using a systematic factorial design: 2 tree species X 3 burn levels X 4 tree size categories X 3 replicates (Total = 72 forest sites). Our focus was to examine saproxylic beetles that actually develop in, and emerge from, fire-killed trees following forest fire. Therefore, from each sampling site, we retrieved a 50-cm bole section from five fire-killed trees that we felled (Total 360 trees). This was carried in early June 2006, one year after fire, and when initial colonization or attack by saproxylic beetles was achieved naturally (Boulanger and Sirois 2007). The retrieved bole sections were then enclosed in rearing (emergence) cages placed in a field insectarium, where they were suspended while respecting their natural vertical orientation. The offspring were bred out, and the emerging adults/larvae were collected monthly from June to November in 2006 and 2007 in a vial (with preservative) placed under the bole section. Vouchers are conserved in the insect collection of the Laurentian Forestry Centre. All collected specimens were identified to the lowest taxonomic level (species or genus level) whenever possible, otherwise they were identified only to the family level. The "species"

Our sampling protocol thus enabled us to investigate simultaneously the effects of variation in tree species (black spruce and jack pine), burn severity (low, moderate and high) and tree size (db1=8-12 cm; db2=12-16cm; db3=16-20 cm; db4=20-24 cm) for beta diversity of these saproxylic beetles. The tree size classes were based on diameter at breast height (dbh); and burn severity classes were visually defined following criteria adopted by *Ministère des Ressources naturelles de la Faune du Québec* (MRNF) but with modifications of the classes. These habitat variables have been shown to differentially influence the distribution of saproxylic beetles in burned forests (Saint-Germain *et al.*, 2004; Boulanger *et al.*, 2010); and we expect them to also differentially contribute to the turnover and nestedness component of beta diversity. For example, different tree species might attract or host different beetle species and, consequently, might contribute more to the turnover than nestedness component of beta diversity. The effects of burn severity and tree size class is to generally create a gradient of habitat suitability, the lowest suitability being in small trees and/or severe-burns and highest suitability being in large trees of low severity classes (Boulanger *et al.*, 2010). Such a hierarchical suitability gradient would potentially contribute to the nestedness component of beta diversity. On the other hand, the within habitat class contribution for beta diversity components might differ; for example, there could be more dispersion or turnover withinsmall than within-large tree size classes due to intensive competition in smaller trees that cause segregated species co-occurrences, or some species being biased towards species-poor sites,

We compiled a presence–absence matrix (sites x species) of all saproxylic beetles in each sampling site by pooling the respective data from five fire-killed tree boles. The total number of "species" recorded was 62 species. Nearly half of these species (30 species) were very rare recorded from only one or two sites, and many were primarily common in dead trees in unburnt forest such as those generated by gap dynamics (Boucher, 2011). We

lists of boles were pooled per-site in subsequent beta diversity analysis.

termed as idiosyncratic species (Azeria *et al.*, 2006; 2009b).

**2.2 Saproxylic beetle (emerging) distribution data** 

#### **3.1 Deconstruction or partitioning of beta diversity**

Our first goal was to partition the overall beta diversity pattern of the saproxylic beetles into two components: the "spatial" turnover and nestedness component. We also wanted to examine their relationship to beta diversity values expected *under* and *beyond* random community assembly given a null model. We emphasize that throughout the paper, beta diversity refers to site-to-site composition dissimilarity. We partitioned beta diversity of saproxylic beetles following the framework of Baselga (2010) as: βsor = βsim + βnes; where βsor (Sorensen dissimilarity) represents the total difference in species composition between two sites, and βsim (Simpson dissimilarity) and βnes (nestedness-driven dissimilarity) are its "turnover" and 'nestedness" components, respectively. We computed these components using the function provided by Baselga (2010), as implemented in R version 2.11 (R-Development-Team, 2010).

#### **3.2 Null model analysis of beta diversity**

We used null models to disentangle beta diversity values expected under null distributions and beyond random assembly given constraints set by null models ("delta" beta diversity; also see Chase et al 2011). In principle, the null model analysis for site-to-site beta diversity is methodologically the same as that used for species-to-species co-occurrence patterns in the context of a "biodiversity deconstruction" framework (Azeria *et al.*, 2009a).

First, we computed the beta diversity for the observed data using the Sorensen dissimilarity index (βsor), consistent with beta diversity partitioning framework. Second, null models are applied to randomize the observation data to generate "null" communities (n=1000) for which beta diversity will be calculated. We used two null models: the Fixed-Fixed (FF) and Fixed-Equiprobable (FE) null models. Both null models maintain species richness of sites from the observation matrix. The FF null model also maintains species frequency as in the observation data, while FE null model sample species from the regional species pool equiprobably. We used the function *permatfull*, a wrapper for *commsimulator*, in the R-package Vegan (Oksanen *et al.*, 2010) to generate 1000 null matrices according each null model. For the FF null model, we used the quasi swap algorithm (Miklós & Podani, 2004), which generates matrices that are independent of each other and different from the original matrix.

We computed beta diversity (using the Sorensen dissimilarity index) for each of the 1000 null matrices (βsor-null.mat), from which the beta diversity expected by "random" chance was

Towards a Better Understanding of Beta Diversity: Deconstructing

βsor was higher than βsor-null.mat.

R-package Vegan (Oksanen *et al.*, 2010).

**diversity** 

Composition Patterns of Saproxylic Beetles Breeding in Recently Burnt Boreal Forest 81

true for more dissimilar site pairs. Note that the index is based on dissimilarity; the tests for statistical significance of "more dissimilar than expected" are made by taking the complement of the transformed values (for index values 0.5 to 1, subtract them from 1). Alternatively, to test statistical significance, one may use the proportion of null matrices in which βsor-null.mat is the same or larger (or lower) βsor for obtaining higher (lower) dissimilarity between sites. In the latter case, a comparable index with high βRC values (closer to 1, which was based on shared species) will be the number of matrices for which

To assess the relationships between the different algorithms, we examined the probability

Our second goal was to examine the contribution of a set of habitat variables (tree species, burn severity classes, tree size classes) as well as their interactions for the overall (βsor), turnover (βsim) and nestedness-driven (βnes) beta diversity of saproxylic assemblages. In addition, we examined the spatial effect on composition dissimilarity by considering its correlation with site-to-site geographical distance, and then by considering the dissimilarity between and within the forest burns (four burns). We used a nonparametric, Permutatitional Multivariate Analysis of Variance (PERMANOVA; Anderson, 2001) to test whether categories for each habitat factor differed in their variability in species composition or beta diversity (based on βsor, βsim and βnes). The method computes a pseudo F-ratio, which is the ratio of composition dissimilarity within a treatment (habitat class) to that of between treatments, and then tests its significance by permutation (9999 replicates). Accordingly, significant results might indicate differences in dissimilarity among the classes (difference between treatment's centroid in multivariate space), due to differences in the within-class dispersion (i.e., mean distances of members to their group centroid) or both. The potential role of each is distinguished by running a complementary analysis, the Permutational Analysis of Multivariate Dispersions (PERMDISP; Anderson *et al.*, 2006), which tests whether classes (treatments) differed in their within-treatment dispersion or beta diversity. As an example, consider tree species, Black spruce (BSP) and Jack pine (JPI), as sources of variability. A significant result by PERMANOVA and a non-significant difference by PERMDISP would suggest differences in their across class dissimilarity (BSP ≠ JPI), i.e., the treatments differ in their centroid in multivariate space and not in the within-treatment dispersion (BSP-BSP = JPI-JPI dispersion). When significant results are also obtained by PERMDISP, one may run pairwise tests to examine which of the classes had higher dispersion (particularly when dealing with more than two classes). We performed PERMANOVA and PERMDISP using the functions *adonis* and *betadisper*, respectively, in the

We also applied multivariate regression analysis on distance matrices (MRM) (Zapala & Schork, 2006; Lichstein, 2007), which might help model the beta diversity variation of within- and between habitat classes, as well as the difference between them, simultaneously. MRM is essentially an ANOVA-like analysis performed on the site-to-site dissimilarity matrix expanded into a vector and modelled against the corresponding contrast of classes of habitat (for each habitat variable) changed into vector. For example, when considering the effect of tree species, the contrasts [within black spruce (BSP-BSP), within Jack pine (JPI-JPI) and between the two species (BSP-JPI)] are unfolded into vector. Note that according to our

computed for βsor-SES (based on the FF and FE null models) and that of βRC metric.

**3.3 Analysis of habitat effects on overall beta diversity and components of beta** 

estimated by computing their mean (βsor-null) and standard deviation (βsor-null.sd). These values effectively represent the beta diversity expected by random chance, given differences in richness among sites (and species frequencies for FF null model). Third, we estimated beta diversity independent of and beyond random chance (βsor-diff) by computing the difference in beta diversity values between the observation data (βsor) and null communities (βsor-null), i.e. βsor-diff = βsor - βsor-null. The difference could also be expressed by effect sizes (or standard deviation units) as βsor-SES = βsor-diff /βsor-null.sd. The βsor-SES measures the number of standard deviations that the observed dissimilarity βsor is above or below the mean index of the dissimilarity obtained in the null distribution (βsor-null). The value of βsor-diff (and βsor-SES) will be positive for sites pairs that are more dissimilar than expected and negative for those less dissimilar (more similar) than expected. For ordinations or other graphical representations, the βsor-SES can be rescaled using the ranging formula as (βsor-SES - βsor-SES.min)/ (βsor-SES.max - βsor-SES.min), where βsor-SES.min and βsor-SES.max are the minimum and maximum βsor-SES values. This will rescale the βsor-SES between 0 (for sites that are more similar) and 1 (for sites that are more dissimilar) in the same manner as applied for species-pairs in Azeria et al. (2009a; 2011).

Our main goal was to link the beta diversity components with that expected under null models. It is important that all the terms are expressed in the same units; therefore we examined how the overall (βsor), and its partition to "spatial" turnover (βsim) and nestedness (βnes) components were related to that expected under the null model (βsor-null) and beyond (βsor-diff).

We also applied the null model recently proposed by Chase et al (2011; also see Chase, 2007) which is a modification of the Raup-Crick metric (βRC) (Raup & Crick, 1979). The null model for computing βRC maintains species richness of site pairs as in the observation data, and species are sampled *proportional* to their frequency, rarity and commonness. The original implementation of the metric (Raup & Crick, 1979) also maintains for species richness, but it does not constrain for species incidence, i.e., it assumes species would be sampled equiprobably. The latter option is also available if desired. We computed βRC using the Rfunction provided by Chase et al (2011). The program computes a βRC that shows the probability that the observed site-to-site dissimilarity is by chance by counting the number of null matrices where observed shared species is less than simulated. The computed probabilities (between 0 and 1) are rescaled by subtracting -0.5 and multiplying by 2 into values between -1 (less dissimilar) and 1 (more dissimilar). Values around 0 are not different from expected by chance (for other details see Chase *et al.*, 2011).

It is worth mentioning the similarities and differences in our implementation of the null model from that of Chase et al. (2011). Our implementation of the FF null model maintains the species frequency *exactly* as in the observed data, while the null model associated for βRC samples species *proportional* to their incidence in the data set. On the other hand, our implementation of the FE null model will be effectively similar to that of βRC when species are sampled equiprobably, as in the original formulation of Raup-Crick metric (Raup & Crick, 1979). The other difference is, the βsor-diff indicates actual differences in dissimilarity values between observed and null communities in beta diversity-units, while βRC expresses the difference (for the corresponding null model) in terms of a probability index. The probability associated with βsor-diff can be estimated by applying inverse-logit (R function *plogis*) or normal (R function *pnorm* ) transformations on its standardized form, the βsor-SES. Transformed values will be closer to zero (negative βsor-SES) for site pairs that are less dissimilar (more similar) than expected by random chance (one tail), and the converse is

estimated by computing their mean (βsor-null) and standard deviation (βsor-null.sd). These values effectively represent the beta diversity expected by random chance, given differences in richness among sites (and species frequencies for FF null model). Third, we estimated beta diversity independent of and beyond random chance (βsor-diff) by computing the difference in beta diversity values between the observation data (βsor) and null communities (βsor-null), i.e. βsor-diff = βsor - βsor-null. The difference could also be expressed by effect sizes (or standard deviation units) as βsor-SES = βsor-diff /βsor-null.sd. The βsor-SES measures the number of standard deviations that the observed dissimilarity βsor is above or below the mean index of the dissimilarity obtained in the null distribution (βsor-null). The value of βsor-diff (and βsor-SES) will be positive for sites pairs that are more dissimilar than expected and negative for those less dissimilar (more similar) than expected. For ordinations or other graphical representations, the βsor-SES can be rescaled using the ranging formula as (βsor-SES - βsor-SES.min)/ (βsor-SES.max - βsor-SES.min), where βsor-SES.min and βsor-SES.max are the minimum and maximum βsor-SES values. This will rescale the βsor-SES between 0 (for sites that are more similar) and 1 (for sites that are more dissimilar) in the same manner as applied for species-pairs in Azeria et

Our main goal was to link the beta diversity components with that expected under null models. It is important that all the terms are expressed in the same units; therefore we examined how the overall (βsor), and its partition to "spatial" turnover (βsim) and nestedness (βnes) components were related to that expected under the null model (βsor-null) and beyond

We also applied the null model recently proposed by Chase et al (2011; also see Chase, 2007) which is a modification of the Raup-Crick metric (βRC) (Raup & Crick, 1979). The null model for computing βRC maintains species richness of site pairs as in the observation data, and species are sampled *proportional* to their frequency, rarity and commonness. The original implementation of the metric (Raup & Crick, 1979) also maintains for species richness, but it does not constrain for species incidence, i.e., it assumes species would be sampled equiprobably. The latter option is also available if desired. We computed βRC using the Rfunction provided by Chase et al (2011). The program computes a βRC that shows the probability that the observed site-to-site dissimilarity is by chance by counting the number of null matrices where observed shared species is less than simulated. The computed probabilities (between 0 and 1) are rescaled by subtracting -0.5 and multiplying by 2 into values between -1 (less dissimilar) and 1 (more dissimilar). Values around 0 are not different

It is worth mentioning the similarities and differences in our implementation of the null model from that of Chase et al. (2011). Our implementation of the FF null model maintains the species frequency *exactly* as in the observed data, while the null model associated for βRC samples species *proportional* to their incidence in the data set. On the other hand, our implementation of the FE null model will be effectively similar to that of βRC when species are sampled equiprobably, as in the original formulation of Raup-Crick metric (Raup & Crick, 1979). The other difference is, the βsor-diff indicates actual differences in dissimilarity values between observed and null communities in beta diversity-units, while βRC expresses the difference (for the corresponding null model) in terms of a probability index. The probability associated with βsor-diff can be estimated by applying inverse-logit (R function *plogis*) or normal (R function *pnorm* ) transformations on its standardized form, the βsor-SES. Transformed values will be closer to zero (negative βsor-SES) for site pairs that are less dissimilar (more similar) than expected by random chance (one tail), and the converse is

from expected by chance (for other details see Chase *et al.*, 2011).

al. (2009a; 2011).

(βsor-diff).

true for more dissimilar site pairs. Note that the index is based on dissimilarity; the tests for statistical significance of "more dissimilar than expected" are made by taking the complement of the transformed values (for index values 0.5 to 1, subtract them from 1). Alternatively, to test statistical significance, one may use the proportion of null matrices in which βsor-null.mat is the same or larger (or lower) βsor for obtaining higher (lower) dissimilarity between sites. In the latter case, a comparable index with high βRC values (closer to 1, which was based on shared species) will be the number of matrices for which βsor was higher than βsor-null.mat.

To assess the relationships between the different algorithms, we examined the probability computed for βsor-SES (based on the FF and FE null models) and that of βRC metric.

#### **3.3 Analysis of habitat effects on overall beta diversity and components of beta diversity**

Our second goal was to examine the contribution of a set of habitat variables (tree species, burn severity classes, tree size classes) as well as their interactions for the overall (βsor), turnover (βsim) and nestedness-driven (βnes) beta diversity of saproxylic assemblages. In addition, we examined the spatial effect on composition dissimilarity by considering its correlation with site-to-site geographical distance, and then by considering the dissimilarity between and within the forest burns (four burns). We used a nonparametric, Permutatitional Multivariate Analysis of Variance (PERMANOVA; Anderson, 2001) to test whether categories for each habitat factor differed in their variability in species composition or beta diversity (based on βsor, βsim and βnes). The method computes a pseudo F-ratio, which is the ratio of composition dissimilarity within a treatment (habitat class) to that of between treatments, and then tests its significance by permutation (9999 replicates). Accordingly, significant results might indicate differences in dissimilarity among the classes (difference between treatment's centroid in multivariate space), due to differences in the within-class dispersion (i.e., mean distances of members to their group centroid) or both. The potential role of each is distinguished by running a complementary analysis, the Permutational Analysis of Multivariate Dispersions (PERMDISP; Anderson *et al.*, 2006), which tests whether classes (treatments) differed in their within-treatment dispersion or beta diversity. As an example, consider tree species, Black spruce (BSP) and Jack pine (JPI), as sources of variability. A significant result by PERMANOVA and a non-significant difference by PERMDISP would suggest differences in their across class dissimilarity (BSP ≠ JPI), i.e., the treatments differ in their centroid in multivariate space and not in the within-treatment dispersion (BSP-BSP = JPI-JPI dispersion). When significant results are also obtained by PERMDISP, one may run pairwise tests to examine which of the classes had higher dispersion (particularly when dealing with more than two classes). We performed PERMANOVA and PERMDISP using the functions *adonis* and *betadisper*, respectively, in the R-package Vegan (Oksanen *et al.*, 2010).

We also applied multivariate regression analysis on distance matrices (MRM) (Zapala & Schork, 2006; Lichstein, 2007), which might help model the beta diversity variation of within- and between habitat classes, as well as the difference between them, simultaneously. MRM is essentially an ANOVA-like analysis performed on the site-to-site dissimilarity matrix expanded into a vector and modelled against the corresponding contrast of classes of habitat (for each habitat variable) changed into vector. For example, when considering the effect of tree species, the contrasts [within black spruce (BSP-BSP), within Jack pine (JPI-JPI) and between the two species (BSP-JPI)] are unfolded into vector. Note that according to our

Towards a Better Understanding of Beta Diversity: Deconstructing

*proportional* to their regional frequency in the observed data.

Composition Patterns of Saproxylic Beetles Breeding in Recently Burnt Boreal Forest 83

the beta diversity consequent upon and beyond richness gradients respectively, our results suggest that both components might be conditioned by species frequency (commonness and rarity) patterns. This also implies that the appropriate null model for beta diversity analysis *sensu lato* should preserve not only observed species richness and but also the species frequency. The importance for null models to also control for species commonness and rarity has been emphasized in the modified Raup-Crick metric (βRC) proposed by Chase *et al*. (2011). They modified the original Raup-Crick probabilistic index that samples species equiprobably (Raup & Crick, 1979) by implementing a null model that samples species

Fig. 1. The relationship of overall betadiversity (βsor ), and its partitions, the "spatial" turnover (βsim) and nestedness (βnes) components, with beta diversity values expected under

null distributions (βsor-null) and beyond-null model expectations (βsor-diff). The null communities were generated using fixed-fixed null model (FF null model).

construction, the within comparisons (BSP-BSP and JPI-JPI) are regarded as distinct values or factors and indexed accordingly in the explanatory variable. This re-indexing of the pairwise factors is done for each explanatory variable considered here, i.e., tree species, burn severity, tress size classes and forest burns where the sampling sites were located. The significance of MRM is tested by permutation of the distance vector. Note that this provides a computationally efficient way to compute an equivalent permutation test if the raw data of species composition were permuted and distance was computed afterwards (Hayden *et al.*, 2009). We will focus mainly on the results of pairwise tests among the "new classes" (e.g., BSP-BSP, BSP-JPI and JPI-JPI) of the MRM with that of PERMANOVA and PERMDISP to further unravel the source of variation for the overall beta diversity and beta diversity components.

#### **4 Results and discussion**

#### **4.1 The relationship of beta diversity components and expectations under and beyond null models**

We present results (Fig. 1 and Fig. 2) that illustrate the relationship between overall, observed beta diversity (βsor), and more specifically its turnover (βsim) and nestedness (βnes) components (sensu Baselga, 2010) with beta a diversity pattern expected by *"random"* (βsornull) and with *deviations beyond* (βsor-diff) random distribution of species among sites using two null models. The first null model preserves both species richness and species incidence (FF null model; Fig. 1), while the second null model preserves only species richness (FE null model; Fig. 2). While the overall beta diversity (βsor) value was generally correlated to that expected by *random* (βsor-null) under FF (Fig.1a) and FE (Fig.1a) null models, a stronger relationship was evident with that of *deviations beyond* null model expectations, i.e., βsor-diff (FF: Fig.1b; FE: Fig. 2b). This indicates that beta diversity patterns of saproxylic beetles show a strong signal of deterministic underlying processes that deviates from random expectations given variations in species richness and species incidence.

A more explicit examination of the beta diversity components pattern via partitioning into its turnover (βsim) and nestedness (βnes) components (sensu Baselga, 2010) revealed that the two components exhibited distinct relationships to beta diversity values expected under (βsor-null) and beyond (βsor-diff) null models. Thus, the turnover component was linearly related to beta diversity deviating *beyond* "random" expectations (Figs 1d and 2d), while the nestedness component of beta diversity was related to the *null* expectations (Figs. 1e and 2e). The converse relationships, i.e., turnover with expectations under a null model (Figs. 1c and 2c) and that of nestedness components with deviations beyond "random" (Figs. 1f and 2f) were not significant. Taken together our results demonstrate quantitatively (via null models) that the beta diversity components (sensu Baselga , 2010) do indeed reflect different extents of dependence on richness gradients. More specifically, we show the independence of the turnover and the dependence of the nestedness component on richness variations. In addition, our result emphasizes that the nestedness-driven component (sensu Baselga, 2010) should be interpreted as a general effect of richness difference on beta diversity.

Moreover, our results indicate that the relationship of turnover and nestedness components with respective null model expectations, i.e., βsor-diff and βsor-null, were stronger when the null model also preserved the species incidence (FF-null model: Fig. 1d, r=0.84; Fig.1e, r=0.73) than when species were assumed to be sampled equiprobably (FE-null model: Fig. 2d, r=0.68; Fig.2e, r=0.44). Thus, while the nestedness and turnover components would reflect

construction, the within comparisons (BSP-BSP and JPI-JPI) are regarded as distinct values or factors and indexed accordingly in the explanatory variable. This re-indexing of the pairwise factors is done for each explanatory variable considered here, i.e., tree species, burn severity, tress size classes and forest burns where the sampling sites were located. The significance of MRM is tested by permutation of the distance vector. Note that this provides a computationally efficient way to compute an equivalent permutation test if the raw data of species composition were permuted and distance was computed afterwards (Hayden *et al.*, 2009). We will focus mainly on the results of pairwise tests among the "new classes" (e.g., BSP-BSP, BSP-JPI and JPI-JPI) of the MRM with that of PERMANOVA and PERMDISP to further unravel the source of variation for the overall beta diversity and beta diversity

**4.1 The relationship of beta diversity components and expectations under and** 

expectations given variations in species richness and species incidence.

should be interpreted as a general effect of richness difference on beta diversity.

Moreover, our results indicate that the relationship of turnover and nestedness components with respective null model expectations, i.e., βsor-diff and βsor-null, were stronger when the null model also preserved the species incidence (FF-null model: Fig. 1d, r=0.84; Fig.1e, r=0.73) than when species were assumed to be sampled equiprobably (FE-null model: Fig. 2d, r=0.68; Fig.2e, r=0.44). Thus, while the nestedness and turnover components would reflect

We present results (Fig. 1 and Fig. 2) that illustrate the relationship between overall, observed beta diversity (βsor), and more specifically its turnover (βsim) and nestedness (βnes) components (sensu Baselga, 2010) with beta a diversity pattern expected by *"random"* (βsornull) and with *deviations beyond* (βsor-diff) random distribution of species among sites using two null models. The first null model preserves both species richness and species incidence (FF null model; Fig. 1), while the second null model preserves only species richness (FE null model; Fig. 2). While the overall beta diversity (βsor) value was generally correlated to that expected by *random* (βsor-null) under FF (Fig.1a) and FE (Fig.1a) null models, a stronger relationship was evident with that of *deviations beyond* null model expectations, i.e., βsor-diff (FF: Fig.1b; FE: Fig. 2b). This indicates that beta diversity patterns of saproxylic beetles show a strong signal of deterministic underlying processes that deviates from random

A more explicit examination of the beta diversity components pattern via partitioning into its turnover (βsim) and nestedness (βnes) components (sensu Baselga, 2010) revealed that the two components exhibited distinct relationships to beta diversity values expected under (βsor-null) and beyond (βsor-diff) null models. Thus, the turnover component was linearly related to beta diversity deviating *beyond* "random" expectations (Figs 1d and 2d), while the nestedness component of beta diversity was related to the *null* expectations (Figs. 1e and 2e). The converse relationships, i.e., turnover with expectations under a null model (Figs. 1c and 2c) and that of nestedness components with deviations beyond "random" (Figs. 1f and 2f) were not significant. Taken together our results demonstrate quantitatively (via null models) that the beta diversity components (sensu Baselga , 2010) do indeed reflect different extents of dependence on richness gradients. More specifically, we show the independence of the turnover and the dependence of the nestedness component on richness variations. In addition, our result emphasizes that the nestedness-driven component (sensu Baselga, 2010)

components.

**4 Results and discussion** 

**beyond null models** 

the beta diversity consequent upon and beyond richness gradients respectively, our results suggest that both components might be conditioned by species frequency (commonness and rarity) patterns. This also implies that the appropriate null model for beta diversity analysis *sensu lato* should preserve not only observed species richness and but also the species frequency. The importance for null models to also control for species commonness and rarity has been emphasized in the modified Raup-Crick metric (βRC) proposed by Chase *et al*. (2011). They modified the original Raup-Crick probabilistic index that samples species equiprobably (Raup & Crick, 1979) by implementing a null model that samples species *proportional* to their regional frequency in the observed data.

Fig. 1. The relationship of overall betadiversity (βsor ), and its partitions, the "spatial" turnover (βsim) and nestedness (βnes) components, with beta diversity values expected under null distributions (βsor-null) and beyond-null model expectations (βsor-diff). The null communities were generated using fixed-fixed null model (FF null model).

Towards a Better Understanding of Beta Diversity: Deconstructing

normal distribution.

Composition Patterns of Saproxylic Beetles Breeding in Recently Burnt Boreal Forest 85

to βRC metric, we estimated the probability associated with βsor-diff values by applying a inverse-logit transformation on the standard effect size of their deviations (standard deviation units), i.e., βsor-SES (also see Azeria *et al.*, 2009a; 2011). Accordingly, the βsor-SES computed using the FF-null model was linearly related to βRC, albeit the data points for βRC were slightly lower from the one-to-one line (Fig. 3a). This slight deviation is expected given that our FF-null model *preserves* observed species frequency while the null model for βRC samples species *proportional* to their respective observed frequencies. On the other hand, values of βRC were higher than the values computed for βsor-SES under FE-null model, which samples species equiprobably (Fig 3b). When the two metrics were computed using a comparable null models that sampled species equiprobably, the βsor-SES and βRC were highly correlated (Fig. 3c). The same results were obtained when βsor-SES was transformed by

Fig. 3. Relationship between dissimilarity/beta diversity values beyond null expectations as expressed in terms of standard effect size (βsor-SES) and the Raup-Crick metric (βRC). The βsor-SES values were transformed using logistic function. In all null models, the species richness of sites is maintained as in the observed. The species frequency was preserved in computing the βsor-SES (FF-null model: a,c), or equiprobably (FE-null model: b,d) (see method). The null models for computing the βRC sampled species proportional to observed frequencies (βRC:

a,b) or equiprobably (\*βRC: c,d) (for details of the method see Chase et al., 2001).

Fig. 2. The relationship of overall betadiversity (βsor ), and its partitions, the "spatial" turnover (βsim) and nestedness (βnes) components, with beta diversity values expected under null distributions (βsor-null) and beyond-null model expectations (βsor-diff). The null communities were generated using fixed-equiprobable null model (FE null model).

The utility of null models for disentangling beta diversity patterns independent of richness variation has only recently been emphasized (Anderson *et al.*, 2011; Chase *et al.*, 2011). Chase *et al*. (2011) have highlighted that the βRC dissimilarity metric (as expressed in probability between 0 and 1, or as rescaled from -1 to 1; see Chase et al 2011) reflects the beta diversity or dissimilarity patterns independent of richness variations (also see Chase, 2007). We also found that βRC was generally correlated withthe turnover component of beta diversity βsim (r =0.73) as well as to our computation of beta diversity beyond null expectations βsor-diff (r =0.90). However, it should be emphasized that the βRC dissimilarity metric is a *probability*  index rather than a direct measure of the departure from random chance in terms of "beta diversity units" sensu stricto (βsim or βsor-diff). To obtain a dissimilarity measure comparable

Fig. 2. The relationship of overall betadiversity (βsor ), and its partitions, the "spatial" turnover (βsim) and nestedness (βnes) components, with beta diversity values expected under

The utility of null models for disentangling beta diversity patterns independent of richness variation has only recently been emphasized (Anderson *et al.*, 2011; Chase *et al.*, 2011). Chase *et al*. (2011) have highlighted that the βRC dissimilarity metric (as expressed in probability between 0 and 1, or as rescaled from -1 to 1; see Chase et al 2011) reflects the beta diversity or dissimilarity patterns independent of richness variations (also see Chase, 2007). We also found that βRC was generally correlated withthe turnover component of beta diversity βsim (r =0.73) as well as to our computation of beta diversity beyond null expectations βsor-diff (r =0.90). However, it should be emphasized that the βRC dissimilarity metric is a *probability*  index rather than a direct measure of the departure from random chance in terms of "beta diversity units" sensu stricto (βsim or βsor-diff). To obtain a dissimilarity measure comparable

null distributions (βsor-null) and beyond-null model expectations (βsor-diff). The null communities were generated using fixed-equiprobable null model (FE null model). to βRC metric, we estimated the probability associated with βsor-diff values by applying a inverse-logit transformation on the standard effect size of their deviations (standard deviation units), i.e., βsor-SES (also see Azeria *et al.*, 2009a; 2011). Accordingly, the βsor-SES computed using the FF-null model was linearly related to βRC, albeit the data points for βRC were slightly lower from the one-to-one line (Fig. 3a). This slight deviation is expected given that our FF-null model *preserves* observed species frequency while the null model for βRC samples species *proportional* to their respective observed frequencies. On the other hand, values of βRC were higher than the values computed for βsor-SES under FE-null model, which samples species equiprobably (Fig 3b). When the two metrics were computed using a comparable null models that sampled species equiprobably, the βsor-SES and βRC were highly correlated (Fig. 3c). The same results were obtained when βsor-SES was transformed by normal distribution.

Fig. 3. Relationship between dissimilarity/beta diversity values beyond null expectations as expressed in terms of standard effect size (βsor-SES) and the Raup-Crick metric (βRC). The βsor-SES values were transformed using logistic function. In all null models, the species richness of sites is maintained as in the observed. The species frequency was preserved in computing the βsor-SES (FF-null model: a,c), or equiprobably (FE-null model: b,d) (see method). The null models for computing the βRC sampled species proportional to observed frequencies (βRC: a,b) or equiprobably (\*βRC: c,d) (for details of the method see Chase et al., 2001).

Towards a Better Understanding of Beta Diversity: Deconstructing

\*\*\*p <0.001; \*\* <0.01; \* <0.05; **§** <0.10

Composition Patterns of Saproxylic Beetles Breeding in Recently Burnt Boreal Forest 87

Source of variability df *SS F- value SS F- value SS F- value*  Tree species (Ts) 1 0.911 **8.546\*\*\*** 0.706 **11.327\*\*\*** -0.030 -1.835 Burn severity (Bs) 2 0.717 **3.362\*\*\*** -0.058 -0.466 0.486 **14.956\*\*\***  Tree size/dbh class (Dc) 3 1.502 **4.696\*\*\*** 0.334 1.789 0.578 **11.869\*\*\***  Forest burns "zone" 3 0.484 *1.513***§** 0.394 *2.109* **§** 0.035 0.726 Ts x Bs 2 0.369 *1.729* **§** 0.281 *2.253* **§** -0.033 -1.026 Ts x Dc 3 0.317 0.992 0.181 0.969 0.082 1.683 Bs x Dc 6 0.555 0.868 0.240 0.642 0.213 *2.184* **§** Ts x Bs x Dc 6 0.684 1.069 0.300 0.804 0.050 0.510

Residuals 44 4.691 2.741 0.714 Total 70 10.230 5.119 2.094

indicated by bold typeface, and those indicated in italics are marginal.

Table 1. PERMANOVA table of saproxylic assemblages indicating the effect of habitat variables on the overall beta diversity (βsor) and its "spatial" turnover (βsim) and nestednessdriven (βnes) components. Note that βsor = βsim + βnes. Significance of the pseudo F-ratio was tested using a permutation test (9999 permutations); significant results Pr (>F-value) are

composition differentiation between burns were not significant (Tables 1; Fig. 5h).

The influence of tree species on the turnover component of beta diversity indicates that there is some level of differentiation in community ensembles between black spruce and jack pine. This pattern is expected among saproxylics that exhibit host-tree specificity (Allison *et al.*,2004; Janssen et al. in press) Such distinct habitat preferences or suitability of trees for component species can lead to segregated species distributions (Azeria *et al.*, 2010). It was intriguing that the within-treatment dissimilarity was higher for jack pine than black spruce;

Taken together our results suggest that the total beta diversity (βsor) pattern of saproxylic beetles was driven by differences between tree species and between tree size classes, as well as variation in within-treatment dissimilarity among burn severity classes and to some extent among sites within forest burns. The effects of these habitat attributes on overall beta diversity may be through influences on its "spatial" turnover (dissimilarity by species replacement) and/or nestedness (dissimilarity by richness variation) components. It is crucial that the two components of beta diversity are disentangled for a proper understanding of the most likely distinct underlying mechanisms (Baselga, 2008; 2010). Our results indicate that the turnover (βsim) and nestedness (βnes) components of beta diversity of saproxylic beetles are indeed dependent upon different habitat attributes (Tables 1 &2; Figs. 4 & 5). The turnover component of beta diversity was primarily driven by tree species, which showed a significant composition differentiation between black spruce and jack pine (PERMANOVA, Table 1; Fig. 4b). In addition, there was a significant difference in the within-treatment turnover among tree species: turnover was higher within jack pine than within black spruce (PERMDISP, Table 2; Fig. 4b; Fig. 5e). In fact, the withintreatment species turnover for jack pine was to the same extent as that found between jack pine and black spruce (Fig. 5e). The within forest-burn turnover was also lower for the two north burns (F1, F2) than for the southern (F3, F4) forest burns (Tables 2; Fig. 5h), but

Total (βsor ) Turnover (βsim) Nestedness (βnes)

These results underscore that a clear distinction should be made between the *"actual"* values of turnover beyond null model (βsim and βsor-diff) and the *probabilities* associated with these values as estimated from null model, using the βsor-SES and βRC metrics. As demonstrated above, a more direct link to the turnover component of beta diversity (βsim) is obtained by beta diversity beyond "null" distribution βsor-diff. It is noteworthy that βsor-diff will be negative when observed dissimilarity is less than expected (sites were more similar), and a proper rescaling should be made (e.g., subtract the minimum value) in subsequent analysis that require positive values (e.g. ordinations or PERMANOVA). Certainly, the null model derived βsor-SES (this study, also see Azeria *et al.*, 2009a) and βRC metrics are applicable and important measures for studying beta diversity independent of richness variation (Chase, 2007; Chase *et al.*, 2011). However, caution should be exercised as the values can become biased for site pairs that have extremely low species richness relative to the regional species pool. For example, although the difference between the observed and null expectation βsor-diff is small, the variation of the null expectations βsor-sd might be so small that βsor-SES become inflated (for related caveats with βRC and other issue see Chase *et al*., 2011). Although, the inverselogit/normal transformation should minimize the bias (see Azeria *et al*, 2009a), caution should still be exercised when using the value in subsequent analysis. In addition, our results offer interesting qualitative comparisons among the null model based indices of beta diversity and how constraints imposed on species incidence would influence the values.

#### **4.2 Effect of habitat factors for overall and components of beta diversity**

We found a significant effect of tree species, burn severity, and tree-size class on overall beta diversity (βsor) of saproxylic beetles (PERMANOVA Table 1, Fig. 4a, d, and g). In addition, we found a marginal effect of interaction terms between tree species and burn severity. The effect of tree species and tree-size class was primarily due to compositional difference between treatments (location of treatment in multivariate response) but not due to differences in the within-class dispersion (PERMDISP, Table 2). In contrast, the influence of burn severity was primarily due to differences in the within-class dispersion, which was lower for low-severity burn (more homogeneous) than that of moderate- and high-severity burns (PERMDISP, Table 2). Results from multivariate regression analysis on distance matrices (MRM) provide a concise summary of the simultaneous effect of within- and between-treatment on the overall beta diversity pattern (Fig. 5a, b, c and d). The MRM also showed some trends that were not evident through the use of PERMANOVA and PERMIDISP. For example, the beta diversity or composition dissimilarity within jack pine (JPI-JPI) was similar to that found between jack pine and black spruce (BSP-JPI) (Fig 5a).

Overall, the effect of geographical distance on beta diversity patterns was only marginal, and when detected it was due to differences in the within-burns dissimilarity (lower for F1, forest burn in the north, than the others, Fig 5d) rather than between-burns differences. In other words, there was no increase in composition dissimilarity of saproxylic beetles in burned forests with increasing site-to-site geographical distance (βsor: r= 0.032; βsim: r= 0.040; βnes: r= -0.017). Thus our results do not provide support for the "distance decay of similarity" hypothesis (Nekola & White, 1999). It seems that saproxylic beetles might be good dispersers due to the ephemeral nature of their habitats (Boulanger *et al.*, 2010) and thus may not be strongly limited by dispersal, at least at the scale of our study (up to 200 km). Baselga (2010) has shown that composition dissimilarity of longhorn beetles increase with geographical distance measured across larger scales (up to 3000 km) across Europe.


Towards a Better Understanding of Beta Diversity: Deconstructing Composition Patterns of Saproxylic Beetles Breeding in Recently Burnt Boreal Forest 87

\*\*\*p <0.001; \*\* <0.01; \* <0.05; **§** <0.10

86 Research in Biodiversity – Models and Applications

These results underscore that a clear distinction should be made between the *"actual"* values of turnover beyond null model (βsim and βsor-diff) and the *probabilities* associated with these values as estimated from null model, using the βsor-SES and βRC metrics. As demonstrated above, a more direct link to the turnover component of beta diversity (βsim) is obtained by beta diversity beyond "null" distribution βsor-diff. It is noteworthy that βsor-diff will be negative when observed dissimilarity is less than expected (sites were more similar), and a proper rescaling should be made (e.g., subtract the minimum value) in subsequent analysis that require positive values (e.g. ordinations or PERMANOVA). Certainly, the null model derived βsor-SES (this study, also see Azeria *et al.*, 2009a) and βRC metrics are applicable and important measures for studying beta diversity independent of richness variation (Chase, 2007; Chase *et al.*, 2011). However, caution should be exercised as the values can become biased for site pairs that have extremely low species richness relative to the regional species pool. For example, although the difference between the observed and null expectation βsor-diff is small, the variation of the null expectations βsor-sd might be so small that βsor-SES become inflated (for related caveats with βRC and other issue see Chase *et al*., 2011). Although, the inverselogit/normal transformation should minimize the bias (see Azeria *et al*, 2009a), caution should still be exercised when using the value in subsequent analysis. In addition, our results offer interesting qualitative comparisons among the null model based indices of beta diversity and

how constraints imposed on species incidence would influence the values.

**4.2 Effect of habitat factors for overall and components of beta diversity** 

We found a significant effect of tree species, burn severity, and tree-size class on overall beta diversity (βsor) of saproxylic beetles (PERMANOVA Table 1, Fig. 4a, d, and g). In addition, we found a marginal effect of interaction terms between tree species and burn severity. The effect of tree species and tree-size class was primarily due to compositional difference between treatments (location of treatment in multivariate response) but not due to differences in the within-class dispersion (PERMDISP, Table 2). In contrast, the influence of burn severity was primarily due to differences in the within-class dispersion, which was lower for low-severity burn (more homogeneous) than that of moderate- and high-severity burns (PERMDISP, Table 2). Results from multivariate regression analysis on distance matrices (MRM) provide a concise summary of the simultaneous effect of within- and between-treatment on the overall beta diversity pattern (Fig. 5a, b, c and d). The MRM also showed some trends that were not evident through the use of PERMANOVA and PERMIDISP. For example, the beta diversity or composition dissimilarity within jack pine (JPI-JPI) was similar to that found between jack pine and black spruce (BSP-JPI) (Fig 5a). Overall, the effect of geographical distance on beta diversity patterns was only marginal, and when detected it was due to differences in the within-burns dissimilarity (lower for F1, forest burn in the north, than the others, Fig 5d) rather than between-burns differences. In other words, there was no increase in composition dissimilarity of saproxylic beetles in burned forests with increasing site-to-site geographical distance (βsor: r= 0.032; βsim: r= 0.040; βnes: r= -0.017). Thus our results do not provide support for the "distance decay of similarity" hypothesis (Nekola & White, 1999). It seems that saproxylic beetles might be good dispersers due to the ephemeral nature of their habitats (Boulanger *et al.*, 2010) and thus may not be strongly limited by dispersal, at least at the scale of our study (up to 200 km). Baselga (2010) has shown that composition dissimilarity of longhorn beetles increase with geographical distance measured across larger scales (up to 3000 km) across Europe.

Table 1. PERMANOVA table of saproxylic assemblages indicating the effect of habitat variables on the overall beta diversity (βsor) and its "spatial" turnover (βsim) and nestednessdriven (βnes) components. Note that βsor = βsim + βnes. Significance of the pseudo F-ratio was tested using a permutation test (9999 permutations); significant results Pr (>F-value) are indicated by bold typeface, and those indicated in italics are marginal.

Taken together our results suggest that the total beta diversity (βsor) pattern of saproxylic beetles was driven by differences between tree species and between tree size classes, as well as variation in within-treatment dissimilarity among burn severity classes and to some extent among sites within forest burns. The effects of these habitat attributes on overall beta diversity may be through influences on its "spatial" turnover (dissimilarity by species replacement) and/or nestedness (dissimilarity by richness variation) components. It is crucial that the two components of beta diversity are disentangled for a proper understanding of the most likely distinct underlying mechanisms (Baselga, 2008; 2010).

Our results indicate that the turnover (βsim) and nestedness (βnes) components of beta diversity of saproxylic beetles are indeed dependent upon different habitat attributes (Tables 1 &2; Figs. 4 & 5). The turnover component of beta diversity was primarily driven by tree species, which showed a significant composition differentiation between black spruce and jack pine (PERMANOVA, Table 1; Fig. 4b). In addition, there was a significant difference in the within-treatment turnover among tree species: turnover was higher within jack pine than within black spruce (PERMDISP, Table 2; Fig. 4b; Fig. 5e). In fact, the withintreatment species turnover for jack pine was to the same extent as that found between jack pine and black spruce (Fig. 5e). The within forest-burn turnover was also lower for the two north burns (F1, F2) than for the southern (F3, F4) forest burns (Tables 2; Fig. 5h), but composition differentiation between burns were not significant (Tables 1; Fig. 5h).

The influence of tree species on the turnover component of beta diversity indicates that there is some level of differentiation in community ensembles between black spruce and jack pine. This pattern is expected among saproxylics that exhibit host-tree specificity (Allison *et al.*,2004; Janssen et al. in press) Such distinct habitat preferences or suitability of trees for component species can lead to segregated species distributions (Azeria *et al.*, 2010). It was intriguing that the within-treatment dissimilarity was higher for jack pine than black spruce;

Towards a Better Understanding of Beta Diversity: Deconstructing

Composition Patterns of Saproxylic Beetles Breeding in Recently Burnt Boreal Forest 89

Fig. 4. A two-dimensional MDS representation of site-to-site dissimilarity showing the effect of tree species (a-c), burn severity (d-f) and tree-size class (g-i) on saproxylic beetle's total beta diversity(a,d,g), and when partitioned into "spatial" turnover (b,e,h) and nestednessdriven (c,f,i) components. The ellipse encloses 80% of the dispersion of habitat classes from their respective group centroid as accounted by the first two axes. Habitat classes are Tree species: BSP=Black spruce and JPI=Jack pine; Burn severity: L=Low, M=Moderate and H=High; Tree size/dbh (diameter at breast height in cm) classes: db1=8-12, db2=12-16,

The nestedness-driven beta diversity (βnes), however, was influenced by gradients of burn severity and tree size/dbh, but not by tree species (PERMANOVA, Table 1, Fig 4f, i). There was significant nestedness-driven composition dissimilarity between high severity burns and both low- and moderate-severity burns (Fig 4f and Fig. 5j); and between small-sized trees (Db1) and large-sized trees (Db3 and Db4) (Fig 4i; Fig 5k). In addition, there was a significant difference in the within-class dissimilarity for burn severity; βnes was higher for the high-severity burns and moderate-severity burns than in the low-severity burn class (PERMDISP; Table 2, Fig 4f and Fig 5j). Finally, although PERMIDISP indicated a nonsignificant result for tree species, the MRM analysis indicated that the nestedness-driven dissimilarity was higher within black spruce than jack pine (Fig. 5i), which was opposite to

db3=16-20 and db4=20-24.

that observed for turnover component (Fig. 5e).

this might be related to jack pine being less common (although widely distributed) than black spruce in the landscape. The high turnover within the southern forest burns compared to the northern forest burns may be related to their higher heterogeneity in terms of composition and structure. Notably, the dissimilarity within and between the southern burns (F3 and F4) was of the same extent as that observed with respect to the northern forest burns (F1 and F2).


\*\* P<0.01; \* <0.05; **§** <0.10

Table 2. Summary of PERMDISP analysis examining differences in the within-treatment dispersion/dissimilarity based on metrics of overall, turnover and nestedness components of beta diversity. Significant differences are indicated in bold letters; while bold-italics indicate marginally significant results (9999 permutations). Treatments or classes of habitat attributes considered are- Tree species: BSP=Black spruce and JPI=Jack pine; Burn severity: L=Low, M=Moderate and H=High; Tree size/dbh (diameter at breast height in cm) classes: Db1=8-12, Db2=12-16, Db3=16-20 and Db4=20-24. Forest burns were located north to south: F1-F2-F3/F4.

this might be related to jack pine being less common (although widely distributed) than black spruce in the landscape. The high turnover within the southern forest burns compared to the northern forest burns may be related to their higher heterogeneity in terms of composition and structure. Notably, the dissimilarity within and between the southern burns (F3 and F4) was of the same extent as that observed with respect to the northern forest

Groups 1 0.029 2.697 0.065 **5.802**\* 0.014 1.594

Groups 2 0.085 **4.583**\* 0.033 1.496 0.077 **5.991\*\*** 

Residual 69 0.751 0.770 0.587 Pairwise tests BSP=JPI BSP <JPI BSP=JPI

Residual 68 0.634 0.755 0.438

Residual 67 0.694 0.712 0.443

Residual 67 0.579 0.674 0.496

Table 2. Summary of PERMDISP analysis examining differences in the within-treatment dispersion/dissimilarity based on metrics of overall, turnover and nestedness components of beta diversity. Significant differences are indicated in bold letters; while bold-italics indicate marginally significant results (9999 permutations). Treatments or classes of habitat attributes considered are- Tree species: BSP=Black spruce and JPI=Jack pine; Burn severity: L=Low, M=Moderate and H=High; Tree size/dbh (diameter at breast height in cm) classes: Db1=8-12, Db2=12-16, Db3=16-20 and Db4=20-24. Forest burns were located north to south:

Db1=Db2=Db3=

Db4

tests L < (M=H) L =M=H L < (M**§**=H)

Groups 3 0.068 *2.185***§** 0.047 1.486 0.043 *2.168§*

Groups 3 0.111 **4.279**\*\* 0.112 **3.702**\* 0.053 *2.377§*

Pairwise tests F1<(F3=F4)=F2 (F1=F2) *§* < F3=F4 F1=F2=F3=F4

Total (βsor ) Turnover (βsim) Nestedness (βnes) *SS F value SS F value SS F value*

Db1=Db2=Db3=

Db4 Db1=Db2=Db3=Db4

burns (F1 and F2).

Tree species (Ts)

Burn severity (Bs)

Tree size/dbh class (Dc)

Forest burns "zone" 

\*\* P<0.01; \* <0.05; **§** <0.10

F1-F2-F3/F4.

Source of variability df

Pairwise

Pairwise tests

Fig. 4. A two-dimensional MDS representation of site-to-site dissimilarity showing the effect of tree species (a-c), burn severity (d-f) and tree-size class (g-i) on saproxylic beetle's total beta diversity(a,d,g), and when partitioned into "spatial" turnover (b,e,h) and nestednessdriven (c,f,i) components. The ellipse encloses 80% of the dispersion of habitat classes from their respective group centroid as accounted by the first two axes. Habitat classes are Tree species: BSP=Black spruce and JPI=Jack pine; Burn severity: L=Low, M=Moderate and H=High; Tree size/dbh (diameter at breast height in cm) classes: db1=8-12, db2=12-16, db3=16-20 and db4=20-24.

The nestedness-driven beta diversity (βnes), however, was influenced by gradients of burn severity and tree size/dbh, but not by tree species (PERMANOVA, Table 1, Fig 4f, i). There was significant nestedness-driven composition dissimilarity between high severity burns and both low- and moderate-severity burns (Fig 4f and Fig. 5j); and between small-sized trees (Db1) and large-sized trees (Db3 and Db4) (Fig 4i; Fig 5k). In addition, there was a significant difference in the within-class dissimilarity for burn severity; βnes was higher for the high-severity burns and moderate-severity burns than in the low-severity burn class (PERMDISP; Table 2, Fig 4f and Fig 5j). Finally, although PERMIDISP indicated a nonsignificant result for tree species, the MRM analysis indicated that the nestedness-driven dissimilarity was higher within black spruce than jack pine (Fig. 5i), which was opposite to that observed for turnover component (Fig. 5e).

Towards a Better Understanding of Beta Diversity: Deconstructing

**5. Conclusion** 

Composition Patterns of Saproxylic Beetles Breeding in Recently Burnt Boreal Forest 91

variance (relative to average) in species richness exhibited within severe burns. We suspect that an interaction effect between burn severity and tree-size might play a role, but this is

The above results clearly underscore that the turnover and nestedness components of saproxylic beetles are driven by different post-fire habitat legacies. Disentangling causal factors influencing beta diversity patterns of saproxylic beetles can have important conservation implications, such as in post-fire salvage logging operations where maintaining saproxylic diversity can concern management plans. Indeed, concern over ecological consequence of post-fire salvage logging on biodiversity and ecosystem function is a pressing management issue (Lindenmayer *et al.*, 2008), and saproxylic beetles constitute key component of burned forest ecosystems (Grove, 2002; Cobb *et al.*, 2010). Knowledge about the factors influencing beta divesity components can be crucial in setting salvage logging practices that will conserve also these essential groups. For example, given the species composition turnover exhibited between jack pine and black spruce, conserving the totality of saproxylic beetle species will require a management approach that maintains a mosaic of both tree species in the landscape. In addition, higher turnover within-jack pine than black spruce (and within moderate severity burns) may require special considerations to conserve all species occupying that habitat. On the other hand, the nestedness driven component will perhaps require emphasis on habitat attributes that increase species richness, e.g., large trees of lower severity burns. These decisions could also have implications for the structure of forest stands that are to be set aside for protective purposes. For example, in old-growth stands with numerous large-diameter trees (e.g., Db4=20-24), a smaller number of trees might suffice to capture the totality of species given the low differentiation of both turnover and nestedness component of beta diversity within large trees (given all other attributes are considered). However, if the available forest stands contain only moderately sized trees (e.g., Db2=12-16), then protecting more trees during salvage-logging might be required (given high turnover), although variation of nestednessdriven beta diversity was to the same extent as that of larger trees. These examples are just to illustrate the implication of disentangling the turnover and nestedness components and underlying factors for management (also see Azeria *et al.*, 2006; 2009b; Baselga, 2010).

Ecologists face a continuing challenge to disentangle and explain the species 'turnover' from composition dissimilarity that is driven by variation in species richness among sites. The recently proposed beta diversity partitioning framework and null model approaches make a significant step forward in meeting this challenge. We demonstrate the explicit quantitative relationship of the beta diversity components as depicted in the partitioning framework with that expected *under* and *beyond* random assembly of communities by using null models. Our results indicate that the turnover component indeed reflects the beta diversity deviating *beyond* that expected from "random" assembly given species richness variations, while the nestedness component of beta diversity was related to the *null* expectations. In addition, the beta diversity components were conditioned by variation in species frequency; this also implies that null models for beta diversity studies should perhaps preserve both the observed richness and species frequency patterns. We concur with others (Baselga 2010, Anderson et al. 2011) that the distinction between the two beta diversity components is important, because they were driven by different underlying causes (habitat factors) and this has implications for post-fire management in the boreal forest. Additional studies that directly incorporate habitat and

difficult to isolate statistically from their independent effects in PERMDISP analysis.

Fig. 5. A simultaneous analysis of beta diversity patterns (βsor, βsim, βnes) differences within- (gray bars) and between-habitat (black bars) classes according tree species (a,e,i), burn severity (b,f,j), tree-size class (c, g, k) and burns (d, h, i) on saproxylic beetle's total beta diversity(a-d), and when partitioned into "spatial" turnover (e-h) and nestedness-driven (i-l) components. On the y-axis are mean values of dissimilarity (beta diversity) for the corresponding habitat contrasts. Habitat codes are as in Fig. 4.

The effects of burn severity and tree size gradients on the nestedness component of beta diversity indicate that hierarchical suitability of the attributes is causing richness-driven composition differences between respective treatments. Generally, large trees of low burn severity are more suitable to saproxylic beetles than small trees of high burn severity, because they provide suitable oviposition conditions and high quality food resources (i.e., subcortical tissues) for larvae feeding directly on the bark and wood of trees (Allison *et al.*, 2004; Saint-Germain *et al.*, 2004). Indirectly, these trees may also be more suitable to predators of these larvae (Kenis *et al.*, 2004). Indeed, large-sized and/or less severly burned treess have been shown to support more saproxylic beetle species than small-sized and/or very severely burned trees (Azeria *et al.*, 2010; Boulanger *et al.*, 2010). It is interesting, and counter intuitive, that a similar trend was also observed in the within-class nestedness-driven composition dissimilarity among burn severity classes; thus, nestedness-driven dissimilarity was higher for high-severity than low-severity burns. This might be explained by the high

variance (relative to average) in species richness exhibited within severe burns. We suspect that an interaction effect between burn severity and tree-size might play a role, but this is difficult to isolate statistically from their independent effects in PERMDISP analysis.

The above results clearly underscore that the turnover and nestedness components of saproxylic beetles are driven by different post-fire habitat legacies. Disentangling causal factors influencing beta diversity patterns of saproxylic beetles can have important conservation implications, such as in post-fire salvage logging operations where maintaining saproxylic diversity can concern management plans. Indeed, concern over ecological consequence of post-fire salvage logging on biodiversity and ecosystem function is a pressing management issue (Lindenmayer *et al.*, 2008), and saproxylic beetles constitute key component of burned forest ecosystems (Grove, 2002; Cobb *et al.*, 2010). Knowledge about the factors influencing beta divesity components can be crucial in setting salvage logging practices that will conserve also these essential groups. For example, given the species composition turnover exhibited between jack pine and black spruce, conserving the totality of saproxylic beetle species will require a management approach that maintains a mosaic of both tree species in the landscape. In addition, higher turnover within-jack pine than black spruce (and within moderate severity burns) may require special considerations to conserve all species occupying that habitat. On the other hand, the nestedness driven component will perhaps require emphasis on habitat attributes that increase species richness, e.g., large trees of lower severity burns. These decisions could also have implications for the structure of forest stands that are to be set aside for protective purposes. For example, in old-growth stands with numerous large-diameter trees (e.g., Db4=20-24), a smaller number of trees might suffice to capture the totality of species given the low differentiation of both turnover and nestedness component of beta diversity within large trees (given all other attributes are considered). However, if the available forest stands contain only moderately sized trees (e.g., Db2=12-16), then protecting more trees during salvage-logging might be required (given high turnover), although variation of nestednessdriven beta diversity was to the same extent as that of larger trees. These examples are just to illustrate the implication of disentangling the turnover and nestedness components and underlying factors for management (also see Azeria *et al.*, 2006; 2009b; Baselga, 2010).

## **5. Conclusion**

90 Research in Biodiversity – Models and Applications

Fig. 5. A simultaneous analysis of beta diversity patterns (βsor, βsim, βnes) differences within- (gray bars) and between-habitat (black bars) classes according tree species (a,e,i), burn severity (b,f,j), tree-size class (c, g, k) and burns (d, h, i) on saproxylic beetle's total beta diversity(a-d), and when partitioned into "spatial" turnover (e-h) and nestedness-driven (i-l)

The effects of burn severity and tree size gradients on the nestedness component of beta diversity indicate that hierarchical suitability of the attributes is causing richness-driven composition differences between respective treatments. Generally, large trees of low burn severity are more suitable to saproxylic beetles than small trees of high burn severity, because they provide suitable oviposition conditions and high quality food resources (i.e., subcortical tissues) for larvae feeding directly on the bark and wood of trees (Allison *et al.*, 2004; Saint-Germain *et al.*, 2004). Indirectly, these trees may also be more suitable to predators of these larvae (Kenis *et al.*, 2004). Indeed, large-sized and/or less severly burned treess have been shown to support more saproxylic beetle species than small-sized and/or very severely burned trees (Azeria *et al.*, 2010; Boulanger *et al.*, 2010). It is interesting, and counter intuitive, that a similar trend was also observed in the within-class nestedness-driven composition dissimilarity among burn severity classes; thus, nestedness-driven dissimilarity was higher for high-severity than low-severity burns. This might be explained by the high

components. On the y-axis are mean values of dissimilarity (beta diversity) for the

corresponding habitat contrasts. Habitat codes are as in Fig. 4.

Ecologists face a continuing challenge to disentangle and explain the species 'turnover' from composition dissimilarity that is driven by variation in species richness among sites. The recently proposed beta diversity partitioning framework and null model approaches make a significant step forward in meeting this challenge. We demonstrate the explicit quantitative relationship of the beta diversity components as depicted in the partitioning framework with that expected *under* and *beyond* random assembly of communities by using null models. Our results indicate that the turnover component indeed reflects the beta diversity deviating *beyond* that expected from "random" assembly given species richness variations, while the nestedness component of beta diversity was related to the *null* expectations. In addition, the beta diversity components were conditioned by variation in species frequency; this also implies that null models for beta diversity studies should perhaps preserve both the observed richness and species frequency patterns. We concur with others (Baselga 2010, Anderson et al. 2011) that the distinction between the two beta diversity components is important, because they were driven by different underlying causes (habitat factors) and this has implications for post-fire management in the boreal forest. Additional studies that directly incorporate habitat and

Towards a Better Understanding of Beta Diversity: Deconstructing

Chicoutimi, Canada.

*Management*, 260, 1114-1123

diversity. *Ecosphere*, 2, 19-28

10, 337-343

Washington, DC

University, Princiton, NJ.

Composition Patterns of Saproxylic Beetles Breeding in Recently Burnt Boreal Forest 93

Bergeron, Y., Gauthier, S., Flannigan, M. & Kafka, V. (2004) Fire regimes at the transition between

Boucher, J. (2011) *Impacts de la coupe de récupération après feu sur les coléoptéres associés aux* 

Boulanger, Y., Sirois, L. & Hébert, C. (2010) Distribution of saproxylic beetles in a recently

Chase, J. (2007) Drought mediates the importance of stochastic community assembly.

Chase, J. M. (2010) Stochastic community assembly causes higher biodiversity in more

Chase, J. M., Kraft, N. J. B., Smith, K. G., Vellend, M. & Inouye, B. D. (2011) Using null

Cobb, T. P., Hannam, K. D., Kishchuk, B. E., Langor, D. W., Quideau, S. A. & Spence, J. R.

Connor, E. F. & Simberloff, D. (1983) Interspecific competition and species co-occurrence patterns on islands: null models and the evaluation of evidence. *Oikos*, 41, 455-465 Gotelli, N. J. (2001) Research frontiers in null model analysis. *Global Ecology and Biogeography*,

Gotelli, N. J. & Graves, G. R. (1996) *Null Models in Ecology*, edn. Smithsonian Institution,

Grove, S. J. (2002) Saproxylic insect ecology and the sustainable management of forests.

Harrison, S., Ross, S. J. & Lawton, J. H. (1992) Beta diversity on geographic gradients in

Hayden, D., Lazar, P. & Schoenfeld, D. (2009) Assessing statistical significance in microarray

Hubbell, S. P. (2001) *The unified neutral theory of biodiversity and biogeography*, edn. Princeton

Huston, M. A. (1994) *Biological diversity: the coexistence of species on changing landscapes*.

Hylander, K., Nilsson, C., Gunnar Jonsson, B. & Göthner, T. (2005) Differences in habitat quality explain nestedness in a land snail meta community. *Oikos*, 108, 351-361. Janssen, P., Hébert, C. & Fortin, D. (in press). Biodiversity conservation in old-growth boreal

Kenis, M., Wermelinger, B. & Grégoire, J. (2004) Research on parasitoids and predators of

Koleff, P., Gaston, K. & Lennon, J. (2003) Measuring beta diversity for presence–absence

Lande, R. (1996) Statistics and partitioning of species diversity, and similarity among

Legendre, P., Borcard, D. & Peres-Neto, P. R. (2008) Analyzing or explaining beta diversity?

forest: black spruce and balsam fir snags harbour distinct assemblages of saproxylic

Scolytidae–a review. *Bark and Wood Boring Insects in Living Trees in Europe, a Synthesis* (ed. by F. Lieutier & K.R. Day & A. Battisti & J.-C. Grégoire & H.F. Evans),

experiments using the distance between microarrays. *PLoS One*, 4.

*Proceedings of the National Academy of Sciences*, 104, 17430-17434

productive environments. *science*, 328, 1388-1391

*Annual Review of Ecology and Systematics*, 33, 1-23

Britain. *Journal of Animal Ecology*, 61, 151-158

Cambridge University Press, Cambridge, UK

pp 237-290. Kluwer Academic Publishers, Dordrecht

beetles. Biodiversity and Conservation.

data. *Journal of Animal Ecology*, 72, 367-382

multiple communities. *Oikos*, 76, 5-13

Comment. *Ecology*, 89, 3238-3244

mixedwood and coniferous boreal forest in northwestern Quebec. *Ecology*, 85, 1916-1932

*brûlis en forêt boréale: une dynamque temporaelle*. Université du Québec à Chicoutim,

burnt landscape of the northern boreal forest of Québec. *Forest Ecology and* 

models to disentangle variation in community dissimilarity from variation in α-

(2010) Wood feeding beetles and soil nutrient cycling in burned forests: implications of post fire salvage logging. *Agricultural and Forest Entomology*, 12, 9-18

spatial effects into null model analysis such as using habitat- and spatially-constrained null models are needed in order to increase our understanding of the factors that control beta diversity patterns across space (also see Chase *et al*., 2011).

#### **6. Acknowledgements**

We thank Chantiers Chibougamau and Barrette Chapais Ltée for logistical support, and all of our field assistants that helped in conducting the field work and in sample sorting. We are indebted to Y. Dubuc and G. Pelletier for specimen identification. We thank J. Hodson for his comments and revising the English. This study was supported by the *Fonds de recherche sur la nature et les technologies* and the *Ministère des Ressources naturelles et de la Faune*  through the *Programme de recherche en partenariat sur l'aménagement et l'environnement forestiers*-II, by the iFor consortium (Université Laval), by Environment Canada, by the Canadian Forest Service, and by *la Fondation de l'Université du Québec à Chicoutimi*.

#### **7. References**


spatial effects into null model analysis such as using habitat- and spatially-constrained null models are needed in order to increase our understanding of the factors that control beta

We thank Chantiers Chibougamau and Barrette Chapais Ltée for logistical support, and all of our field assistants that helped in conducting the field work and in sample sorting. We are indebted to Y. Dubuc and G. Pelletier for specimen identification. We thank J. Hodson for his comments and revising the English. This study was supported by the *Fonds de recherche sur la nature et les technologies* and the *Ministère des Ressources naturelles et de la Faune*  through the *Programme de recherche en partenariat sur l'aménagement et l'environnement forestiers*-II, by the iFor consortium (Université Laval), by Environment Canada, by the

Allison, J., Borden, J. & Seybold, S. (2004) A review of the chemical ecology of the

Anderson, M., Crist, T., Chase, J., Vellend, M., Inouye, B., Freestone, A., Sanders, N., Cornell,

Anderson, M. J. (2001) A new method for non parametric multivariate analysis of variance.

Anderson, M. J., Ellingsen, K. E. & Mcardle, B. H. (2006) Multivariate dispersion as a

Atmar, W. & Patterson, B. D. (1993) The measure of order and disorder in the distribution of

Azeria, E. T., Bouchard, M., Pothier, D., Fortin, D. & Hébert, C. (2011) Using biodiversity

Azeria, E. T., Carlson, A., Pärt, T. & Wiklund, C. G. (2006) Temporal dynamics and nestedness of an oceanic island bird fauna. *Global Ecology and Biogeography*, 15, 328-338 Azeria, E. T., Fortin, D., Hébert, C., Peres-Neto, P., Pothier, D. & Ruel, J. C. (2009a) Using

Azeria, E. T., Fortin, D., Lemaître, J., Janssen, P., Hébert, C., Darveau, M. & Cumming, S. G.

Azeria, E. T. & Kolasa, J. (2008) Nestedness, niche metrics and temporal dynamics of a metacommunity in a dynamic natural model system. *Oikos*, 117, 1006-1019 Baselga, A. (2008) Determinants of species richness, endemism and turnover in European

Baselga, A. (2010) Partitioning the turnover and nestedness components of beta diversity.

and select indicator species. *Diversity and Distributions*, 15, 958-971

H., Comita, L. & Davies, K. (2011) Navigating the multiple meanings of diversity: a

deconstruction to disentangle assembly and diversity dynamics of understorey plants along post fire succession in boreal forest. *Global Ecology and Biogeography*, 20, 119-133

null model analysis of species co-occurrences to deconstruct biodiversity patterns

(2009b) Fine scale structure and cross taxon congruence of bird and beetle assemblages in an old growth boreal forest mosaic. *Global Ecology and Biogeography*, 18, 333-345 Azeria, E. T., Ibarzabal, J. & Hebért, C. (2010) Importance of habitat and biotic interactions

for assembly of Saproxylic beetles emerging from logs: A test using habitatconstrained null models. *Ecological Society of America (ESA) 95th Annual Meeting*,

Canadian Forest Service, and by *la Fondation de l'Université du Québec à Chicoutimi*.

Cerambycidae (Coleoptera). *Chemoecology*, 14, 123-150

measure of beta diversity. *Ecology Letters*, 9, 683-693

species in fragmented habitat. *Oecologia*, 96, 373-382

roadmap for the practicing ecologist. *Ecology Letters*, 14, 19-28

diversity patterns across space (also see Chase *et al*., 2011).

**6. Acknowledgements** 

**7. References** 

*Austral Ecology*, 26, 32-46

Pittsburgh, PA

longhorn beetles. *Ecography*, 31, 263-271.

*Global Ecology and Biogeography*, 19, 134-143


**Part 2** 

**History of Biodiversity** 

