Long-term effect of forest harvesting on boreal species assemblages under climate change

Logging is the main human disturbance impacting biodiversity in forest ecosystems. However, the impact of forest harvesting on biodiversity is modulated by abiotic conditions through complex relationships that remain poorly documented. Therefore, the interplay between forest management and climate change can no longer be ignored. Our aim was to study the expected long-term variations in the assemblage of bird and beetle communities following modifications in forest management under different climate change scenarios. We developed species distribution models to predict the occurrence of 87 species of birds and beetles in eastern Canadian boreal forests over the next century. We simulated three climate scenarios (baseline, RCP4.5 and RCP8.5) under which we varied the level of harvesting. We also analyzed the regional assemblage dissimilarity by decomposing it into balanced variations in species occupancy and occupancy gradient. We predict that forest harvesting will alter the diversity by increasing assemblage dissimilarity under all the studied climate scenarios, mainly due to species turnover. Species turnover intensity was greater for ground-dwelling beetles, probably because they have lower dispersal capacity than flying beetles or birds. A good dispersal capacity allows species to travel more easily between ecosystems across the landscape when they search for suitable habitats after a disturbance. Regionally, an overall increase in the probability of occupancy is projected for bird species, whereas a decrease is predicted for beetles, a variation that could reflect differences in ecological traits between taxa. Our results further predict a decrease in the number of species that increase their occupancy after harvest under the most severe climatic scenario for both taxa. We anticipate that under severe climate change, increasing forest disturbance will be detrimental to beetles associated with old forests but also with young forests after disturbances.

: The studied scenarios. We considered three of them for reference (Climate=Baseline, Har-vest=NoHarvest; Climate=RCP 4.5, Harvest=NoHarvest; Climate=RCP 8.5, Harvest=NoHarvest) and compared them to six other scenarios. Under each climate scenario (Baseline, RCP 4.5 and RCP 8.5) we compared NoHarvest to Harvest 0.5 and Harvest, which makes six comparisons in total.   Table 2: Natural habitats used (if less than 50% of the pixel's surface is in disturbance).
Habitat Condition Dense conifer -vegetation exceeds more than 50% of the surface -canopy density is higher than 60% of the surface -coniferous stands occupy more than 75% of the surface Open conifer -vegetation exceeds more than 50% of the surface -canopy density is less than 60% of the surface -coniferous stands occupy more than 75% of the surface Mixed-wood -vegetation exceeds more than 50% of the surface -neither coniferous nor broad-leaf trees account for more than 75% of the surface -the vegetation surface with trees is greater than the surface without trees Open habitat -vegetation exceeds more than 50% of the surface -the vegetation surface with trees is less than the surface without trees Other -Non-vegetation (trees, shrubs, grasses, bryophytes) that exceeds more than 50% of the surface Therefore, we used for ground-dwelling beetles 20 pixels (231 m resolution) centered on the focal pixel (i.e., 140 21 pixels in total) to calculate the frequencies of the land cover classes with a radius of 2.5 pixels (∼0.58 km).

141
For flying beetles and birds, we used 37 pixels with a radius ∼0.81 km to calculate land cover frequencies. 142 We also considered the distance to the nearest burned areas as well as the stand age as potential predictor 143 variables.

144
2.5 Association between bird species and habitat type 145 We considered three land cover classes to calculate species winning conditional probability (WCP), i.e., an 146 increase in the probability of species occupancy from uncut to post-harvest in 2100, conditional to species 147 habitat (see Eq. 2). We calculated WCP for (1) generalist bird species; (2) species associated with early-to-148 mid succession forest (EMSF), which included the following land covers: wetland, deciduous, mixed-wood 149 (young to mid-age), and coniferous habitat (young to mid-age); and (3) species associated with late succession 150 forest (LSF), which included the following land covers: mixed-wood (mature to old), deciduous (mature to 151 old) and coniferous habitat (mature to old). Stand age was classified as follows: young (stand age < 30 152 years), mid-age (stand aged between 30 and 60 years), mature (stand aged between 60 and 80 years), and old 153 (stand age greater than 80 years). This classification was used only for birds because information regarding 154 beetles' associated land cover was unavailable at the moment of the study. Bird-associated land cover has 155 been based on habitat associations from Gauthier and Aubry (1995) and Robert (2019). trap. This ratio provides an index of beetle proficiency to fly or to walk on the ground. As the ratio provides 164 continuous estimates, we kept only beetles with a ratio higher than 2 as flying beetles and beetles with a 165 ratio lower than 3 as ground-dwelling beetles. Beetle species with ratios between 2 and 3 were removed from 166 the study. intercept was added to include any differences between sampling years based on three full potential models 174 deferring only in the fixed effects of the models. We partitioned the modelling strategy into several steps 175 starting from data pre-processing to the species occurrence probability projections. SDMs were estimated 176 only for species recorded at least at 1% and 5% of sampling sites for birds and beetles, respectively. We used 177 land cover variables extracted from LANDIS forest landscape simulations for the simulated scenarios (see 178   Table S1 for predictor variables description). We also standardized the predictor variables just before any

184
We suggested three full potential models, the first one contained only linear forms of the predictor 185 variables. However, the second and third full models included also quadratic transformations of the covariates 186 (see Table S2). The best models were selected for each full model based on a 10-fold cross-validation procedure where ROP s denotes the average of the ROP over all the species under the scenario s and is given by and where N species , N pixel and P s,i,j respectively represent the number of species, the number of pixels in 194 the study area, and the probability of occurrence of the species i for the scenario s at the cell j.

195
Winning conditional to habitat or to relation to disturbance probability was calculated as We split the Bray-Curtis dissimilarity index into two additive components: where BC_grad Sim,Ref is the occupancy gradient and BC_bal Sim,Ref is the species balanced variation in occupancy. The latter was defined as To quantify the main source of assemblage dissimilarity, we performed a decomposition of the regional Bray- selected for the projection with AU C ≥ 0.7 (Fig. 2a-b-c) respectively. Stand age was the most frequently 215 selected predictor variable for the two taxa (Fig. 2d). Moreover, ground-dwelling beetles were more closely 216 associated with old undisturbed forests than flying beetles. However, the response of flying beetles was much 217 greater to the distance from burned forests than ground-dwelling-beetles, which showed a weaker response 218 that decreased with time after wildfire (Fig. 2d).   beetles, flying beetles had lower dissimilarity in assemblages regardless of climate scenarios (Fig. 3). However, 234 the dissimilarity in species assemblage over harvest gradient decreased with climate change (see the difference 235 in Bray-Curtis dissimilarity between baseline and RCP 8.5 in Fig. 3). This result was also observed in the 236 structure of land cover under RCP8.5 compared to RCP4.5 and baseline by varying harvest intensity (see 237 Fig S4).
Figure 3: Regional Bray-Curtis dissimilarity measures (Eq. 3) computed based on the difference in species regional occupancy between reference and simulated scenario. At each climate change scenario (Baseline, RCP 4.5, RCP 8.5), we computed the dissimilarity between Harvest 0.5 and Harvest vs. NoHarvest in 2100. Under Baseline, we compared respectively Baseline-NoHarvest to Baseline-Harvest 0.5 and Baseline-NoHarvest to Baseline-Harvest. Under RCP 4.5 we compared respectively RCP 4.5-NoHarvest to RCP 4.5-Harvest 0.5 and RCP 4.5-NoHarvest to RCP 4.5-Harvest. Finally, under RCP 8.5 we compared respectively RCP 8.5-NoHarvest to RCP 8.5-Harvest 0.5 and RCP 8.5-NoHarvest to RCP 8.5-Harvest. The abbreviation "None" represents the case of species associated to closed forests.

Occupancy analysis 239
We observed a change in the regional occupancy between the scenarios of forest harvesting (Harvest 0.5 and 240 Harvest vs. NoHarvest) under all three climate change scenarios. More precisely, we observed an increase 241 in regional occupancy in the bird community from nonharvest to harvest (passage from 13% to 15% in 242 the regional occupancy from nonharvest to harvest scenarios under RCP 4.5). The opposite behaviour was 243 observed in the beetle community (for instance a decrease in regional occupancy from 31% to 27% was 244 observed for flying beetles by comparing the same scenarios). These changes in the overall occupancy of 245 birds, and the decrease for beetles from NoHarvest to Harvest cases was visible over 114,118 km 2 of the 246 study area (Fig. 4) and also regionally (Fig. 5a). Post-harvesting changes in species assemblages were mainly 247 caused by species turnover (BC_ratio Sim,Ref < 0.5) for the two taxa with a strongest turnover for ground-248 dwelling than flying beetles. Under each climate scenario, the occupancy of bird species associated with young 249 habitats together with generalist species increased significantly contrary to species associated with mature 250 and old forests, when we compared harvested to unharvested scenarios (Fig. 5b). Under baseline and RCP 251 4.5, 50% of bird species associated with late succession forest increased in occupancy from no harvested to 252 harvested scenarios (Fig. 6a), whereas 82% of species associated with early to mid succession forest increased 253 in occupancy (see Fig S2 for the variation at the species scale). By comparing RCP 8.5 to RCP 4.5, the 254 proportion of winner species (from non-harvested to harvested landscapes) that were associated with early 255 to mid-succession forests decreased from 82% to 73%. Under the same conditions, species associated with 256 late succession forests increased from 50% to 63% (Fig. 6a). When we compared non-harvested to harvested 257 scenarios for flying beetles, species that were associated with burned stands had the weakest increase in 258 occupancy compared to species associated with harvested stands or closed forests (Fig. 5c).

259
The percentage of winning flying beetles associated with closed forests decreased from 53% to 47% from 260 baseline to RCP 8.5 under a forest harvest scenario (Fig. 6b). Winning species of flying beetles associated 261 with harvested stands even decreased from 42% to 33% under the same conditions. Among flying beetles 262 associated with burned forests, the percentage of winning species went from 40% to 0%. Regarding ground-263 dwelling beetles, the percentage of winning species associated with closed forests decreased from 64% to 57% 264 under the harvest scenario under climate change (baseline to RCP 8.5) (Fig. 6c).  It is important to underline that the effects of forest harvesting were studied under each climate scenario.

322
We showed that species turnover was the dominant source of dissimilarity in species assemblages for the 323 two studied taxa but ground-dwelling beetles registered the highest rate. This might be due to their lower 324 moving capacity. They cannot escape easily new conditions generated by disturbances. As seen earlier, they interactions are critical, and they might be the first to disappear in such turnovers.

334
In comparison, beetle communities that fly to colonize burned trees after wildfire are much more struc- level. On the long-term, this lowers landscape heterogeneity as well as biodiversity at the landscape level.

352
As observed for ground-dwelling beetles, changes in occupancy are also most visible for predators in flying 353 beetles. While the proportion of "losers" increases from 50 to 63% between climate change scenarios for 354 predatory ground dwelling beetles, it increases from 44 to 77% for predatory flying beetles. The scenarios 355 tested integrate burned forests interacting with increased harvesting and the resulting landscapes suggest 356 that an erosion in the flying beetle food chain is likely to impact all trophic levels.

357
Our results were based only on the most common species, the response of rare specialists and threatened 358 species might decrease the regional occupancy with both harvest activities and climate conditions (Cadieux 359 et al., 2019). Therefore, it will be interesting to analyze the response of the rare species to disturbance. One

362
In conclusion, forest harvesting activities are predicted to modify the composition of beetle and bird 363 species assemblages under all climate scenarios. We also indicated that the consequences of maintaining the 364 current level of forest harvesting under the worst climate change scenario could modify the species-habitat 365 association. Ideally, we aim to generate a situation where both early and late successional stages can coexist 366 by managing forests. Interestingly, we observed opposite responses in occupancy between birds, which was 367 generally positive, and beetle, which was negative, in our study area. The perfect separation between the 368 two taxonomic groups regarding the change in occupancy could reflect the impact of ecological traits such 369 as movement capacity and diet differences between boreal species. In any case, forest management must be 370 adapted to achieve conservation objectives under all climate scenarios. Our simulations show that harvest 371 activities must be reduced. Even when it maximizes biodiversity, harvest could be harmful to taxa with low 372 moving capacity (beetles in our case) which characterize old-growth forests.

373
Supporting information 374 File S1. The studied beetle and bird species.

375
Table S1 Description of the used predictors.
376 Table S2 The potential 3 full models.