Aim

Climate change is expected to influence boreal bird communities significantly, notably through changes in forest habitat (composition and age structure), in the coming decades. How these changes will accumulate and interact with anthropogenic disturbances remains an open question for most species.

Location

Northeastern Alberta, Canada.

Methods

We used the LANDIS-II forest landscape model to project changes in forest landscapes, and associated bird populations (72 passerine species), according to three climatic scenarios (baseline, RCP 4.5, RCP 8.5) and three forest harvesting scenarios of differing intensity.

Results

Both forest harvesting and climate-related drivers were projected to have large impacts on bird communities in this region. As a result of climate-induced increases in fire activity as well as decreased conifer productivity, our simulations projected that an important proportion of Alberta’s boreal forests would transition to treeless habitat (i.e. grass-, or shrub-dominated vegetation) while many conifer-dominated stands would likely be replaced by broadleaf tree cover. Consequently, the abundance of bird species associated with open and deciduous habitats were projected to increase. With a strong anthropogenic climate forcing scenario (RCP 8.5), sharp declines in abundance of coniferous trees were also projected, particularly in mature and old forest stands, triggering major declines for bird species associated with coniferous and mixedwood forest types.

Main Conclusions

As the most comprehensive simulation of climate change and harvesting impacts on avian habitats in the North American boreal region to date, our study reveals the importance of considering key habitat characteristics like forest age structure and composition through forest landscape modeling, and identifies 18 bird species particularly sensitive to climate change. Our simulations suggest that a change in forest management practices could play an important role in the conservation of boreal bird species vulnerable to climate change. The intensive forest harvesting simulated accelerated declines in bird abundance compared to a “no harvesting” scenario.

We used LANDIS-II (Scheller & Mladenoff, 2004) to project future forest attributes within our study area. In LANDIS-II, the forest landscape is represented by a grid of interacting cells within which stand-level forest processes (tree establishment, growth, competition and mortality) occur while landscape-level processes, such as tree seed dispersal and forest disturbances including fire and harvesting, generally affect multiple cells in a spatially, interactive manner. In our experiment, we set cell resolution to 250 m (6.25 ha) and simulations were run at 10-year time steps across all activated extensions. Forest composition and structure in each cell were initialized using forest properties derived from the Alberta Biodiversity Monitoring Institute (ABMI) cover products (as of the year 2010; ABMI 2012) and the Canadian National Forest Inventory (NFI; nfi.nfis.org) and combined with stand age cohort data derived from provincial forest inventory plots. Each of these cells was then assigned to a ‘‘landtype’’ with homogeneous soil (Mansuy, Thiffault, Paré, Bernier, Guindon, Villemaire, et al., 2014) and climate conditions.

A modified version of the LANDIS-II Biomass Succession extension v 3.1 (Scheller et al., 2004) was used to simulate forest succession. The Biomass Succession extension emulates succession at the stand (cell) level by simulating the recruitment and growth of tree cohorts (not individual trees). It permits multiple cohorts of tree species to establish and interact within a cell through resource (i.e., “growing space” sensu Scheller et al. 2004) limitations based on species-specific traits. The succession of each cell is driven by these stand-level interactions, in addition to disturbance history and seed source availability. Specific parameters that define basic life-history traits are assigned to all species (see Table 1 for a full listing). To account for the effects of climate change, the forest gap model PICUS (version 1.5; http://picus.boku.ac.at) was used to develop the dynamic tree species- and landtype-specific parameters required to operate LANDIS-II. PICUS is an individual tree-based, spatially explicit forest ecosystem model that simulates germination, establishment, growth, and mortality of individual trees on 100 m² patches of forest area (see Boulanger, Taylor, et al. (2016) for more details). Hence, three dynamic inputs of the Biomass Succession extension, namely maximum biomass (maxAGB; g.m^-2), maximum aboveground net primary productivity (maxANPP; g.m^-2.yr^-1), and species establishment probability (SEP), were derived from PICUS by running monospecific stand simulations for each combination of species, climate conditions and landtype. Those parameters were allowed to change during the course of the subsequent LANDIS simulations to represent the effect of climate on each species’ potential growth. For a complete description of the calibration, the validation procedures regarding these parameters, and a description of how these dynamic inputs were derived from the outputs of PICUS, refer to Tremblay et al. (2018).

We considered two natural disturbance agents, namely wildfires and drought. Wildfire accounts for the majority of areas naturally disturbed in the study area (Tymstra, Wang, & Rogeau, 2005) and is widely recognized to have major impacts on Canada’s forest landscapes (Volney & Hirsch, 2005; Price et al., 2013). Fire simulations were carried out using the LANDIS-II Base Fire extension (He & Mladenoff, 1999), which simulates stochastic fire events dependent upon fire ignition, initiation and spread. Fire regime data (annual area burned, fire occurrence, and mean fire size) were first compiled into ‘‘fire regions’’ corresponding to the Canadian Homogeneous Fire Regime (HFR) zones (Boulanger et al. 2014). Baseline and future fire regime parameters within each fire region were calibrated with models developed by Boulanger et al. (2014) and they were updated to account for changing climate conditions under the different RCP scenarios (Gauthier et al., 2015).

Drought-induced mortality was simulated by first modelling species-specific mortality curves according to the climate moisture index (CMI), calculated as the difference between annual precipitation and potential evapotranspiration (see Brecka 2018 for more details). Species-specific mortality was retrieved from undisturbed permanent sample plots located within the Boreal Plains ecozone. This was used to construct non-linear exponential regression models predicting the 10-yr proportion of species biomass killed according to decadal CMI values (Chen, Luo, Reich, Searle, & Biswas, 2016). Using the same climate datasets described above, we projected future CMI values under all landtypes under each climate scenario and each 30-year period (i.e. 2011-2040, 2041-2070; 2071-2100). Future species-specific drought-related mortality was then projected using future CMI values and drought-related mortality models. Projected drought-related species-specific mortality was then included in the LANDIS-II simulations by removing biomass accordingly using the Biomass Harvest extension (v3.0; Gustafson, Shifley, Mladenoff, Nimerfro & He, 2000). Drought-induced mortality was applied equally to all tree age cohorts.         

Forest harvesting was simulated using the Biomass Harvest extension. Only clearcut harvesting was simulated as this logging strategy is most frequently used in the study area (https://alpac.ca/index.php/forest-sustainability/forest-planning). Only stands in upland areas and that comprised cohorts older than 60 years old were allowed to be harvested. When harvested, clearcutting was simulated to remove of all age cohorts present except for the 0 – 10 year age cohort. Mean harvested patch size and total harvested area were summarized by forest management units. Harvesting parameters were held constant throughout the simulations. Three harvesting scenarios were simulated according to a gradient of harvesting pressure, from no harvesting (no harvesting), to clearcutting with intensity similar to current management practices (baseline harvesting - applied to 0.3% of the harvestable upland area per year; ABMI, 2017a), to high-intensity clearcutting (high harvesting - applied to 0.6% of the harvestable upland area per year).

Simulations were run for three climate scenarios (baseline, RCP 4.5 and RCP 8.5) as well as under the three harvesting scenarios. Five replicate simulations were run for 200 years, starting in the year 2000, with 10-year time steps. Except for scenarios involving the baseline climate, fire regime parameters were allowed to change in 2010, 2040, and 2070 according to the average climate corresponding to each forcing scenario. Dynamic growth and establishment parameters (SEP, maxANPP and maxAGB) as well as drought mortality were allowed to change according to climate scenarios following the same schedule but only for upland areas. Indeed, our current understanding of the vulnerability of peatland systems to climate change is very limited (e.g., Schneider, Devito, Kettridge & Bayne, 2016). As such, lowland pixels were kept as “active” to allow fire spread and seed dispersal, but growth parameters were kept constant. A similar simulation strategy was used by Stralberg, Wang, et al. (2018) in this area. As a result, future forest landscape, as well as bird community results, were reported for uplands only.

Boreal songbird community

To represent the boreal bird community, we selected passerines with breeding ranges that overlapped with the study region. We further limited this selection to 72 songbirds that were adequately modeled within Northern Alberta by excluding species that were too rare (number of detections < 5 x degrees of freedom in models) or for which model goodness of fit was low (AUC < 0.6). These predictive models were based on point count data, including surveys from the North American Breeding Bird Survey (BBS; pwrc.usgs.gov/bbs/), Boreal Avian Modelling Project (BAM; borealbirds.ualberta.ca), and the ABMI (abmi.ca). The models were built following the methodology outlined in Ball, Sólymos, Schmiegelow, Haché, Schieck, & Bayne (2016) and Sólymos, Azeria, Huggard, Roy & Schieck (2019). Land cover associations were based on the dominant landcover (native vegetation and human footprint) type within a 150-m radius buffer around the points. Native vegetation classes included deciduous, mixedwood, white spruce, pine, black spruce forest stands, treed fen, shrub, grass/herb, graminoid fen, marsh, and swamp cover types. Ages of forest stands (area-weighted average age at the year of the survey) originating from natural disturbances or forest harvesting were also assessed within the 150-m radius buffers. Survey counts were modelled by Poisson generalized linear models with a logarithmic link. We used the QPAD approach (Sólymos, Matsuoka, Bayne, Lele, Fontaine, Cumming, Stralberg, Schmiegelow, & Song 2013) to account for differences in sampling protocol and covariate effects on detectability via offsets in the generalized linear models. This approach standardizes the estimates to reflect density (number of singing individuals per ha) within the different land cover type and stand age categories.

As a result of these models, ABMI provides expected density for the selected bird species for each cover type, and data are available via the `cure4insect` R extension package (R Core Team 2019, Sólymos, Allen, Azeria, White, ABMI & BAM 2018; see model summaries at abmi.ca/data). ABMI bird habitat models were built only for species that had enough samples and were adjusted for the species detection distance (see ABMI (2017b) for modelling details). We used expected density per species in conjunction with projected forest cover types to estimate expected bird abundance. Aboveground biomass density (t/ha) for each tree species as well as stand age and stand origin as projected by LANDIS-II were used to derive habitat types (forest cover types) using the same classification scheme used to define ABMI forest cover types (see above). Stand-scale (250 m) forest cover information was derived for each LANDIS-II simulation run at a 10-yr time step. We also associated each bird species to the following habitat types (deciduous, mixedwood, coniferous and treeless) and stand age classes (<30yr (young); 30-60yr (closed); 60-80yr (mature); >80yr (old)).

Analyses

Cumulative impacts of harvesting and climate change were assessed by comparing temporal trends of tree species aboveground biomass and songbird abundance under each climate and harvesting scenario over time. Outputs from the five simulation replicates were averaged. Trends were assessed using simulations where all disturbances (fire, drought and harvesting) were considered. The impact of climate change and harvesting on simulated songbird abundance was calculated as the percentage of change in simulated songbird abundance relative to the proportion obtained under the baseline harvesting and baseline climate scenario (hereafter referred to the “reference scenario”) following:

((ProjAbund_t / RefAbund_t) - 1)*100                                         [1]

where RefAbund_t is the abundance of a bird species under the reference scenario (baseline climate and baseline harvesting), ProjAbund_t is the projected abundance of the same species for the given future time period, and t is time in years (Cadieux et al., 2019). This method was used to get a direct assessment of the effects of climate change and harvesting while controlling for forest succession.

Additional simulations were conducted to assess the importance of a selected driver of change on bird abundance. Simulations were conducted following a three-way factorial design according to harvesting (no harvesting, 0.3% and 0.6%), fire (baseline fire, projected fire) and climate change effects on dynamic biomass succession inputs and drought. The relative contribution of each factor was assessed by estimating the variance of songbird abundance that it explained using omega-squared values (ω²) calculated following a 3-way factorial ANOVA. We calculated ω² for each driver of change, at each time step, as:

ω² = [SS_effect – (df_effect)*(MS_error)] / [MS_error + SS_tot]                             [2]

where SS_effect is the sum of squares related to the driver of change (the effect), df_effect is the degree of freedom of the effect, MS_error is the mean square of the error and SS_tot is the total sum of squares. ANOVA and ω² calculations were performed separately for each RCP scenario. Based on the relative importance of drivers on bird species, we defined species sensitive to climate change as those declining > 25% according to climate-sensitive drivers (either Fire or Growth) under RCP 8.5 and the current harvest scenario at 2100.