Data & codes - From forest stand decline to salvage logging: cascading impacts on saproxylic beetle diversity

1. Sampling design Our sampling design is the result of the opportunistic combination of three different forest contexts: two in France (the Silver fir-dominated forests Abies alba Mill.) in French Pyrenees, hereafter in text "fir-dominated forests", and the oak-dominated forests (mix of Quercus petraea and Quercus robur) in Loire valley, and one in Germany (the Norway spruce-dominated forests (Picea abies) in Bavaria, hereafter in text "spruce-dominated forests" (Figure 1, Table 1). They originated from three different projects: CLIMTREE for the 2 montane coniferous contexts (spruce context were also funded via BIOKLIM project; Bässler et al., 2009; Thorn et al., 2017), while the oak-dominated forest context was funded via CANOPEE project (Sallé et al., 2020). While they all share the same methodology to describe the forest structure—i.e. deadwood and TreMs—, they show differences in sampling design and beetle trapping methods (see Table 1). The three projects shared the common objective of studying the consequences of forest decline on the biodiversity (with an important focus on saproxylic beetles). Along a gradient of stand decline, we selected 56 study plots in fir-dominated forests (28 in Aure valley (Central Pyrenees; 854 to 1570 m a.s.l.; 42°51'46.8"N, 0°36'08.9"E) and 28 on Sault plateau (Eastern Pyrenees; 705 to 1557.3 m a.s.l.; 42°50'58.7"N, 2°00'41.3"E)), 15 plots in oak-dominated forests (6 in Orleans forest (107 to 174 m a.s.l.; 48°00'28"N, 2°09'39"E) and 9 in Vierzon forest (120 to 190 m a.s.l.; 47°15'53.9"N, 2°05'55.4"E)), and 28 plots in spruce-dominated forests (660 to 1352 m a.s.l.; 49°1'24.96"N 13°38'5.222"E; Figure 1). Fir and spruce forests included 12 and 9 plots respectively which have been salvage logged (Figure 1, Table 1), while the most valuable oak trees were regularly harvested during the decline process, which reduced the volume of deadwood in this context. In fir and oak forests, plots were set up in managed forests and the forests around were also predominantly managed. Plots in spruce forests were set up in the Bavarian Forest National Park, which is characterized by large areas of unmanaged disturbance areas but also by intervention sites with salvage logging in the buffer and former development zones of the National Park (see Müller et al., 2010; Thorn et al., 2016). Stand decline in fir and oak forests were mainly the result of drought events (the 2003 drought being a major event in the occurrence of these declines), which mainly resulted in patch dynamics (Sallé et al., 2020). In spruce forests, the decline studied was the result of windstorms and bark beetle outbreaks, which mainly resulted in stand-replacing dynamics (Thorn et al., 2017). The most recent impactful event for the spruce forests was the effect of the windstorm Kyrill in 2007, from which most of the decline in our studied plots was the results, compounded by the subsequent bark beetle outbreaks. 2. Environmental conditions To assess the dendrometric characteristic of each plot, we implemented a sampling protocol using relascope with 1/50 ratio (factor 1). We recorded the type (i.e. living tree, standing dead tree and lying deadwood), the tree species, the diameter at breast height (dbh) on living trees and standing dead trees higher than 4 m, and at mid-height on standing trees lower than 4 m and on lying deadwood), and the decay stage for deadwood (from 1 = fresh deadwood with full cover of bark, to 4 = soft wood with no more remaining barks). In addition, for each recorded tree, we noted the tree-related microhabitats (TreMs) on living and standing dead trees, based on the 47 types described by Larrieu et al. (2018). Mean plot area was approximately of 0.3 ha and depended on the dbh of the largest trees and their location from the center. Finally, we converted plot data into 1 ha densities by allocating the coefficient Nd to every measured tree and deadwood piece (Pardé & Bouchon, 1988) We calculated the coefficient Nd as Nd = 𝜋 × 108×((arctan(1/50))/(𝜋 × dbh))2. Such allocation is important for the better representation of tree diameter classes (small trees being under evaluated with relascopic sampling protocol [and vice versa with larger trees]), and result comparability (standardization at the scale of 1 ha). Field data were collected in 2017 in the fir and spruce forests, and in 2020 in oak forests. In parallel, we used eco-morphological traits related to each ligneous substrates (i.e. living trees and deadwood) and TreM, to describe their changes. For ligneous substrates, we directly used field measurements i.e. type of substrate (as an ordinal scale: 1 = "lying deadwood", 2 = "standing deadwood", 3 = "living trees"), decay stage (from 1 = "hard deadwood fully covered by bark”, to 4 = "very soft wood out of bark"), and diameter classes (see Bouget et al., 2024). In addition, we used a list of two eco-morphological traits characterizing each TreM found in the field i.e. (i) type of substrate bearing the TreM (as an ordinal scale: 1 = "lying deadwood", 2 = "standing deadwood", 3 = "living trees") and (ii) its association with deadwood (whether the TreM contains decaying deadwood, mould, or neither, hereafter named “saproxylic TreMs”; for additional information, see Bouget et al., 2024; data updated). We used the quantitative eco-morphological traits from ligneous substrates and TreMs to compute community weighted means (CWM), functional dispersion (FDis) for each of them, and overall functional dispersion, using FD R-package (Laliberté et al., 2014; for more information, see Bouget et al., 2024). We used the sum of weight Nd as abundance measure. We also calculated the volume of the different types of deadwood, i.e. standing and lying, low and high decay level (i.e. level ≥ 3 on a scale from 1 to 4), small and large (dbh ≥ 40 cm), and deadwood diversity combining all these features with tree species (Siitonen et al., 2000). Volumes were calculated using length and diameter for all logs and standing deadwood of less than 4 m tall, and we used Schaeffer’s volume equations for living trees and standing deadwood of more than 4 m tall—using the equation volume = 0.000138 × (𝜋 × dbh)2 – 0.0066 × 𝜋 × dbh + 0.0955 for coniferous trees, and volume = 0.000203 × (𝜋 × dbh)2 – 0.0097 × 𝜋 × dbh + 0.1194 for broadleaf trees. In addition, we calculated the total TreM diversity (based on the 47 types described by Larrieu et al., 2018), and the density per hectare of trees bearing fruiting bodies of saproxylic fungi, tree injuries and exposed sapwood, rot-holes, and saproxylic TreMs known to be important for saproxylic beetles (Larrieu et al., 2018; Stokland et al., 2012). To account for the two successive processes of tree decline and death, we measured two gradients—i.e. stand decline and mortality. First, in all forest contexts, we quantified plot-level “stand mortality” as the proportion of standing dead tree basal area relative to the total basal area of both standing dead and living trees. Second, we quantified the plot-level “stand decline” using ARCHI protocol (Drénou, 2025). This method uses tree architecture to classify tree on a decline gradient, from healthy fully functioning tree to dead tree, including three decline stages in between—i.e. resilient, stressed, and irreversible death. We applied the method on the 20 closest dominant trees from the center of each plot in fir- and oak-dominated forests. ARCHI protocol has not been employed in spruce-dominated forests and therefore has not been calculated in this context (Table 1). We used the proportion of dying trees—i.e. stressed and irreversible death categories, therefore not adding dead trees—as a measure of “stand decline”. The use of both metrics provides a dynamic perspective where trees are first declining due to the disturbances and are ultimately dead. Along the manuscript, we specifically refer to stand decline and mortality as it is used by foresters, only focusing on trees. 3. Beetle sampling, identification and characterization We sampled saproxylic beetles using different types of traps in each study context. In fir- and spruce-dominated forests, we used flight-interception traps (©Polytrap). They consisted in crossed pair of transparent plastic panes (40 × 60 cm) above a funnel conducting into a container filled with an unbaited preservative (50% propylene glycol and 50% water with a drop of detergent). There were two interception-traps in fir forests, each at least 20 m from the other, placed around the center of each study plot (Cours et al., 2022), while there was only one interception-trap in spruce forest, placed at the center of each study plot (Cours et al., 2021). The traps were suspended roughly 1.5 m above the ground and were sampled every month from mid-May to mid-September 2016 in spruce forests and 2017 in fir forests. In oak forests, we used a combination of two types of traps in each plot, two green unbaited multi-funnel traps and one black multi-funnel trap baited with volatiles (a blend of cerambycid pheromones, ethanol and alpha-pinene, hooked to the trap, see Roques et al., 2023) to sample different families of saproxylic beetles (mostly Buprestidae beetles with the two unbaited green multi-funnel traps, and both Curculionidae and Cerambycidae species with the baited black multi-funnel trap; Sallé et al., 2020). The baited traps have an attractive effect decreasing with distance and mostly help capturing nearby individuals (Turchin & Odendaal, 1996). As in fir and spruce forests, collectors were filled with unbaited preservative. Each trap was installed in a separate tree at least 50 m from each other, around the plot center. In oak forests, the project objective was to sampled canopy fauna specifically (CANOPEE project). Therefore, these traps were suspended in the tree canopy, approximately 15 m above the ground. The traps were sampled every month from mid-May to the end of August 2020. All the captured saproxylic beetles were identified at the highest possible taxonomic level—i.e. 87% at the species level, 1% at the genus level and 12% at the family level. We retained the genera in which 100% of regional species are known to be saproxylic, and families in which at least 75% of regional species are known to be saproxylic (Bouget et al., 2019). Since we could not identify most of Staphylinidae species, and since their exclusion is not considered to have a significant impact on the analysis of the response of saproxylic beetles (Parmain et al., 2015), we excluded them from the data sets. We characterized each saproxylic beetle species by seven functional traits: the preferred (i) deadwood diameter and (ii) stage of decay, (iii) canopy closure, (iv) the mean body size, (v) their trophic regime, (vi) the used TreM, and (vii) whether they are flower-visitors or not. The three first indices were niche traits based on species occurrence data, i.e. class of deadwood diameter and stage of decay in which species were recorded (see Gossner et al., 2013; Janssen et al., 2017). Mean body size was based on direct measurements from the FRISBEE database (Bouget et al., 2019), as well as the larval feeding guild, the larval microhabitat, and the adult floricolous behavior. For further analysis, we specifically focused on “Xylophagous” and “Saproxylophagous” regarding the larval feeding guild, and on “Xylofungicolous” and “Cavicolous” species regarding the TreM used by the larva (Table 2). 4. Statistical analysis Data analysis was conducted with R software 4.5.1 (R Core Team, 2025). Regarding the response variables, we used the saproxylic beetle community matrices and the trait databases—i.e. the preferences in deadwood diameter, stage of decay and canopy closure, and the mean body size—to calculate the functional indices. First, we removed vagrant species that are only associated with broadleaf tree species in both coniferous-dominated forests (and vice versa for oak-dominated forests), based on Bouget et al. (2019). The functional richness (FRic), Rao’s entropy index (RaoQ), and evenness (FEve) were calculated using the FD R-package, implementing abundance community matrices (Laliberté et al., 2014). We also calculated the CWM and FDis of each of the quantitative traits used for the functional indices (i.e. preferences in deadwood diameter, stage of decay and canopy closure, and the mean body size). To isolate non-random assembly processes from the potential sampling effects associated with species richness, we calculated Standardized Effect Sizes (SES) for all functional diversity indices (FRic, RaoQ, FEve, FDis). We employed a null model approach, generating 999 randomized indices by randomly swapping trait values between species across each context species pool. The SES for each index was calculated as follows: SES = (Iobs − mean(Inull))/sd(Inull); where Iobs represents the observed index value, and mean(Inull) and sd(Inull) represent the mean and standard deviation of the null distribution, respectively. This standardization ensures that the resulting values reflect functional assembly (convergence or divergence) independent of taxonomic richness. In addition, we calculated the taxonomic diversity and abundance for the entire communities, but also for each ecological guilds (Bouget et al., 2019), defined by larval feeding guild (i.e. xylophagous, saproxylophagous, and mycetophagous), larval microhabitat (i.e. xylofungicolous and cavicolous), and adult behavior (i.e. floricolous, see Table 2). Finally, we used iNEXT R-package (Hsieh & Chao, 2020) to calculate standardized species richness—i.e. Hill number q = 0—therefore disentangling from the effect of abundance. We also included the calculation of Shannon diversity index for all species—i.e. Hill number q = 1. Due to inherent variations in sampling completeness across contexts and bird guilds, we standardized our comparisons using community-specific coverage levels (Table S2). Following the guidelines of Chao et al. (2014), we set the target coverage for each group at a level that required, at most, an extrapolation to double the observed sample size, thereby ensuring the reliability of our richness estimates. Consequently, one plot was excluded in spruce-dominated forests from the final analysis due to insufficient sampling completeness—i.e. sample coverage = 33% (see Table S1). Finally, to calculate the cascading effects of stand decline and mortality, and of salvage logging on environmental conditions and saproxylic beetles, we employed piecewise Structural Equation Modelling (SEM), using the piecewiseSEM R-package, version 4.1.2 (Lefcheck et al., 2020). In our SEMs, hierarchical relationships were the direct effects of stand decline and mortality, and salvage logging on trees (i.e. ligneous substrates) and subsequently on TreMs, which all may had effects on saproxylic beetles—i.e. the cascading pathway. We constructed SEM for each forest context, to account for the different sampling schemes. The models implemented in the SEM were a mix of linear mixed models, for which we log-transformed the variables that required it, and generalized linear mixed models for the abundance and species richness of saproxylic beetles, using respectively the negative binomial (after testing for residual over-dispersion, using “check_overdispersion” function from performance R-library, Lüdecke et al., 2020) and Poisson distribution families. We also added a random variable using site—i.e. Aure and Sault in fir context, and Orléans and Vierzon in oak context; Figure 1—, and the plot altitude as fixed covariable when significant in model to account for the within-site biogeographic context. We tested each model using the function “check_model” from performance R-package (Lüdecke et al., 2020), to test for residual normality, predictor collinearity, and potential outliers. We applied a selection of best SEM by the selection of predictors in each model that reduce overall SEM AIC by at least two. We divided the analysis in two datasets: (i) a first dataset for testing both the effects of stand decline (only in fir- and oak-dominated forests) and mortality, therefore not using salvage logged plots, and (ii) a second dataset for testing the effect of salvage logging, therefore using only salvage logged plots and unsalvaged plots from the most severe mortality class in each coniferous context—i.e. 10-45% in fir-dominated forests and 70-100% in spruce-dominated forests (Figure 1). In models, salvage logging has been used as a binary variable—i.e. absence–occurrence of salvage logging. In spruce-dominated forests, since we only tested the effect of either stand mortality or salvage logging—i.e. only one tested effect in each part of analysis—on ligneous substrates, we used the 5% α-error p.value to determine the significant effects. In fir- and oak-dominated forests, we used a 2.5% α-error p.value to test the effects of both stand decline and mortality effects as a strict adjustment of the p.value (0.05/2). For TreMs, we tested either stand mortality and decline, or salvage logging, as well as the ligneous substrate predictors (i.e. 14 or 15 tested effects in total). Likewise for biodiversity, we tested for the potential effects of all predictors (i.e. 23 or 24 tested effects). Therefore, we respectively adjusted the α-error p.value to 0.0036 or 0.0033 (as the result of 0.05/14 and 0.05/15), and 0.0022 or 0.0021 (as the result of 0.05/23 and 0.05/24). In addition, we also tested for the direct effect of stand decline and mortality, and salvage logging on the abundance of the species summing at least 10 individuals and occurring in at least 10% of the plots in each of the three contexts (Figure 5). We used the same model structure as in the SEM for the species abundance—i.e. the tested predictor, a random site factor in fir- and oak-dominated forests, and the negative binomial distribution family.