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

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,

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 c o m p o n e n t s o f b e t a d i v e r s i t y h a v e d i f f e r e n t 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 k n o w n . I n t h i s s t u d y , w e e x a m i n e t h e i r e f f e ct 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).

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 (49 o 15'-50 o 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).

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.

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).

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 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 = sorsor-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-SESsor-SES.min )/ ( sor-SES.maxsor-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 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.

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 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.

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 www.intechopen.com 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 . 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). 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 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).
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.

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;. 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; 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).

Source of variability df
Total (   ) 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). 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).

Conclusion
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 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).

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.