This study leveraged bird census and forest structure data of the Canton of Zurich, Switzerland (1,729 km²), located in the eastern part of the Swiss plateau (Figure 1). The climate is continental with strong influence from the Atlantic, having cold winters and warm summers, a mean annual temperature of 9.9°C, and mean annual precipitation of 1,179 mm. Forests cover ca. 31% of the surface and span from deciduous to subalpine conifer forests. Coniferous forests are the most common in the region, with dominance of Norway Spruce (Picea abies; 38%) and Silver Fir (Abies alba; 12%); the most common broadleaved tree species is European Beech (Fagus sylvatica), with a coverage of 23% (Amt für Landschaft und Natur, 2020).

The forests in the Canton of Zurich are representative for forests of the central European continental biogeographical region (Sabatini et al., 2021), which dominates Central Europe and the plains and valleys surrounding the Alps. These forests are mostly dominated by the abovementioned species, namely Norway spruce, European Beech and Silver Fir, as well as Scots Pine (Pinus sylvestris) and oak (Quercus sp.; Leuschner and Ellenberg, 2017).

We used data from the breeding bird atlas of the Canton of Zurich 2008 (Weggler et al., 2009), with a total of 1,440 habitat polygons (with sizes ranging from 6.4 to 80.2 ha; mean = 49.24 ± 0.30). We only used polygons, covering both forested and mixed use polygons that included forests. These data were collected from 2006 to 2008, using transect counts conducted by 250 trained volunteer ornithologists. Each polygon was surveyed 5 times between the 20th of March and 30th of June, for a period of 30 min, in the early hours of the day (from sunrise until 10 am), with every bird seen or heard recorded. We analyzed data from 15 forest bird species (Table 1; Keller et al., 2010), covering different habitat preferences and ranging from rare to common in the study region (prevalence values between 0.024 and 0.915). To reduce uncertainty in the data we reduced the original relative abundance information into presence and absence (abundance ≥ 0 converted to presence).

We chose predictors for our ecological niche models relating to forest stand structure, landscape structure, and climate, likely of relevance for the 15 bird species, as described below.

Two sources of forest inventory data were available for this study: the national forest inventory of Switzerland and the forest inventory of the Canton of Zurich (Brassel and Lischke, 2001). The sample densities of the national inventories are 1 × 1 km for the first inventory (Brändli, 1996), and 1.4 × 1.4 km for the second and third inventories (Brassel and Lischke, 2001; Brändli, 2010), using a system of concentric plots for different diameter classes (200 m² for diameter at breast height (dbh) < 12 cm, 500 m² for dbh < 36 cm and an auxiliary plot for ingrowth). The data from the Canton of Zurich has a denser grid, with circular inventory plots at each 80 × 300 m and a fixed sized of 300 m². The main data set used for the simulation was the one from the Canton of Zurich, due to the denser grid, while we used the national forest inventory data to calibrate the dbh-height relationship for not measured heights in the Canton’s data.

We selected forest variables that are likely of ecological relevance for the focal species, and possible to produce during simulations of future forest dynamics and management with the forest simulator SILVA (Pretzsch et al., 2002; Pretzsch et al., 2008). The selected variables (Table 2) were calculated for three scales. These are, the entire stand, each tree species group (broadleaves and conifers), and each of the (most common) individual tree species that are implemented in SILVA: Norway Spruce, Silver Fir, Scots Pine, European Beech, Common and Sessile Oak (Quercus robur and Q. petrea, considered together due to their structural similarities), as well as European Larch (Larix decidua). Additionally, the North American Douglas Fir (Pseudotsuga menziesii), although not very common in the study area, was considered for the future scenarios due to the increasing interest in its use in Europe (Isaac-Renton et al., 2014). The mean value of each predictor among the inventory plots within each forest stand was assigned to the respective stand (Amt für Landschaft und Natur, 2010). As forest stands (delineated on the forest map described below) were considerably smaller than the bird sampling polygons, the variables needed to be compiled to the bird polygon level by calculating the area-weighted mean value of all contained forest stands (for all variables), and the respective standard deviation as a measure of heterogeneity (for the stand level variables only). Tree density is the total number of trees divided by the area of the inventory plot; crown coverage is the ratio between the space occupied by crowns and the plot area; tree volume is calculated as a function of the tree’s basal area, tree height, and species-specific form factors; and tree diversity is the Shannon diversity index calculated from the tree density values for each tree species.

The landscape structure variables were calculated based on forest stand data (Table 2). Percentage forest cover and mean forest patch size were directly derived from a forest map, compiled after an aerial photographic survey, and obtained from the Amt für Landschaft und Natur (2010). The remaining three variables were calculated using the Guidos Toolbox (Soille and Vogt, 2009; Vogt and Riitters, 2017), considering a forest edge width of 50 m. Forest islets are small, isolated “islands” of forests within the open landscape. The four climatic variables (Table 2) chosen were Annual mean temperature, Temperature seasonality (annual range in temperature), Mean temperature of wettest quarter and Precipitation seasonality (annual range in precipitation), and were extracted from the CHELSA “bioclim” variable dataset (Karger et al., 2017).

Forest simulations were performed using the forest simulator SILVA (Pretzsch et al., 2002, 2008) which was developed to support practitioners of sustainable forest management. SILVA is a single-tree-based model that is distance-dependent (tree positions matter) and age-independent. The time scale to be simulated ranges from 5 years up to a rotation period. SILVA is optimized for the most important tree species in central Europe in pure and mixed stands. With local adjustment, e.g., based on inventory data, it can be tuned for about 80% of the central European forests.

The core equation for estimating the site dependent height growth potential is a Chapman-Richards function:

which describes a tree’s potential height h_pot at a given age t dependent on three site-dependent (and species specific) parameters a, k, and p (Kahn, 1994; Pretzsch et al., 2002). Using height-age relationships from long term research data from experimental plots located in the Canton of Zurich (data provided by the Swiss Federal Institute for Forest, Snow and Landscape Research -WSL), new parameters (a, k, and p) were fitted for each tree species and groups abovementioned and the corresponding curves were implemented in SILVA. The long-term research offered the advantage of covering a broad age-height range for all species used in SILVA. The 90% quantile was used to fit the potential growth curves to avoid biases caused by thinning experiments.

SILVA simulates growth of single trees in forested conditions. As the competition among trees is evaluated in a spatially explicit way, the model can cover a broad range of existing, and also novel, silvicultural concepts. SILVA can simulate for even-aged or uneven-aged mixed and monospecific forests. The model is generally applied stand-wise, but, as it is the case in this study, it can also be applied on a landscape level where grid-based forest inventory data are available.

For model initialization, we grouped all structure data from the inventory plots into representative strata based on the main dominant species, diameter distribution, top height and land tenure type. Each stratum was simulated independently to improve computing performance. Once the simulations were ready, the structural information was transferred to each inventory point. In a second step, the structure variables were aggregated for each of homogenous forest patch (∼0.6 ha) as described in section “Environmental Data.”

We investigated the difference in predicted bird occurrence patterns between future potential and realistic forest management scenarios. We selected five different scenarios which cover the main silvicultural regimes applied in Central Europe. Each scenario is described by harvesting intensities and the main types of silviculture. Thus, we can simulate silvicultural regimes ranging from wood production oriented, involving high harvest intensities, to more conservation-oriented silvicultural regimes with a higher emphasis on biodiversity (close-to-nature) or protection (set-asides). The scenarios are (see below for their description): Continuous cover; Multifunctional; Wood production; Habitat tree; and Set-aside (Table 3). Each scenario was simulated for 100 years, in 5-year intervals.

The Continuous cover scenario simulates an intensified version of continuous cover regimes that are widely applied by Central European forest offices. Productive monocultures, including Norway Spruce and Scots Pine stands, are managed using selective thinning including the future crop tree thinning concept (Bücking et al., 2007), while Beech stands are managed according to the shelter-wood concept, with different intensities depending on site quality. Thus, the main objective of this scenario is the sustainable wood production with a focus on forest stability and resilience. As silviculture in the Canton of Zurich is similar to the one applied in Bavaria, we decided to use prescriptions that are already included in SILVA for the Bavarian State Office.

The Multifunctional scenario is designed to balance different forest ecosystem functions and services simultaneously. These include wood production together with groundwater recharge, carbon sequestration and biodiversity. Thus, the predominant objectives of this scenario include a continuous crown cover and an active promotion of high structural heterogeneity and tree species diversity, including a high share of deciduous species. Specifically, silviculture is applied as follows: up to a tree height of 12 m, the stand is simulated through a stem number reduction that must not exceed a total volume of 15 m³/ha. Selective thinning is conducted at tree heights between 12 and 32 m, to remove up to 55 m³/ha in deciduous stands and up to 70 m³/ha in coniferous stands. The target diameter harvest phase starts at tree heights of 32 m and removed between 80 and 144 m³/ha of conifers and 70 m³/ha in deciduous-dominated stands. All the conducted treatments are applied in turns of 10 years (Toraño Caicoya et al., 2018). In this scenario, natural regeneration is the main contributor to regrowth. To improve species diversity, Scots Pine (500 trees/ha) is planted during the regeneration phase in the conifer-dominated stands. In the deciduous-dominated stands, additional to the natural regeneration (Poschenrieder et al., 2018), Scots Pine (250 trees/ha) and Douglas Fir (250 trees/ha) are planted during the regeneration phase to increase the stand’s multifunctionality, and additionally European Beech (6,000 trees/ha) is planted to ensure that the desired species mix is achieved.

The Wood production scenario focuses on the production of wood, where the volume of harvested wood is the first priority. It assumes an increase in the demand of wood products in Europe within the next 100 years. In this scenario, stem reduction is applied on deciduous tree species, removing up to 25 m³/ha, while no stem reduction is applied to conifer-dominated stands. In the deciduous stands, a selective thinning is applied at a tree height of 12–17 m, and thinning from below for a tree height of 17–30 m, removing up to a total of 25 m³/ha. The final felling is conducted by a light selective thinning and a target diameter felling of trees with a dbh diameter between 20 and 200 cm and a removal of up to 500 m³ of the standing volume. In the conifer stands, strong selective thinning and target diameter felling is applied in two height phases. During the first phase from 12 to 19 m, up to 60 m³/ha are removed, targeting conifers with a dbh diameter between 40 and 200 cm and broadleaves with diameters > 5 cm. During the second phase, up to 120 m³/ha are removed for diameters between 50 and 200 cm for all species. The final felling phase aims at trees with a height > 31 m, applying a light selective thinning and a target diameter felling with a removed standing volume of up to 500 m³/ha and a diameter between 5 and 200 cm for deciduous trees and a diameter between 20 and 200 cm for conifers. In addition to the natural regeneration (Poschenrieder et al., 2018), Norway Spruce (4,000 trees/ha) and Douglas Fir (100 trees/ha) are planted during the regeneration phase.

The Habitat tree scenario is a variation of the multifunctional scenario. This scenario aims at increasing the number of large trees remaining in the stand. For all tree species, residual large trees are promoted, however, with a preference for broadleaves and Oak in particular. To achieve this, the regeneration phase is modified as follows: for stands with a broadleaved species as dominant species, the maximum diameter during the target diameter felling is reduced for any possible admixture of oak and broadleaves until 60 cm, and leaving 20% of the volume for the oak and 15% for the rest of the broadleaved species standing. The intensity of the selective thinning for oak (being a light demanding species) is intensified to favor its growth. For the stands with a conifer species as dominant species, the multifunctional scenario remains unaltered for the conifer admixtures, but no treatment is applied to any broadleaved admixture in the regeneration phase. In summary, in broadleaved stands around 15–20% of the broadleaves are allowed to grow until they reach their natural mortality, favoring trees with dbh diameter > 15 cm, and in conifer stands all broadleaves in the admixture are left without management.

The Set-aside scenario, in which no management was conducted and the forest dynamics was driven solely by natural mortality and regeneration, was used as the reference scenario. The bird responses in the other scenarios were compared to the bird responses the Set-aside scenario.

The bird occurrence data were fitted to the environmental data using Boosted Regression Trees (BRTs; Friedman et al., 2000; Elith et al., 2008). BRTs combine decision trees and boosting in a single algorithm. While decision trees relate a response to their predictors through recursive binary splits, boosting combines many simple models to give improved predictive performance. This way, BRTs work as an additive regression model where many individual decision tree models are fitted in a forward stagewise fashion (Elith et al., 2008). Thus, they are very flexible and capable of dealing with complex (non-linear) responses, and are robust to missing data. Also, they include a stochastic element in the boosting process, which by randomly selecting subsamples (bag fraction) to fit the data at each iteration, makes them robust against overfitting (Friedman, 2002). For this reason, BRTs are highly performing and are a primary choice for ecological niche (or habitat) modeling (e.g., Elith et al., 2006).

We used BRTs as implemented in R version 3.6.2 (R Development Core Team, 2020), using code from the package “gbm” version 3.1.5 (Ridgeway, 2019). We ran the BRTs using the recommended bag fraction of 0.5 (Elith et al., 2008). The BRT parameterization is done by defining three main parameters: the maximum number of splits in each decision tree (tree complexity); the number of trees to use (model complexity); and the model regularization (or shrinkage) parameter (Elith et al., 2008). The ideal parameter set depends on the particular data used, so we first defined the number of trees to be used for each case via cross-validation by maximizing the model log likelihood. The tree complexity and shrinkage parameters were defined via a heuristic search of every possible pair combination, again in a cross-validation procedure, by optimizing each model’s area under the receiver operating characteristics curve (AUC; Swets, 1988).

Initially, we fitted one model per species using all 23 predictors considered. In a second step, we simplified the models, based on the initial models’ variable contributions, and by ensuring that highly correlated variables (r ≥ 0.7) were excluded, thus avoiding further modeling problems (Dormann et al., 2013). In the final (reduced) models, we included six forest structural, two landscape structural, and two climatic variables per model, and in this way gave emphasis to forest structural variables (whose effect is investigated), while acknowledging the role of landscape and climate (Table 4). The final BRT models were reiterated using Generalized Linear Models (GLMs; Hosmer and Lemeshow, 2000) to assess consistency on the species responses (when using a more parsimonious modeling approach). The GLMs were developed by testing to include linear and quadratic terms for all of the ten predictors included in the final models for each species, and then running a multi-model selection (Burnham and Anderson, 2004) based on the Akaike Information Criterion (AIC; Akaike, 1974). This analysis was performed in the R package “MuMIn” version 1.43.17 (Barton, 2020).

All models were checked for ecological plausibility, and evaluated using a 10-fold cross-validation procedure, optimizing AUC. The BRT models were used to interpret the species responses to the environment, and all models with an AUC value ≥ 0.7 were used to project the species occurrences onto the future forest management scenarios. Both the BRTs and GLMs are described according to the ODMAP protocol for reporting species distribution models (Zurell et al., 2020; see Supplementary Information 1).

The BRT models achieved cross-validated AUC values between 0.706 and 0.914. Forest structural parameters were generally highly influential, contributing the most for eight species models, and the second most for the remaining models (Table 5). The GLM model performances ranged from AUC values of 0.622–0.918, while seven of the 15 models did not achieve a sufficient AUC value (0.7) to be further used for prediction or inference.

Within the forest structural variables, the most commonly selected variable (in the BRT models) was the variation (standard deviation) in stand volume (selected in N = 12 models), followed by (mean) volumes of spruce (N = 10 models), broadleaved species (N = 9 models) and oak (N = 8 models). Several bird species showed clear responses to forest structural parameters, including the Middle Spotted Woodpecker and the Golden Oriole (associated with areas with high oak volumes), the Lesser Spotted Woodpecker (high broadleaved and oak volumes), the Wood Warbler (high oak and larch volumes), the Common Firecrest (high spruce volumes), the Crested Tit (low broadleaved volumes), the Spotted Nutcracker (selecting areas with high crown cover) and the Hawfinch (associated with areas with intermediate variation in tree density). We note that in the seven GLM models with acceptable results (AUC higher or equal to 0.7), the species responses identified were qualitatively similar to those of the BRT models—see Supplementary Informations 2, 3 for further details on BRT variable importance and the individual species responses to the environmental variables (partial dependency plots and variable contributions).

The distinct trajectories of forest structural changes among the scenarios are visualized via the development of the standing volume (Figure 2, left) and the tree species diversity (Shannon index; Figure 2, right). The largest final standing volume was achieved by the Set-aside scenario (∼770 m³/ha), followed by Continuous cover (∼600 m³/ha). Multifunctional and Habitat tree followed similar trends, with a small reduction after the first simulation period and a constant increase afterward. The Wood production scenario showed the highest cyclic fluctuations due to the dominance of rotation forestry. In contrast, final tree species diversity was highest in the Set-aside scenario, closely followed by Continuous cover. Both progressively recovered from a drastic decrease after the first simulation periods. Multifunctional and Habitat followed again similar trends, showing the most stable trajectories. The largest changes were observed in the Wood production scenario, showing the trends of rotation forestry, with drastic decreases of diversity after each final harvest.

Bird species showed variable responses to the different management scenarios (Figure 3). Coal Tit, Common Chiffchaff and Goldcrest had almost no change in their mean probabilities of occurrence across the study area and the time period of analysis. In contrast, Black Woodpecker, Great Spotted Woodpecker, Lesser Spotted Woodpecker, Common Firecrest, and Short-toed Treecreeper showed a consistent decrease in their occurrence probabilities in relation to the starting point. The occurrence of Golden Oriole seemed to benefit from short-term management actions (first 25-year time period) of most scenarios (apart from the Wood production scenario), although their mean occurrence probabilities in all scenarios decreased in the intermediate and long-term (50–100 years). Middle Spotted Woodpecker, Crested Tit, Eurasian Treecreeper, and Hawfinch had variable responses depending on the management regimes. While Middle Spotted Woodpecker benefited particularly from the scenarios that promoted the establishment of large oak trees within the forest (i.e., Habitat tree and Multifunctional), Crested Tit and the Eurasian Treecreeper rather benefited from the Wood production scenario. Hawfinch, while also benefitting from the two scenarios that promote oak trees, had a fluctuating response to the Wood production management scenario. Wood Warbler and Spotted Nutcracker generally benefited from all analyzed forest management scenarios, although to a different degree. While the vulnerable Wood Warbler greatly benefitted from the oak tree promoting scenarios, the Spotted Nutcracker benefitted from the Continuous cover and the Set-aside scenarios.