Effects of global change on animal biodiversity in boreal forest landscape: an assemblage dissimilarity analysis

Despite an increasing number of studies highlighting the impacts of climate change on boreal species, the main factors that will drive changes in species assemblages remain ambiguous. We quantify two climate-induced pathways based on direct and indirect effects on species occupancy and assemblage dis-similarity under different harvest management scenarios. The direct climate effects illustrate the impact of climate variables while the indirect effects are reflected through the changes in land cover composition. To understand the main causes in assemblage dissimilarity, we analyze the regional and the latitudinal species assemblage dissimilarity by decomposing it into balanced variation in species occupancy and occurrence and occupancy and occurrence gradient. We develop empirical models to predict the distribution of more than 100 bird and beetle species in the Côte-Nord region of Québec over the next century. Our results show that the two pathways are complementary and alter biodiversity, mainly caused by balanced variation in species occupancy and occurrence. At the regional scale, both effects have an impact on decreasing the number of winning species. Yet, responses are much larger in magnitude under mixed climate effects (a mixture of direct and indirect effects). Regional assemblage dissimilarity reached 0.77 and 0.69 under mixed effects versus 0.09 and 0.10 under indirect effects for beetles and birds, respectively, between RCP8.5 and baseline climate scenarios when considering harvesting. Therefore, inclusion of climatic variables considers aspects other than just those related to forest landscapes, such as life cycles of animal species. Latitudinally, assemblage dissimilarity increased following the climate conditions pattern. Our analysis contributes to the understanding of how climate change alters biodiversity by reshaping community composition and highlights the importance of climate variables in biodiversity prediction.

To predict species occurrence, we used climate and land cover variables that were grouped in two classes 136 of models: Climate-Habitat-Based (CHBMs) and Habitat-Only-Based (HOBMs). These two model classes 137 were designed to study climate-induced effects, as follows: CHBMs for the mixed effects; and HOBMs for 138 the indirect climate effects. Initially, we generated 22 potential climate variables at a 250-m resolution using  Table. 1 in Supporting Information for the description of all 141 potential predictor variables). BioSIM simulates daily maximum and minimum temperatures, precipitation, 142 water deficit, mean daily relative humidity and wind speed by matching georeferenced sources of daily 143 weather data to spatially georeferenced points. BioSIM uses spatial regression to adjust weather data for 144 differences in latitude, longitude, and elevation between the sources of weather data and each field location 145 (for more details, see (Boulanger et al., 2018a)). In our case, the spatially referenced points were 15 000 146 points that were randomly located across the entire Province of Québec, whereas weather data were daily 147 data originating from discrete weather stations that were located within the province. We generated climate 148 variables at a 250 m scale by spatially interpolating data from the 15 000 random points using kriging and 149 elevation as a drifting variable. Land cover maps from the Canadian National Forest Inventory (NFI) were 150 used to generate land cover variables based on k-nearest-neighbor interpolation at a 250-m resolution that  varied between 400 m for beetles to 1000 m for birds. Consequently, we used a matrix of 20 pixels centred 162 on the focal pixel (21 pixels in total) to calculate frequencies of the land cover (see Table.1 and Table. 2 in   163 Supporting Information for descriptions of the 10 land cover variables that were used in the study).

177
The forest landscapes were simulated using the spatially explicit raster-based forest landscape model

188
Only stands that contained tree cohorts that were greater than 60-years-old were allowed to be harvested.   199 We used two classes of models, i.e., Habitat-Only-Based models (HOBMs) and Climate-Habitat-Based mod-200 els (CHBMs), to study indirect climate effects and mixed climate effects, respectively. The mixed effects 201 represented a combination of direct and indirect effects (see filled arrows in Fig. 3). We divided our modeling 202 strategy into five steps (see Fig. 4). We aimed to study the effect of climate change under different forest 203 harvesting scenarios. Therefore, we compared baseline climate scenario with RCP 4.5 and RCP 8.5 under 204 two forest harvesting scenarios (NoHarvest and Harvest), hence the use of two baseline scenarios (see Step 205 5 in Fig. 4). 206 The modelling procedure that is described here was repeated for the two model classes (CHBMs and

2.4.1
Step 1: Data pre-processing 213 We used the following procedure for preparing the two databases: (1) we removed all of the sites that were 214 strongly affected by the spruce budworm outbreak; and (2) we included only the more common species, with 215 Figure 4: The modeling strategy that was used under Habitat-Only-Based models (HOBMs) and Climate-Habitat-Based models (CHBMs).
Step 1 was dedicated to database preparation before starting the crossvalidation. In Step 2, we performed the cross-validation and computed the performance metrics in Step 3. In Step 4, we calibrated SDMs only for the selected species AU C ≥ 0.7 and computed the prediction maps for the six specified scenarios. The last step was devoted to the assemblage analysis. a minimum record of 1% and 5% presence among sites for birds and beetles, respectively. These percentages 216 made sure that we only modeled species that occurred at ≥ 24 sites for birds and ≥ 14 sites for beetles. Generalized linear model (GLM) with linear and quadratic terms for each variable; we ranked the predictors 221 according to their importance, using the Akaike information criterion (AIC); we removed the highly correlated 222 variables (|r| > 0.60). To reduce the separation in the regressions, we removed any predictor with a standard 223 deviation value greater than 50 through a stepwise procedure. 224 We started with a generalized linear mixed model (GLMM) (package 'lme4', (Bates et al., 2015)) with a 225 random intercept to account for differences between sampling years. We developed six and three full potential 226 models differing only in the fixed effects for CHBMs and HOBMs model classes, respectively (see Table 3 in 227 Supporting Information). We used interaction terms between the best temperature variable (where the AIC  Step 4: Calibration and prediction 247 We estimated the parameters of the final model for the selected species based on the full dataset using 248 the same procedure that was described for cross-validation. We subsequently used the simulated predictor 249 variable maps for the six study scenarios and computed the occurrence probability maps. We used the following indices to compare species assemblages between scenarios: where N species , N pixel and P s,i,j represented respectively the species number, the number of pixels in the 258 study area and the occurrence probability of the species i for the scenario s at the cell j. dissimilarity measure, regional (RJ) and spatial (SJ) Jaccard indices computed as follows: , RJ performance increased with the number of species that were analyzed (see Fig. 5I). Furthermore, we added 269 two situations with weaker binarization (Bad and Medium cases in Fig. 5) to illustrate that a gap can be 270 generated between the two indices that is mainly due mainly to binarization error.  The BC_ratio and its components were illustrated by using occurrence probabilities. The same formal-290 ism was adapted regionally through the regional occupancy probabilities. respectively (see Fig. 6). .   Table 5 to Table 9 in the Supporting Information contained the percentage of change in the 312 Figure 6: Percentage of the predictor variables that were included in the species regressions for the two taxa (See Table.    3.2 Regional species assemblage change 314 An increase in assemblage dissimilarity was observed when comparing RCP4.5 and RCP8.5 to the baseline for 315 both taxa. Based on CHBMs, regional dissimilarity from RCP4.5 to RCP8.5 increased respectively by 0.20 and 0.14 for birds and beetles under no harvest (Fig. 7B-C). Also, based on HOBMs, regional dissimilarity 317 from RCP4.5 to RCP8.5 increased respectively by 0.07 and 0.08 for birds and beetles under no harvest. This 318 regional dissimilarity was mainly incurred through balanced variation in species occupancy for both taxa 319 (BC_ratio < 0.5) (see Fig. 7A). 320 We also observed an increase in assemblage dissimilarity from HOBMs to CHBMs for both taxa. Inclusion 321 of the climate variables reshaped assemblage structure strongly under the two forest management levels (see 322 differences in the Jaccard dissimilarity index between HOBMs and CHBMs, Fig. 7B-C). The inclusion of climate variables produced a latitudinal gradient in projections of assemblage dissimilarity. 325 We used CHBMs to analyze how latitudinal changes in temperature and other climate variables would affect 326 assemblage structure. A clear increasing pattern was observed in assemblage dissimilarity heading north, 327 especially for beetle species. For birds, a slight increase in dissimilarity was observed, compared to that of 328 beetles (see the difference between Fig. 9A and Fig. 9C).

339
The present study reveals several intricate effects of global change on animal species assemblages. We showed 340 that the two climate-induced pathways that were studied increased dissimilarity in species assemblage mostly 341 by inducing species turnover, under both forest harvesting scenarios. Of the two pathways, the models 342 integrating both climate and habitat effects strongly modified the assemblage composition more than did 343 the models based only on habitat effects. In fact, the difference in magnitude between the two effects was 344 due to the mismatch that was generated by rapid climate variation compared to slower vegetation change  Table. 6-9 350 in Supporting Information). Nevertheless, we observed almost an opposite feedback between the two taxa 351 regarding changes in the regional occupancy. This could result from physiological differences between taxa 352 in their individual response to climate.

353
Our analysis predicts pronounced variations in assemblage composition for both birds and beetles by  https://cran.r-project.org/web/packages/lme4.