Effects of global change on bird and beetle populations 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 study how species community composition would change following anthropogenic and natural disturbances. We determine the main drivers of assemblage dissimilarity for bird and beetle communities.


| INTRODUC TI ON
Global climate warming affects the functionality of ecosystems by modifying forest composition and biomass, which in turn has repercussions for biodiversity and species assemblages (Hillebrand et al., 2010;Kelly & Goulden, 2008;Pachauri et al., 2014;Thuiller et al., 2011;Zhang & Liang, 2014). In fact, many studies have anticipated that anthropogenic radiative forcing will alter boreal biodiversity and ecosystems (Cadieux et al., 2020;Pachauri et al., 2014;Tremblay et al., 2018). For example, global-scale predictions have shown that the potential high emissions of greenhouse gases would lead mainly to negative effects on biodiversity (Pachauri et al., 2014). By contrast, other results have been published anticipating an increase in biodiversity during this century, which is referred to the northern biodiversity paradox (Berteaux et al., 2010;Matthews et al., 2004;Morin & Thuiller, 2009b). This paradox suggests that the ecological niches of some 'southerly' species would increase in size due to expansion beyond the northern periphery of their ranges (Berteaux et al., 2010).
In Canada, temperature has increased by 1.7°C since 1948, twice as fast as the global average (Bush & Lemmen, 2019). This increase in temperature could lead to the northward migration of thermophilous hard-wood tree species to the detriment of boreal conifers, particularly mid-to-late-successional species Boulanger & Puigdevall, 2021;Duveneck et al., 2014).
Moreover, climate change is expected to influence wildfire activity directly (Boulanger & Puigdevall, 2021), which would favour pioneer and fire-adapted boreal tree species .
Significant changes in species composition are expected within the transition zone between boreal and temperate biomes, where several tree species are currently reaching their thermal limits (Boulanger & Puigdevall, 2021;Brice et al., 2020). Global warming could also drive the occurrence of more extreme climatic events, including severe droughts (Kumar, 2013;Masson-Delmotte et al., 2018), thereby reducing the productivity of several boreal tree species through increasing metabolic respiration (Girardin et al., 2016). Climate-induced changes in insect outbreak regimes, notably those in spruce budworm (SBW; Choristoneura fumiferana), can lead to severe defoliation and death of firs (Abies spp.) and spruces (Picea spp.) and alter forest successional pathways (Labadie et al., 2021;Pureswaran et al., 2015).
Global change might influence biodiversity through behavioural, morphological and physiological changes in organisms or might even influence gene flow (Boulanger & Puigdevall, 2021;Matuoka et al., 2020;Micheletti et al., 2021;Naeem et al., 2012;Sergio et al., 2018;Wisz et al., 2013). The implications of global change on biodiversity could be synthesized through direct and indirect effects.
On one hand, the indirect effects are associated with habitat change (e.g. vegetation) that could occur through various events such as natural and anthropogenic disturbances . On the other hand, direct effects are mainly characterized by the influence of climate and meteorological conditions on the organisms (Micheletti et al., 2021). Both types of effects should be considered when studying the impacts of global warming on ecosystem functioning. For example, ecosystem services can depend on the magnitude of both direct and indirect pathways on global biodiversity.
Understanding both types of effects should thus help identify the best conservation actions that are required to maintain valued services. If habitat-based global change is projected to decrease species occupancy, land management actions could be adapted by targeting vegetation restoration. Otherwise, species translocation actions could be adopted if direct changes are anticipated that would cause a decline in species occupancy (Micheletti et al., 2021).
Climate and land-use changes may interact in complex ways to impact changes in wildlife habitats (Bentz et al., 2010;Tremblay et al., 2018). For example, an increase in disturbance rate due to warmer conditions could favour the regeneration of warmadapted, pioneer tree species (Boulanger & Puigdevall, 2021;Brice et al., 2020), which could affect biodiversity through adjustments in species ecological traits. Moreover, the synergistic effects of forest harvest and natural disturbances with direct climate change could alter species assemblages, and, therefore, biodiversity. In this work, we investigated future climate-induced variations in bird and beetle assemblages in Québec's boreal forest under different forest harvesting scenarios. This will enhance our knowledge regarding the interactive effects of climate change and forest use on future biodiversity. We assessed the effects of future climate conditions on these assemblages by comparing forest landscapes that were simulated under two anthropogenic radiative-forcing scenarios, that is, Representative Concentration Pathway (RCP) 4.5 and RCP 8.5 (van Vuuren et al., 2011) with landscapes that are simulated under average historical (baseline) climate conditions. Forest landscapes were simulated through a spatially explicit raster-based forest landscape model to capture key forest ecosystem processes. We further assessed how forest management may affect future assemblages by simulating species occupancy if current forest harvest levels are maintained vs there would be no harvest activities. We analyse the effects of global change in two ways: (1) indirectly, that is, habitatbased change (Boulanger & Puigdevall, 2021;Micheletti et al., 2021;Wisz et al., 2013) and (2) through a mixed effect combining indirect and direct effects, that is, those stemming from climate variables per se (Thuiller et al., 2018). The choice of the two pathways was made to quantify separately their implications on biodiversity in a cascading configuration. More precisely, starting first with habitat effects, through Habitat-Only-Based Models (HOBMs), and adding afterward meteorological conditions, through Climate-Habitat-Based Models (CHBMs). We also determined the main drivers of assemblage dissimilarity following the two pathways.

K E Y W O R D S
anthropogenic disturbance, assemblage dissimilarity, assemblage dissimilarity decomposition, biodiversity modelling, boreal forests, climate change, latitudinal gradient We opted for a species distribution modelling (SDM) approach (Guisan & Thuiller, 2005;Guisan & Zimmermann, 2000;Peterson et al., 2011) to model the single-species occurrence probability based on extensive field surveys of beetles and birds. We analysed the assemblage structure based on continuous occurrence probabilities (probability-based; Guisan & Thuiller, 2005;Grenié et al., 2020) to avoid overprediction risks (Dubuis et al., 2011;Gelfand et al., 2005;Grenié et al., 2020). It has been previously demonstrated that probability-based richness provided a better fit of the actual richness than would the threshold-based approach (Grenié et al., 2020). We computed dissimilarity measures (Albouy et al., 2012;Baselga, 2010Baselga, , 2013 between different climate scenarios under the two forest harvesting levels. In this context, we adapted continuous-based decomposition in the context of occurrence probabilities to detail the two components of Beta-diversity: balanced variation in species abundances and abundance gradient, which are generalizations of turnover and nestedness (Baselga, 2013). Quantifying the main causes of change in biodiversity could be very helpful in assessing the potential underlying determinants because species replacement and nestedness, for example, are two different ways of generating assemblage dissimilarity (Baselga, 2013).

| Study area and occurrence data
The study area is located in the Côte-Nord region of Québec, Canada Miller). Wildfires are the major natural disturbances in northern areas that have yet to be logged (Boucher et al., 2017). The southern part of the study area belongs to the eastern balsam fir-white birch (Betula papyrifera Marsha.) subdomain, mostly dominated by balsam fir and white spruce (Picea glauca [Moench] Voss) mixed with white (paper) birch. Forest harvesting had been the main source of forest disturbance since the late 1990s in this latter area (Bouchard & Pothier, 2011). In Québec, logging affected around 0.8% of public forest annually (Bureau du Forestier en chef, 2010). This part of the territory also has been affected by an outbreak of the SBW that began in 2006 and which is still ongoing. Tree mortality began around 2015 and has subsequently increased.
We used presence-absence data that were collected between 2004 and 2018 to model species distributions. Given that we were mainly focussing on the impacts of fire and forest harvesting, we wanted to remove sites that were located in stands heavily damaged by the SBW outbreak (Ministère des Forêts de la Faune et des Parcs, 2018) by using a cumulative index of defoliation severity from 2007 to 2018 (Labadie et al., 2021). Annual defoliation severity was based on aerial surveys characterizing damage that was incurred by SBW since 2006 (Ministère des Forêts de la Faune et des Parcs, 2018) and was classified between 0 and 3, with 3 indicating the highest level of defoliation. The cumulative severity of the outbreak was obtained by summing the annual severity in Labadie et al. (2021). Sites with cumulative severity values of 10 and above were discarded from analyses. This threshold attests for a strong impact of the SBW over a long period, which modifies stand structure and therefore affects species presence. Indeed, Labadie et al. (2021) showed (see Figure 2 in this paper) that, at this level, the response of shrubs is clearly stronger that of herbaceous plants (no overlap in the 95% CI), which indicates important changes in forest structure that modify habitat of many species.
F I G U R E 1 Study area with the occurrences of the two taxa.
We used the data from the Atlas of Breeding Birds (Atlas des oiseaux nicheurs du Québec, 2018), which were based on species occurrences that were detected using unlimited distance 5-min point counts (Bibby et al., 2000), which were collected during the breeding season (late May to mid-July) between 2010 and 2018. For beetles, we merged different databases that had been collected in (Bichet et al., 2016Janssen et al., 2009;Légaré et al., 2011). In addition, we used data from 54 sites that had been sampled in 2018 in the northern part of the territory, along the northeast principal road going to Labrador. The sampling protocols were characterized by one multidirectional flight-interception trap per site, which sampled flying beetles, and four meshed pitfall traps per site, which sampled epigeic beetles during their peak period of activity (June-August; Bichet et al., 2016;Janssen et al., 2009). For beetles, we used species-level identifications wherever possible; otherwise, we standardized the identification to the genus level (around 92% initial identifications were considered at the species level).

| Predictor variables
To predict species occurrence, we used climate and land cover variables that were grouped in two classes of models: Climate-Habitat-Based (CHBMs) and Habitat-Only-Based (HOBMs). These two model classes were designed to study climate-induced effects, as follows: CHBMs for the mixed effects and HOBMs for the indirect climate effects. Initially, we generated 22 potential climate variables using the BioSim platform (Régnière et al., 2017), including the annual average of temperature, precipitation and water deficit between 2004 and 2018 (see Table S1 for the description of all potential predictor variables). BioSIM simulates daily maximum and minimum temperatures, precipitation, water deficit, mean daily relative humidity and wind speed by matching georeferenced sources of daily weather data to spatially georeferenced points. BioSIM uses spatial regression to adjust weather data for differences in latitude, longitude and elevation between the sources of weather data and each field location (for more details, see Boulanger, Parisien, & Wang, 2018). In our case, the spatially referenced points were 10,000 points that can be located as far as 500 km from the simulated points which include New Brunswick, Ontario and Newfoundland-Labrador weather stations. The simulation of daily weather variables was carried out using BioSIM v11, a stochastic weather generator developed by Régnière and St-Amant (2007).
This tool utilizes georeferenced sources of weather data, that is, data from nearby weather stations, to match with a grid of geographical coordinates (here the 10,000 simulation points), taking into account factors such as latitude, longitude and elevation. The result is a simulation of daily maximum and minimum temperatures, F I G U R E 2 Distribution of the land cover and three climate variables under the simulated management scenarios in 2100. Panels (a, b, f, g, k, l) represent the distribution of land cover over the six study scenarios that were classified into four large cover classes: Natural, which included conifer dense, conifer open, mixed wood and open habitat; Fire and Cut for the land cover disturbed by fire and forest Harvest; and Others for the rest. The bar plots represent the relative frequency of each of the four habitats in the map. Panels (d-f, j-l, m-o) represent the distribution of mean temperature, precipitation and water deficit for the three climate scenarios Baseline, RCP4.5 and RCP8.5, respectively. Latitudes and longitudes are given in the UTM coordinate system. precipitation, mean daily relative humidity and wind speed that is restored from the monthly normals by introducing stochastic variation while maintaining daily temporal autocorrelation. The average elevation of each cell was calculated by intersecting the 30 arcsec Digital Elevation Model with the cell. Multiple spatial linear regression equations are used to calculate the local monthly gradients, used to adjust values for each weather station, by taking into consideration each variable on a monthly basis from 24 nearby stations in the normals database (located up to 500 km from the simulation point). Weather values for each simulated point are then distanceaveraged using the adjusted values of the four nearest weather stations. We generated climate variables by spatially interpolating data from the 10,000 random points using kriging and elevation as a drifting variable. Land cover maps from the Canadian National Forest Inventory (NFI) were used to generate land cover variables based on k-nearest-neighbour interpolation at a 250-m resolution that was referenced to the year 2001 (Beaudoin et al., 2014). To estimate forest composition, we used the relative importance of tree species groups (conifer and deciduous species), treed land and tree canopy-closure maps from these NFI data to generate five natural land cover classes: closed-canopy conifer forest; open-canopy mature conifer forest; mixed forest; open area; and others (Labadie et al., 2021; for more details, see Table S2). We used four of natural land cover categories together with six disturbance classes. For the disturbed stands (by fire or harvest), we subdivided them into three age classes: [0,10],]10,20] and]20,50]. years. In addition, we considered stand age and the distance to the nearest burned area as potential predictor variables. Stand age maps were based on the year 2001 (Beaudoin et al., 2014) with an update according to the registered fire and forest harvest disturbances between 2004 and 2018. Furthermore, it has been mentioned in Bichet et al. (2016) and Zhao et al. (2013) that the influence of the landscape varied between 400 m for beetles to 1000 m for birds. Consequently, we used a matrix of 20 pixels centred on the focal pixel (21 pixels in total) to calculate the relative frequencies of the land cover (see Tables S1 and S2 for description).

| Climate and forest management scenarios
We obtained future climate projections from the Canadian Earth System Model version 2 (CanESM2) for RCP 4.5 and RCP 8.5 for the period 2071-2100, which were further downscaled to a 10km resolution using ANUSPLIN (McKenney et al., 2013). Future monthly normals at each random point that was previously used to assess baseline climate were directly assessed from changes that were observed between the 1981 and 2010 period and future projections in the 10-km cell in which the random point was located.
Daily time series were stochastically generated for each random point from these future monthly normals using BioSIM. Future climate variables at each random point were calculated by averaging these daily values from 30 BioSIM simulations (Boulanger, Parisien, & Wang, 2018). Climate scenarios varied depending upon the mean annual temperature that was expected to increase between 3°C (RCP4.5) and 7.5°C (RCP8.5) throughout the southern boreal region by 2100 (compared to 2000, see Figure 2), while average precipitation was projected to increase between 7% and 10% regionally, with relatively small differences among scenarios (Boulanger & Puigdevall, 2021;Boulanger, Taylor, et al., 2018, see Figure 2).
The forest landscapes were simulated using the spatially explicit raster-based forest landscape model LANDIS-II (Scheller & Mladenoff, 2004), which simulates stand-and landscape-scale processes, including forest succession, seed dispersal and natural (wildfires and spruce budworm outbreaks) and anthropogenic (forest harvest) disturbances. This model has been extensively used in Québec over the last decade and has been thoroughly validated under various forest conditions Boulanger & Puigdevall, 2021;Taylor et al., 2017;Tremblay et al., 2018). Forest landscapes were initialized for the year 2000 conditions using the NFI attribute maps (Beaudoin et al., 2014) and provincial sample plots. Tree growth and regeneration as well as wildfires were climatesensitive in simulations. We simulated two levels of forest harvesting scenarios according to a gradient of harvesting pressure, from no harvesting to high-intensity harvesting, similar to current management practices in Québec (forest harvest applied to 8% of the harvestable upland area per 10 years). The Biomass harvest extension was used to simulate forest harvests (Gustafson et al., 2000). Only stands that contained tree cohorts that were greater than 60 years old were allowed to be harvested. Mean harvested patch size varied between 40 and 150 km 2 . Forest harvest rates were held constant throughout the simulations unless sufficient numbers of stands did not qualify for harvest. In the latter case, forest harvesting contin-

| The modelling strategy
We used two classes of models, that is, Habitat-Only-Based models (HOBMs) and Climate-Habitat-Based models (CHBMs), to study indirect climate effects and mixed climate effects, respectively. The mixed effects represented a combination of direct and indirect effects (see filled arrows in Figure 3). We divided our modelling strategy into five steps (see Figure 4). We aimed to study the effect of climate change under different forest harvesting scenarios. Therefore, we compared baseline climate scenario with RCP 4.5 and RCP 8.5 under two forest harvesting scenarios (NoHarvest and Harvest), hence the use of two baseline scenarios (see Step 5 in Figure 4).
The modelling procedure that is described here was repeated for the two model classes (CHBMs and HOBMs) according to potential predictor variables under each class. We standardized the variables to facilitate model convergence (MacKenzie et al., 2017) prior to model calibration that was based on the corresponding database. Steps 2 and 3 correspond to a cross-validation loop, in which we split the date into 10 folds of relatively equal size, so that at each step, 9 folds were used for training and the one remaining fold was used for testing.

| Step 1: Data preprocessing
We used the following procedure for preparing the two databases: (1) we removed all of the sites that were strongly affected by the spruce budworm outbreak; and (2) we included only the more common species, with a minimum record of 1% and 5% presence among sites for birds and beetles, respectively. These percentages made sure that we only modelled species that occurred at ≥24 sites for birds and ≥14 sites for beetles.

| Step 2: Model estimation
First, we removed the highly correlated variables, based on pairwise Pearson correlation coefficient (r), and retained the 5 most important predictor variables (Zurell et al., 2020). To do so, we fitted a univariate generalized linear model (GLM) with linear and quadratic terms for each variable; we ranked the predictors according to their importance, using the Akaike information criterion (AIC); we removed the highly correlated variables (|r| > 0.60). To reduce the separation in the regressions, we removed any predictor with a standard deviation value greater than 50 through a stepwise procedure.
We started with a generalized linear mixed model (GLMM; package 'lme4', Bates et al., 2015) with a random intercept to account for differences between sampling years. We developed six and three full potential models differing only in the fixed effects for CHBMs and HOBMs model classes, respectively (see Table S3). We used interaction terms between the best temperature variable (where the AIC criterion of the corresponding univariate regression is the lowest among all selected temperature variables) with distance to the nearest burned stand and with stand age variables to include the effect of latitudinal variation.
The best models were selected based on each full model under 10-fold cross-validation by minimizing the AIC criterion (package 'MuMIn', Barton, 2015) under the following conditions: (1)

| Step 4: Prediction
We estimated the parameters of the final model for the selected species based on the full dataset using the same procedure that was described for cross-validation. We subsequently used the simulated predictor variable maps for the six study scenarios and computed the occurrence probability maps.

| Step 5: Assemblage and occupancy analysis
We used the following indices to compare species assemblages between scenarios: Regional occupancy probability (ROP): ROP was calculated as the regional average of the occurrence probabilities for the study area (Bichet et al., 2016). We also used the percentage of change in (1)

F I G U R E 4 Modelling 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 cross-validation. 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 AUC ≥0.7 and computed the prediction maps for the six specified scenarios. The last step was devoted to the assemblage analysis.
RBC Sim,Ref and SBC Sim,Ref,j represented, respectively, the regional Bray-Curtis dissimilarity between the scenarios Sim and Ref, the spatial Bray-Curtis at the pixel j selected randomly.
From Figure 5, we compared the performance of the Jaccard index that was based on our continuous output (SJ OP ) and the traditional Jaccard index that was based on the presence-absence transformation (SJ Incidence ). From the simulation results, the SJ OP yields results that were close to SJ Incidence ; moreover, performance increased with the number of species that were analysed (see Figure 5i). Furthermore, we added two situations with weaker binarization (Bad and Medium cases in Figure 5) to illustrate that a gap can be generated between the two indices that is mainly due mainly to binarization error.
Beta ratio: We separated the Bray-Curtis dissimilarity index into two additive components main drivers in assemblage dissimilarity: (1) the occurrence gradient (BC_grad); and (2) balanced variation in species occurrence (BC_bal). We used these notations instead of the abundance gradient and balanced variation in species abundance (Baselga, 2013), given that we worked with occurrence and occupancy probabilities. We presented a detailed explanation for each component, as follows: Occurrence Probabilities Jaccard dissimilarity. The three levels of binarization concerned the rank choice of the occurrence probabilities that were simulated from a mixture of uniform distributions under constraints. Panels (a-c) represented the three simulated occurrence probabilities and the corresponding binarized maps with a 0.5 threshold for one species and under one scenario. Panels (d-f) showed the frequency of the occurrence probability at the different pixels. Panels (g-i) concerned the absolute difference between Jaccard dissimilarity based on the binarized data and the continuous version by varying the number of species.
in occurrence between the two scenarios is too large and in the same direction in almost all species.
To measure the fraction of each component of Bray-Curtis dissimilarity, we used the BC_ratio given by If BC _ ratio Sim,Ref > 0.5, the assisted community change was caused mostly by the occurrence gradient, whereas a value smaller than 0.5 indicated a dominance of balanced change in occurrence in the assemblage dissimilarity (Albouy et al., 2012). We used the package 'betapart' (Baselga & Orme, 2012) for the assemblage analysis.
The BC_ratio and its components were illustrated by using occurrence probabilities. The same formalism was adapted regionally through the regional occupancy probabilities.
For bird species, stand age was the most frequently selected predictor variable with 43.9% and 71% of selection in CHBMs and HOBMs, respectively (see Figure 6). For beetle species, mean temperature of the warmest month (WarmMTmean) and stand age were the most frequently selected variables under CHBMs and HOBMs with 33.7% and 40.3% of selection, respectively (see Figure 6).

| Occurrence and regional occupancy
The analysis of the occurrence and occupancy helped us to evaluate the direction and magnitude of potential community changes following global change. The difference between the indirect and mixed climate effects was visualized through ∆ROP Sim,Ref (Figure 7a).
The magnitude of changes in occupancy compared with the baseline reference scenario was larger when we included the climate variables (CHBMs) when compared to the models including only habitat variables (HOBMs). The change in occupancy was slightly more significant for birds, compared with beetles under HOBMs. Under CHBMs, we observed that occupancy increased for birds but decreased for beetles when comparing the same radiative-forcing scenarios (RCP4.5 and RCP8.5) to the baseline (the projection maps also demonstrated this result; compare Figure 8a,b,e,f). However, we observed a decrease in the percentage of winner species with climate change for the two taxa, except for birds under HOBMs with no forest harvest (Figure 9; see Figure S1 responses of five variable classes on eight winner and loser species). Moreover, to observe the change at the species level, File S1 contained the percentage of change in

F I G U R E 6
Percentage of the predictor variables that were included in the species regressions for the two taxa (See Table S1

| Regional species assemblage change
An increase in assemblage dissimilarity was observed when comparing RCP4.5 and RCP8.5 to the baseline for 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 forest harvest (Figure 7b,c). Also, based on HOBMs, regional dissimilarity from RCP4.5 to RCP8.5 increased, respectively, by 0.07 and 0.08 for birds and beetles under no forest harvest. This regional dissimilarity was mainly incurred through balanced variation in species occupancy for both taxa (BC_ratio <0.5; see Figure 7a).
We also observed an increase in assemblage dissimilarity from HOBMs to CHBMs for both taxa. Inclusion of the climate variables reshaped assemblage structure strongly under the two forest management levels (see differences in the Jaccard dissimilarity index between HOBMs and CHBMs, Figure 7b,c).

| Latitudinal changes in species assemblage
The observed, compared to that of beetles (see the difference between Figure 10a,c).
Our models predict that beetles would show greater sensitivity to climate variations, given that an increase in dissimilarity was observed even for a medium level of climate change (Figure 10a), that is, RCP 4.5 following changes in temperature between the baseline and the two forcing scenarios. In Figure 11a,b, we depicted spatial change in Tmin and WarmMTmean between RCP8.5 and baseline climate scenarios. Tmin and WarmMTmean were the most frequently selected respective predictor variables for birds and beetles (see Figure 6b for the selection percentage). Furthermore, the observed latitudinal dissimilarity gradient was derived mainly from balanced variation in occurrence (BC_ratio <0.5 in Figure 10b,d).
Yet, the two taxa behaved in a contrasting manner with regard to their latitudinal variation in the occurrence gradient component of BC dissimilarity.

| DISCUSS ION
The present study reveals several intricate effects of global change on animal species assemblages. We showed that the two climateinduced pathways studied increased dissimilarity in species assemblages mostly by inducing species turnover, under both forest harvesting scenarios. Of the two pathways, the models integrating both indirect (habitat) and direct climate effects modified the assemblages more than did the models based only on indirect (habitat) effects. In fact, the difference in magnitude between the two F I G U R E 9 Percentage of winner species per taxon under the four landscapes that were compared and the two model classes (CHBMs and HOBMs). Species i was considered a winner if the regional occupancy probability under the simulated scenario was higher than under the reference scenario (ROP Sim,i > ROP Ref,i  Nevertheless, we observed almost an opposite feedback between the two taxa regarding changes in the regional occupancy. This could result from physiological differences between taxa in their individual response to climate. Our analysis predicts pronounced variations in assemblage composition for both birds and beetles by 2100 that will be caused by climate change under any of the forest harvest levels. Whereas a rise in anthropogenic radiative forcing should increase dissimilarity in species assemblage following each of the two climate pathways, the magnitude of change differed considerably between the two pathways. For example, the minimum and the maximum assemblage dissimilarities were 0.49 and 0.14 under climate-habitat and habitatonly models, respectively, for both taxa. This result could explain the subsequent mismatch between the climate and the biota (Micheletti et al., 2021;Stralberg et al., 2015;Wu et al., 2015). For example, Brice et al. (2020) noted that under climate change the variation in climate conditions would be faster than the capacity of tree species to migrate, a discrepancy that is expected to create a gap between the climatic niche and the observed distribution of species.
Regarding the direction of change at the taxon level, we observed that birds and beetles responded in almost opposite directions. In fact, the increase in bird occupancy that resulted from the mixed combination between direct and indirect climate effects coincided in overall with the northern biodiversity paradox that was emphasized by Berteaux et al. (2010), Matthews et al. (2004), and Morin and Thuiller (2009b), which anticipated an increase in biodiversity in northern protected areas during this century. More precisely, the northern biodiversity paradox suggests that climate change could lead to a potential increase in the abundance of many species for which the low temperatures are currently a limiting factor (Berteaux et al., 2010). For beetles, a substantial decrease in regional occupancies is predicted following an increase in climate change. This outcome aligned with the global scale biodiversity trajectory predicting mostly negative effects of high emissions of greenhouse gas on biodiversity and ecosystem services (Pachauri et al., 2014). Given that beetles are poikilotherms, their internal The authors demonstrated that insects found at middle elevations and on mountain tops were less tolerant to temperature increases than were species located at lower latitudes. The observed latitudinal gradient was also resulted from a gradual change in forest composition linked with a northward range shift of species. A projected change in the composition and structure of boreal forest tree species (Boucher et al., 2017;Schneider et al., 2016) is expected to shift ecological niches and this will change local biodiversity and species interactions (Kerr & Packer, 1998;Morin & Thuiller, 2009a).

F I G U R E 11
Temperature difference maps between RCP8.5 and baseline of annual minimum temperature (a) and mean temperature of the warmest month (b) in °C. Map scale: The dark red and the dark blue represented, respectively, the highest and the lowest values. The horizontal and vertical marginal plots, respectively, show the longitudinal and latitudinal distribution of the temperature differences.
Under mixed climate effects, we observed substantial latitudinal changes in biodiversity with increasing radiative forcing, with an almost complete change in beetle assemblages in the northern portion of the study area. Latitudinal patterns were also noted by Brice et al. (2019), who reported a northward decrease in temporal beta diversity of tree species between past (1970-1980) and present (2000-2016) periods. Our results warn about the possibility that direct climatic effects could increase over the next several decades, particularly for insects. This latitudinal gradient could be a result of the increase in temperature after a change in global climate forcing (Holland & Bitz, 2003 We must mention some points regarding our approach when interpreting the results. (1) Our direct and indirect effect results were based on the simultaneous effect of the meteorological and land cover indicators and did not consider the time-lag effect of the climate factors on vegetation. Therefore, we think that the passage from indirect to mixed climate effects under the case RCP8.5 could cover a good range of the potential future situations under climate change.
(2) The climatic envelope used for calibration of the studied organisms was based on the data present within the study area.
Among the studied species, many have a much larger distribution than that limited to the study area, and thus, future work should intend to broaden the climatic envelope used for calibration by working with a subset of species with the available data.
(3) Our modelling approach ignored the biotic interactions and was based on implicit assumptions of individualistic species responses to climate change.
The biotic interactions could influence the assemblage structure through the species-environmental relationship but also by generating nonrandom species co-occurrence pattern (Wisz et al., 2013).
Therefore, birds for instance could be affected by the decrease in beetle population (and possibly other types of insects) on which they feed. This is particularly important in the context of recent insect declines reported in several parts of the world (Harris et al., 2019;Outhwaite et al., 2022;Wagner, 2020).
In conclusion, based on implicit assumptions regarding individual species responses to climate change, our analyses identified potential repercussions of two climate-induced pathways that were based on direct and indirect effects on the assemblage composition of two taxonomic groups regionally and latitudinally, under two different forest harvesting scenarios. The observed change was derived mainly by species turnover regardless of forest harvest level. Moreover, based on the studied species pool and despite the possibility that other species could arrive, we could expect large differences in occupancy between the two studied taxa, which are explained probably by their ecological trait differences. This could indicate the potential range of change in boreal species concerning novel environmental conditions.

CO N FLI C T O F I NTE R E S T S TATE M E NT
We have no competing interests to declare.

DATA AVA I L A B I L I T Y S TAT E M E N T
Beetles occurrence datasets as well as the used predictor variables for the analysis are available from the datad ryad.org repository: https://doi.org/10.5061/dryad.ghx3f fbsb. Birds datasets were provided by BirdsCanada (https://www.birds canada.org).