Developing fine-grained nationwide predictions of valuable forests using biodiversity indicator bird species.

The use of indicator species in forest conservation and management planning can facilitate enhanced preservation of biodiversity from the negative effects of forestry and other uses of land. However, this requires detailed and spatially comprehensive knowledge of the habitat preferences and distributions of selected focal indicator species. Unfortunately, due to limited resources for field surveys, only a small proportion of the occurrences of focal species is usually known. This shortcoming can be circumvented by using modeling techniques to predict the spatial distribution of suitable sites for the target species. Airborne laser scanning (ALS) and other remote sensing (RS) techniques have the potential to provide useful environmental data covering systematically large areas for these purposes. Here, we focused on six bird of prey and woodpecker species known to be good indicators of boreal forest biodiversity values. We used known nest sites of the six indicator species based on nestling ringing records. Thus, the most suitable nesting sites of these species provide important information for biodiversity-friendly forest management and conservation planning. We developed fine-grained, that is, 96 × 96 m grid cell resolution, predictive maps across the whole of Finland of the suitable nesting habitats based on ALS and other RS data and spatial information on the distribution of important forest stands for the six studied biodiversity indicator bird species based on nesting-habitat suitability modeling, that is, the MaxEnt model. Habitat preferences of the study species, as determined by MaxEnt, were in line with the previous knowledge of species-habitat relations. The proportion of suitable habitats of these species in protected areas (PAs) was considerable, but our analysis also revealed many potentially high-quality forest stands outside PAs. However, many of these sites are increasingly threatened by logging because of increased pressures for using forests for bioeconomy and forest industry based on National Forest Strategy. Predicting habitat suitability based on information on the nest sites of indicator species provides a new tool for systematic conservation planning over large areas in boreal forests in Europe, and a corresponding approach would also be feasible and recommendable elsewhere where similar data are available.

INTRODUCTION

Land‐use and land‐use change have been the most significant negative drivers of natural terrestrial ecosystems during the last 50 years (Intergovernmental Science—Policy Platform on Biodiversity and Ecosystem Services; IPBES, 2019), causing detrimental habitat loss and degradation, and leading to species population declines (IPBES, 2019; Sala et al., 2000). Land‐use–based pressures on biodiversity have further accelerated in the early 21st century (Butchart et al., 2010). In addition to climate change, land‐use change is thus predicted to have a very large negative effect on biodiversity during this century as well (Beyer & Manica, 2020; Sala et al., 2000). Within boreal forests ecosystems, forestry, that is, logging and management of forests to the extraction of wood, is among the main threats for the persistence of forest‐dwelling species (Gauthier et al., 2015).

Logging and management result in the fragmentation of contiguous forest landscapes, structurally homogeneous even‐aged stands, and a reduced volume of dead wood (e.g., Kuuluvainen, 2009). In Europe, logging volumes have increased, and the intensity of logging activities has been at the highest level in northern Europe, particularly during the last 10 years (Ceccherini et al., 2020). Large‐scale logging has generated a considerable loss of old‐growth forests in northern Europe during the last 60 years, causing clear negative impacts on forest biodiversity (Fraixedas et al., 2015; Roberge et al., 2018; Svensson et al., 2019). In Finland and Sweden, almost all forest land outside protected areas (PAs) is subject to systematic silvicultural practices, which include thinning of young stands, clear‐cutting at a stand age of 60–120 years, and then replanting to start the rotation cycle (Nordberg et al., 2013; Virkkala & Toivonen, 1999).

In mitigating the negative effects of forestry and other land use on forest biodiversity, sparsely occurring PAs with limited spatial cover, such as those in southern and central Finland, are not alone sufficient to halt the decline of the species suffering the most from the loss of old‐growth and contiguous forests (Esseen et al., 1997; Roberge et al., 2018; Virkkala et al., 2020). This suggests that new PAs should be established, and more biodiversity‐friendly measures should be increasingly applied in forest management in the unprotected parts of forest landscapes to enhance forest species conservation. Importantly, such actions should be targeted at the forest stands with notable potential conservation values. However, a typical obstacle for this is the shortage of data that both cover systematically large areas and provide local‐scale information on occurrences of multiple rare or declining forest species, as well as the critical habitat factors that they depend on.

In this work, we apply an approach based on the concept of biodiversity indicator species, local forest and land‐use data, and species distribution modeling (SDM) tools, to develop nationwide fine‐grained predictions of forest stands with high conservation values. Using this approach, we aimed to develop important spatial information for forest conservation and management planning at different spatial scales, which is crucially needed to preserve species from the negative effects of intensive forestry (Angelstam et al., 2004; Crosby et al., 2020). With respect to SDMs, numerous different approaches (e.g., Elith et al., 2006) are available and increasingly employed to develop regional or national maps of suitable sites for a given species when only a small proportion of all occurrences of a species is known (Brambilla et al., 2015; Zhang et al., 2019).

There are certain preconditions for our approach to be successful. First, there should be sufficient a priori knowledge of the habitat preferences and distributions of focal indicator species, upon which robust SDMs can be built to develop predictions of the most suitable sites for the species. Second, to be able to produce SDMs across wide areas, we need systematically gathered fine‐grained nationwide data that can provide ecologically relevant predictors for the focal species. For boreal forests, airborne laser scanning (ALS) and other remote sensing (RS) techniques are among the most promising means to construct such environmental data, enabling successful SDMs (Vihervaara et al., 2015). Third, it is imperative that the focal indicator species are both sensitive to the impacts of forestry and are also indicative of potential forest biodiversity values across a wider range of species and taxonomic groups (Burgas et al., 2014; Roberge, Mikusiński, et al., 2008). When these preconditions are met, the derived predictive maps can be used to inform not only local forest management planning, but also to provide information across larger regions and entire countries, thereby allowing detection of sites of national importance (Crosby et al., 2020; Lehtomäki et al., 2009).

Nationwide species occurrence maps based on the habitat suitability model have been produced in several studies, in boreal regions recently for the Siberian flying squirrel Pteromys volans (Jokinen et al., 2019), the Eurasian beaver Castor fiber, and the invasive North American beaver Castor canadensis in Finland (Alakoski et al., 2020). In Sweden, citizen‐science observations were used to predict habitat suitability for the Siberian jay Perisoreus infaustus, a declining species preferring old‐growth forests (Bradter et al., 2018) and potential areas were identified for the expanding wolf Canis lupus population (Eriksson & Dalerum, 2018).

In boreal forests, a number of earlier studies have shown that bird species preferring mature and old‐growth forests have notably declined in Finland in the era of intensified forestry and increased logging (e.g., Fraixedas et al., 2015; Kouba et al., 2020; Virkkala et al., 2020). Moreover, previous studies have shown that certain birds of prey and woodpecker species are useful indicators of overall biodiversity values in forests (Burgas et al., 2014; Martikainen et al., 1998; Roberge, Mikusiński, et al., 2008; Roberge & Angelstam, 2006), this indicator status being based on a preference for mature or old‐growth forest conditions by these species or their ecological impacts, such as woodpeckers providing cavities for secondary cavity‐nesters (Martin & Eadie, 1999; Virkkala, 2006). Here, we develop SDMs for six indicator hawk and woodpecker species based on data on their nest sites to derive maps of sites most suitable for nesting, as such predictions can provide particularly useful indications of important locations for biodiversity‐friendly forest management and conservation planning (Björklund et al., 2020; Burgas et al., 2014; Martikainen et al., 1998; Pakkala et al., 2014; Roberge, Mikusiński, et al., 2008; Roberge & Angelstam, 2006; Virkkala et al., 1993).

Our ultimate aim and primary objective was to apply our modeling approach systematically across the whole country, including all forested areas in Finland, and thereby develop fine‐grained national (Figure 1) information and predictive maps of the locations of suitable nesting forests and important areas for our six biodiversity indicator bird species. Based on their biodiversity indicator role, the predicted most suitable nest sites of these bird species provide important information for conservation and management planning at a national scale. As a secondary objective, we examined how large a proportion of the predicted most suitable areas for our indicator bird species is situated in PAs, and how much is outside the PA network, in locations that might provide important sites to complement the PA network and enhance the preservation of boreal biodiversity.

FIGURE 1

Location of the study area (Finland) in northern Europe, and protected area network in Finland

Our indicator species included three hawk species, the European honey buzzard (Pernis apivorus), the northern goshawk (Accipiter gentilis), and the common buzzard (Buteo buteo), and three woodpecker species, the white‐backed woodpecker (Dendrocopos leucotos), the lesser spotted woodpecker (Dryobates minor) and the Eurasian three‐toed woodpecker (Picoides tridactylus) (see Björklund et al., 2015; Burgas et al., 2014; Martikainen et al., 1998; Pakkala et al., 2014; Roberge, Mikusiński, et al., 2008; Roberge & Angelstam, 2006; Virkkala et al., 1993). All three hawk species and the white‐backed woodpecker are regarded as red‐listed (threatened or near‐threatened; Hyvärinen et al., 2019), in addition to which the lesser spotted woodpecker and the Eurasian three‐toed woodpecker have shown considerable long‐term declines in Finland (Väisänen et al., 1998).

MATERIAL AND METHODS

Species data

The locations of the nest sites of the focal species were extracted from the bird ringing data gathered and maintained by the Finnish Natural History Museum (LUOMUS). In our study, we used ringing records from the years 1990–2017 from all nest sites in Finland. However, the ringing data we rather small for the lesser spotted woodpecker. Thus, for this species, we supplemented the ringing data with confirmed nesting observations from the Third Finnish Breeding Bird Atlas carried out in 2006–2010 (Valkama et al., 2011). Only nesting sites that were georeferenced with horizontal accuracy of ≤100 m were included in the study.

In order to avoid erroneously linking nesting data to forest locations that no longer contain suitable mature forest, we excluded nest sites that had been substantially logged after the last ringing event. More precisely, we excluded the nesting site if 30% or more of the forest in the 96 × 96 m square delimited around the nest was logged, as the indicator bird species only breed in mature or old‐growth forests. The logged sites were detected from the database maintained by the Finnish Forest Centre (FFC), complemented by the yearly forest loss information downloaded from the Global Forest Watch database (Hansen et al., 2013). Because these logging announcements have been systematically available only after 2000 and yearly forest loss in the Global Forest Watch database has been recorded from 2001 onwards, we also excluded forest stands under 30 years old. After these exclusions of logged nest sites, the number of separate breeding sites used in modeling was as follows: the European honey buzzard, n = 416; the northern goshawk, n = 4,219; the common buzzard, n = 2,070; the white‐backed woodpecker, n = 182; the lesser spotted woodpecker, n = 123; and the Eurasian three‐toed woodpecker, n = 204. These sites were assigned into one of the grid cells in the 96 × 96 m lattice system covering the whole study area, based on site coordinates.

Forest, land‐cover/land‐use, and climate data

For each grid cell in the 96 × 96 m lattice system, we measured environmental predictor variables representing four variable groups: (1) forest‐stand characteristics, (2) land cover/land use within the 96‐m grid cells, (3) land cover/land use at the broader landscape scale, and (4) two variables showing the geographic variation in summer and winter thermal conditions.

Here, data for forest‐stand characteristics were compiled from three national sources, FFC, Metsähallitus Parks & Wildlife (MPW) and the multi‐source national forest inventory carried out by the Natural Resources Institute Finland (LUKE). We applied these data in two main steps. In the first step, we used a 16 × 16 m lattice system covering the whole of Finland, providing the best possible basis for measuring fine‐resolution forest‐stand characteristics. When constructing the 16‐m resolution forest data, the three national data sources were used in a hierarchical order; that is, for each 16‐m grid cell we used the data that were considered qualitatively the most accurate, or that had the best spatial coverage. Where this was not possible, the next most accurate data were used. The priority order of the forest data sources was as follows: (1) Data on forests collected by FFC, based on ALS (Pukkala, 2019), covered 53% of forests. For the privately owned commercial areas not yet covered by ALS data, we used forest‐stand and compartment data provided by FFC (20% of forest area); (2) data on state‐owned and private PAs not available to ALS and collected by MPW, including detailed forest‐stand and compartment data recorded from the PAs (7% of forest area; Figure 1); and (3) data from multisource national forest inventory developed by LUKE for the whole of Finland based on field survey site data and satellite images (Tomppo et al., 2008), used here for the areas not covered by the previous data sets (20% of forest area).

In the second step of forest data processing, we aggregated variables describing the ecological characteristics of forests stands from the 16‐m resolution data used in the original databases aggregated (6×) to a coarser 96 × 96 m lattice system. Values of the following variables were calculated as the mean of the corresponding values in the 36 16‐m pixels constituting the 96 × 96 m cell: (1) forest‐stand volume (m3/ha); (2) mean trunk diameter at breast height (DBH, cm); (3) mean basal area (m2/ha); (4) tree height (m); (5) dominant tree age (years); and (6) volume of deciduous trees (m3/ha) (see Figure 2, Table 1). The coarser 96‐m data were then used in the modeling of bird nest sites. We focused on these six variables because the data on these variables systematically and with sufficient quality covered the whole country and because these variables provide the most relevant variables ecologically for modeling birds among the ALS data‐based variables (see Vihervaara et al., 2015).

FIGURE 2

Values of total stand volume (a, m3/ha), stand age (b, years), deciduous stand volume (c, m3/ha), and distribution of forest site types (d, 1–4, see section Forest, land cover/land use and climate data in Material and methods) in Finland

TABLE 1

Variables and data used in the analyses in the stand (96 × 96 m) and landscape scale (500 m or 1 km) and original resolution of the different data sets, which all were aggregated to 96‐m resolution in the analyses

Variable	Stand scale	Landscape scale	Data used	Original resolution
Stand volume	+	+	ALS, RS, forest‐stand data	16 m
Stand dbh	+	−	ALS, RS, forest‐stand data	16 m
Basal area	+	+	ALS, RS, forest‐stand data	16 m
Stand height	+	−	ALS, RS, forest‐stand data	16 m
Dominant tree age	+	−	ALS, RS, forest‐stand data	16 m
Deciduous tree volume	+	−	ALS, RS, forest‐stand data	16 m
Main tree species	+	+	ALS, RS, forest‐stand data	16 m
Forest site type	+	−	ALS, RS, forest‐stand data	16 m
Shoreline forest	+	+	CORINE land cover	20 m
Open, treeless mires	+	+	CORINE land cover	20 m
Transitional woodland/shrub on peatland	+	+	CORINE land cover	20 m
Forest on peatland	+	+	CORINE land cover	20 m
Marshlands	+	+	CORINE land cover	20 m
Agricultural areas	+	+	CORINE land cover	20 m
Urban areas	+	+	CORINE land cover	20 m
Water areas	+	+	CORINE land cover	20 m
January mean temperature	+	−	Mean 1981–2010	50 m
GDD	+	−	Mean 1981–2010	50 m

Abbreviations: +, data used; −, data not used; ALS, airborne laser scanning; dbh, diameter breast height; GDD, growing degree days; RS, other remote sensing (satellite).

Two additional variables were measured for the forest‐stand characteristics: (7) dominant tree species (Scots pine [Pinus sylvestris], Norway spruce [Picea abies], or deciduous trees), as recorded across the 16‐m pixels within a 96‐m cell; and (8) dominant forest site type, measured as the median value of the site types recorded for the 16‐m cells using four broadly defined site type classes (Table 1). The four site type classes were ordered according to their productivity as (1) herb‐rich forests and low‐herb heath forests (Oxalis–Myrtillus forest site type), (2) mesic heath forests (Myrtillus forest site type), (3) submesic heath forests (Vaccinium forest site type), and (4) dry heath forests (Calluna forest site type), lichen‐rich pine forests (Cladina forest site type) and rocky areas (Figure 2).

The second group of predictor variables included eight land‐cover/land‐use variables. All these variables were measured from the CORINE Land Cover 2018 database available for the whole of Finland at a resolution of 20 m. The eight CORINE land‐cover/land‐use types were the following: (1) shoreline forests (i.e., 20 × 20 m forest pixels along the shoreline with crown cover >30%); (2) open, treeless mires; (3) transitional woodland/shrub on peatland, crown cover 10%–30% (i.e., sparsely wooded pine mire); (4) forest on peatland, crown cover >30% (i.e., wooded pine mire or spruce mire); (5) marshlands; (6) agricultural areas; (7) urban areas; and (8) water areas (Table 1). Similarly to the case of the forest structure variables, these land‐cover variables were selected based on a priori knowledge of the local and landscape‐level habitat preferences of the six study species, that is, tailoring the selection of CORINE variables to match with the ecological requirements of the species, thereby facilitating the construction of robust SDMs.

The values for these variables were first converted from the 20‐m grid data to the 16‐m lattice system using nearest‐neighbor resampling. Next, a percentage cover (0%–100%) value was recorded for each of the eight land‐cover/land‐use types in all the aggregated 96‐m grid cells.

The third group of predictor variables included three forest structure and quality variables and eight CORINE land‐cover/land‐use variables, recorded at the landscape scale around the breeding sites. For the three hawk species, these variables were measured at a scale of 1 km, and for the three woodpecker species at a scale of 500 m, surrounding nest sites. These two scales were selected based on the known differences in the bird home ranges; that is, birds of prey have larger home ranges than woodpeckers.

The three forest structure variables recorded at the landscape scale were dominant tree age (years), forest‐stand volume (m3/ha) and forest site type (the four broad productivity classes, Table 1). We used the mean across the 96‐m grid cells within the 500 × 500 m, or 1 × 1 km square to represent landscape‐level tree age and volume, and the median for the dominant site type. For the eight CORINE variables, a moving window adjusted to either 500‐m2 or to 1‐km2 surrounding the nest site was employed to calculate their cover as a percentage value (0%–100%).

The two climatic variables were the mean temperature of January (T Jan, °C) and growing degree days (GDD5, °C), which yields the annual temperature sum of days with mean temperature above 5°C. We focused on these two climatic variables as they have been shown to be among the most important determinants of boreal bird species distributions in earlier studies (Virkkala et al., 2013). We used the fine‐grained climate data developed in Heikkinen et al. (2020), which provides the monthly temperature and precipitation data, averaged for the period from 1981 to 2010 at a resolution of 50 m, based on data sourced from 313 meteorological stations (European Climate Assessment and Dataset [ECA&D]) (Klok & Tank, 2009). These data reflect the local impacts of topoclimate across the whole study area. We resampled the 50‐m resolution data on T Jan and GDD5 into the 16 × 16 grid with ArcGIS, and then aggregated these to the 96‐m resolution by calculating the mean.

We used ArcGIS Desktop version 10.5.1 and QGIS 3.4.6 version Madeira in our spatial analyses.

Bird species modeling

We used MaxEnt modeling software (Elith et al., 2011; Merow et al., 2013; Phillips et al., 2006; Radosavljevic & Anderson, 2014) to develop predictions and maps of suitable nesting forest stands for the six bird species. Here, the nest sites of the bird species were treated as presence‐only data, comparable to species records derived from a natural history museum database, which rarely includes information on species’ absences in surveyed locations.

In the model building and in the variable‐selection process, we applied the framework outlined in Williams et al. (2012), where ecologically similar predictor variables are first grouped into two or more subgroups and tested for their importance, using backward elimination to exclude variables that do not explain significant levels of variation in the response variable. Significant variables from each subgroup are then combined and again tested for their importance via backward elimination and associated tests in MaxEnt. Here, we used this approach to reduce the number of predictor variables within three of the subgroups described above: (1) 8 forest‐stand structure and quality variables, (2) 8 land‐cover/land‐use variables recorded at the 96‐m focal (nesting) site, and (3) 11 landscape‐level measures of forest stands and land cover/land use, recorded with a 500‐m and 1‐km buffer around the nest site for woodpecker and hawk species, respectively. In addition, (4) the two local topoclimatic variables, January mean temperature and growing degree days (GDD5) (see Heikkinen et al., 2020) were included in the final predictive model and tested there.

Within the subgroups, we first built a model for each species, using all variables within the subgroup. The variables were tested using step‐by‐step backward elimination of the least important variables, primarily based on the known ecological importance of variables and training gain values calculated by MaxEnt, further guided by the percent contribution and permutation importance statistics. The variable importance testing was executed based on the 10‐fold cross‐validation at each step during the backwards variable elimination process. We also examined the individual variable response curves and compared them to a priori understanding of the species’ environmental relationships. The variable exclusion process was continued until all remaining variables showed a percent contribution and permutation importance of 1 or more (cf. Brambilla et al., 2015; Williams et al., 2012).

Next, the retained variables from the three subgroups, and the two climatic variables, were combined and jointly examined for their importance using the same process, starting from a full model with all variables. This provided the final MaxEnt model for each species, based on which nesting‐habitat suitability was predicted for the whole of Finland. The performance of the final predictive models was assessed using the area under the receiver operating characteristic curve (AUC) statistics, both for the training and test data, as calculated by MaxEnt (Phillips et al., 2006).

Model building also included decisions on correlations between predictor variables. Considering the robustness of MaxEnt against the impacts of high correlations based on the regularization method (Elith et al., 2011), we followed Fourcade et al. (2014) and excluded only the most highly (Pearson correlation test ≥0.9) correlated variables during the process of model building. However, exceptions to this rule were made when selecting the forest‐stand–level variables, some of which were strongly correlated with one another. Thus, following Clark et al. (2014), we allowed two highly correlated forest‐stand variables to be retained here in the model if they showed both considerable ecological significance and high performance in the MaxEnt jackknife tests.

By contrast, we included only one of the two climatic variables, January mean temperature or GDD5, in the final predictive model for each species. These variables were highly correlated with each other (Pearson correlation = 0.89), and also showed broadly redundant spatial trends. Where possible, the selection was done based on MaxEnt‐test values. However, if the test values were inconclusive, we selected the ecologically more relevant of the two based on earlier modeling results (Virkkala et al., 2013).

In all the models, default values were used, except that we used a bias grid in MaxEnt to correct for survey bias in the distribution of ringing sites. The ringing data of the six bird species were heavily concentrated in southern Finland, although all species (except the white‐backed woodpecker) also occur, albeit more sparsely, in northern Finland. As MaxEnt compares the environmental conditions of nesting sites to a randomly sampled set (called background points) across the study area (the whole of Finland), not accounting for survey bias can lead to areas in northern Finland being erroneously designated as unsuitable, when they are in fact poorly surveyed for these species (Phillips et al., 2009). To account for this, we built a kernel density layer (in R v. 3.6.1, with package MASS v. 7.3‐51.4) using the nest locations of the six species, supplemented with nest locations of three predominantly northern bird species (the rough‐legged buzzard [Buteo lagopus], the merlin [Falco columbarius], and the hawk owl [Surnia ulula]), and a 50‐km kernel distance (see Elith et al., 2011; Kramer‐Schadt et al., 2013; Kujala et al., 2015). The additional ringing data included 761, 514, and 640 nest locations of the rough‐legged buzzard, the merlin, and the northern hawk owl, respectively, during the years 1990–2017. In MaxEnt, the bias grid (the density layer) is used to weight the sampling of background points (Phillips et al., 2009).

As the final step of MaxEnt modeling, we predicted the nesting‐habitat suitability for each study species, ranging from 0 to 1, using MaxEnt's cloglog output format. To ensure maximal predictability, the final models were built using all available records for each species.

As MaxEnt models are based on presence‐only data, absolute distribution sizes cannot be derived from the predicted maps (Guillera‐Arroita et al., 2015). Therefore, species distribution patterns from MaxEnt predictions can only be examined in relative terms. We calculated the relative size of distribution of each species by summing the predicted values from the final model and dividing them by the number of grid cells in the study area. However, before estimating the relative distribution size and level of protection for each species, we set predicted values below a species‐specific threshold value, that is, maximum training sensitivity plus specificity (maxTSS; Liu et al., 2013) value provided by MaxEnt, to zero. The proportion of cells with a relative likelihood >0.8 (and as a comparison also >0.7 and >0.9) situated in PAs was calculated for each species.

RESULTS

Key habitat and land‐cover variables explaining species’ nesting sites

Forest‐stand–level characteristics at the nest site was the most important variable group in the MaxEnt models. Looking at the three top variables with the greatest percent contribution from each species model (altogether 18 variables), 11 were stand variables at the nest site level (1–3 variable[s]/species), 2 were land‐cover variables at the nest site, 3 were landscape‐level variables, and 2 were climatic variables (Table 2, Appendix S1: Figure S1). Accordingly, based on permutation performance of the 3 most important variables in the final models, there were 11 stand variables (1–2 variable[s]/species) at the nest‐site level, 1 land‐cover variable at the nest site, 4 landscape‐level variables, and 1 climatic variable (Table 2).

TABLE 2

Variables included in the final models for each species, both and permutation importance (%) are presented

Environmental variable	European honey buzzard		Northern goshawk		Common buzzard		White‐backed woodpecker		Lesser spotted woodpecker		Eurasian three‐toed woodpecker
	Contrib. (%)	Permut. (%)	Contrib. (%)	Permut. (%)	Contrib. (%)	Permut. (%)	Contrib. (%)	Permut. (%)	Contrib. (%)	Permut. (%)	Contrib. (%)	Permut. (%)
A
Stand volume	48.5 (+)	38.8	32.7 (+)	53.7	64.3 (+)	26.9	–	–	1.5 (+)	4.9	23.0 (+)	1.3
Stand dbh	–	–	2.7 (+)	0	–	–	4.0 (+)	10.9	–	–	–	–
Stand height	7.1 (+)	9.7	–	–	4.2 (+)	5.7	1.8 (+)	9.7	–	–	34.1 (+)	16.0
Dominant tree age	–	–	31.7 (+)	22.5	1.9 (+)	16.3	–	–	4.3 (+)	21.9	2.3 (+)	6.4
Deciduous tree volume	5.1 (+)	22.4	–	–	3.0 (+)	9.1	54.3 (+)	21.8	29.2 (+)	24.3	–	–
Main tree species	16.3 (S)	5.0	–	–	–	–	2.7 (D)	1.3	3.1 (D)	1.4	15.8 (S)	3.3
Forest site type	2.2 (1)	2.2	0.9 (2)	4.2	–	–	0.9 (1)	1.1	–	–	0.5 (2)	2.5
B
Shoreline forest	–	–	–	–	–	–	–	–	2.8 (+)	3.5	–	–
Agricultural areas	–	–	5.8 (−)	3.8	–	–	–	–	–	–	4.1 (−)	13.7
Urban areas	6.3 (−)	4.0	3.7 (−)	1.8	5.4 (−)	1.1	–	–	–	–	–	–
Water areas	–	–	2.1 (−)	3.9	5.1 (−)	4.8	–	–	–	–	–	–
C
Stand volume of forest compartments	–	–	1.0 (+)	9.0	1.2 (+)	3.8	0.6 (+)	3.6	–	–	4.3 (+)	34.4
Site type of forest compartments	5.1 (1)	11.3	–	–	1.6 (1)	6.0	12.0 (1)	8.5	8.8 (1)	4.6	1.3 (2)	1.5
Shoreline forest	–	–	–	–	–	–	2.1 (+)	2.7	4.4 (+)	1.1	–	–
Forest on peatland	–	–	–	–	–	–	2.0 (−)	7.8	5.3 (−)	11.4	–	–
Marshlands	–	–	–	–	–	–	5.3 (+)	3.4	16.2 (+)	5.7	–	–
Agricultural areas	5.9 (−)	5.7	–	–	1.4 (−)	1.7	–	–	–	–	–	–
Urban areas	–	–	6.6 (−)	1.1	5.0 (−)	8.2	–	–	13.0 (+)	17.9	9.0 (−)	6.8
Water areas	–	–	–	–	1.8 (−)	3.9	6.2 (+)	0.9	9.1 (+)	2.5	1.2 (−)	4.4
D
January mean temp.	–	–	12.7 (+)	0	–	–	8.2 (+)	28.2	2.4 (+)	0.7	4.3 (+)	9.7
GDD	3.4 (+)	0.8	–	–	5.0 (+)	12.7	–	–	–	–	–	–

Notes: A = stand‐level variables at the nest site, B = land‐cover/land‐use variable at the nest site, C = stand or land‐cover/land‐use variable at the landscape level, D = climatic variable. The direction of the effect of a variable [(+) or (−)] is presented. S = Norway spruce, D = deciduous tree. The most important forest‐site type: (1) herb‐rich and low‐herb heath forests, (2) mesic heath forests. The three most important variables are in bold, both in percent contribution and in permutation importance.

Abbreviations: Contrib., percentage of contribution; Permut., permutation importance; dbh, diameter at breast height; GDD, growing degree days.

Stand volume had the highest percentage contribution in the MaxEnt models in the European honey buzzard, the northern goshawk, and the common buzzard, and the second highest in the Eurasian three‐toed woodpecker (Table 2; Appendix S1: Figure S1). In addition, deciduous tree volume contributed the most in the models of the white‐backed woodpecker and the lesser spotted woodpecker (Table 2; Appendix S1: Figure S1). Stand height had the highest contribution in the model of the three‐toed woodpecker, and it was also included in the models of the honey buzzard, the common buzzard, and the white‐backed woodpecker. The main tree species was included among the three most important variables contributing in the models of the honey buzzard and the three‐toed woodpecker (Table 2, Norway spruce, Appendix S1: Figure S1).

The most important variables that contribute to the importance of permutation differed somewhat from those of percentage contribution. Similarly, as in the percentage contribution, stand volume was the most important variable based on permutation importance for the honey buzzard, the northern goshawk, and the common buzzard, and the deciduous tree volume for the lesser spotted woodpecker (Table 2). However, in the white backed woodpecker the January mean temperature was the most important variable, and deciduous tree volume was the second most important. Deciduous tree volume was also the second most important variable in the honey buzzard. In the three‐toed woodpecker, permutation importance differed most from the percentage contribution, with stand volume of forest compartments at the landscape level being the most important variable (Table 2). However, stand height, which was the most important variable in the percentage contribution, was regarded as the second most important variable in the three‐toed woodpecker based on permutation importance.

Predictions based on MaxEnt models

The model discrimination performance, as measured by MaxEnt's test AUC statistics, was clearly better than random (i.e., 0.5; see Phillips et al., 2009) for all six species, varying from 0.72 (the northern goshawk) to 0.97 (the white‐backed woodpecker) (Table 3). In general, woodpecker species had higher test AUC values (varying between 0.93–0.97) than hawk species (0.72–0.85). The maps of predicted nesting suitability showed clear spatial differences between species (Figure 3).

TABLE 3

Values of test AUC from the final model (SD in parentheses), the relative size of species distribution, and the proportion of high predicted values in PAs

	European honey buzzard	Northern goshawk	Common buzzard	White‐backed woodpecker	Lesser spotted woodpecker	Eurasian three‐toed woodpecker
Test AUC from the final model	0.845 (0.032)	0.717 (0.009)	0.828 (0.008)	0.971 (0.014)	0.926 (0.031)	0.932 (0.020)
Relative size of distribution	0.159	0.166	0.150	0.016	0.040	0.042
Proportion in PAs	0.107	0.319	0.109	0.124	0.078	0.263

Notes: Relative size of distribution was measured as the sum of predicted values from the final model, divided by the number of grid cells in the study area. The proportion of cells with the relative likelihood >0.8 situated in protected areas is presented for each species. The maximum predicted values for each species was 1.0.

Abbreviations: AUC, area under the curve; PA, protected area.

FIGURE 3

Predicted habitat suitability of (a) the European honey buzzard, (b) the northern goshawk, (c) the common buzzard, (d) the white‐backed woodpecker, (e) the lesser spotted woodpecker, and (f) the Eurasian three‐toed woodpecker in Finland

Hawk species had the largest areas projected to be suitable for nesting, the relative prevalence varying between 0.150 and 0.166 (Table 3). Suitable nesting sites for the European honey buzzard, the northern goshawk, and the common buzzard were predicted to occur much more evenly than for woodpecker species (Table 3, Figure 3). High MaxEnt predictions (>0.8) for the honey buzzard occurred particularly in the inlands in southern Finland and for the common buzzard along the southern coast (Figure 3). Suitable nesting sites for the goshawk were predicted over large areas in southern and north‐central Finland, and also sporadically in northernmost Finland (Figure 3).

The white‐backed woodpecker clearly had the smallest projected suitable area (0.016; Table 3). The most suitable nesting sites, based on model predictions >0.8, for the white‐backed woodpecker were scarce and concentrated in southeastern Finland along the lake shores where mature deciduous forests preferred by the species are situated (Figure 3). Projected nest sites for the lesser spotted woodpecker and the three‐toed woodpecker were distributed across the whole country, however, occurring rather sporadically (projected relative prevalence of 0.040 and 0.042, respectively; Table 3). Lesser spotted woodpecker habitats were concentrated in deciduous forests along lake shores, rivers, and coasts and in certain urban areas where marshlands are commonly located (Figure 3). Suitable habitats for the three‐toed woodpecker occurred sporadically, particularly in old‐growth coniferous forests largely located in PAs (Figures 1, 2, 3, 4). Highly suitable areas for nesting occurred in southwestern Finland, and in eastern and in northeastern (Kainuu–Kuusamo) Finland, primarily in PAs (Figures 3 and 4).

FIGURE 4

Predicted habitat suitability for the Eurasian three‐toed woodpecker in northeastern Finland (Kainuu–Kuusamo) in relation to protected areas

The proportion of cells with the highest nesting‐habitat suitability sites (>0.8) in PAs was highest for the northern goshawk (0.319) and the Eurasian three‐toed woodpecker (0.263), whereas for other species this varied between 0.078 and 0.124 (Table 3), indicating that a large proportion of suitable nesting forests occurred outside PAs. When varying relative likelihood values (>0.7 and >0.9), it can be observed that the higher the threshold value, the higher the proportion of suitable nesting forests was included in the PAs (Appendix S1: Table S1). This suggests that PAs indeed harbor highly suitable nesting forests for the six species relatively more often, particularly for the northern goshawk and the Eurasian three‐toed woodpecker.

DISCUSSION

We predicted suitable nesting habitat for the three hawk and three woodpecker species over large areas based on ringing and nest data. We were able to develop forest‐stand–level information on the distribution of important sites, and their underlying key environmental drivers, for our indicator species across extensive areas. This approach provided a way to determine both the nesting habitats of focal species and the valuable habitats of forest biodiversity in general over large boreal forest areas in northern Europe.

Ecological determinants of the study species’ nesting habitats

Our study shows that particularly stand volume and stand height related to high forest age are important variables for the high nesting‐habitat suitability of these indicator bird species. Decline of old‐growth forests and large trees are the most important causes of threat for the red‐listed forest species in Finland (Hyvärinen et al., 2019). Four of the indicator bird species were related to coniferous forests, which grow much older than deciduous forests (e.g., Virkkala & Toivonen, 1999) preferred by the two woodpecker species (white‐backed woodpecker and lesser‐spotted woodpecker). Urban areas, agricultural areas, or water areas had negative effect on the occurrence of species preferring coniferous forests, whereas for species preferring deciduous forests, water areas were positively affected, probably because of the location of mature deciduous forests along shorelines. Thus, our indicator species show versatile nesting‐habitat suitability patterns that are important in preserving boreal forests and their biodiversity.

Earlier studies have shown that the European honey buzzard, the northern goshawk, and the common buzzard all prefer old spruce forests in their breeding sites, and that the proportion of old spruce forests declined around their nest sites with increasing distance (Björklund et al., 2015; Byholm et al., 2020). The breeding success of the northern goshawk has also been shown to be highest in nest sites with a high proportion of old spruce forests in a 200‐m radius from the nest. In a study carried out in southern Finland, habitat suitability of goshawks increased with increasing volume of Norway spruce with small quantities of broad‐leaved trees (Björklund et al., 2020). These forests are subject to considerable logging pressure because about 10% of the predicted suitable habitats of goshawks was clear‐cut within a period of a few years within the studied region (Björklund et al., 2020; see also Byholm et al., 2020). It is very likely that clear‐cutting of suitable goshawk forests will also negatively affect other forest species with conservation importance. For example, goshawk nest sites have been shown to host more birds and polypores than reference sites, thus exhibiting a high level of biodiversity value in these sites (Burgas et al., 2014).

As noted earlier, the habitat preferences of the European honey buzzard, the northern goshawk, and the common buzzard overlap, but also show some differences. In our study, stand volume at the nest site was clearly the most important variable defining nesting‐habitat suitability for all the three hawk species, followed by tree age for the goshawk, and the main tree species (spruce) for the honey buzzard. The presence of herb‐rich and low‐herb heath forests at the landscape scale and the high volume of deciduous trees around nesting sites were important predictors for the European honey buzzard and common buzzard, respectively. These can partly be explained by the fact that the common buzzard often breeds near forest–field edges where deciduous trees are common and the European honey buzzard chooses areas close to water areas where more productive forests are typically located (Björklund et al., 2015).

Stand volume, stand height, and main tree species (spruce) were the main variables defining nesting‐habitat suitability in the Eurasian three‐toed woodpecker. The three‐toed woodpecker prefers old‐growth spruce forests with a significant abundance of dead wood (Lõhmus et al., 2010; Roberge, Angelstam, et al., 2008) because it eats wood‐boring insects (Fayt, 1999) and excavates its nest cavity typically in a dead or dying tree (Pakkala et al., 2018). The density of the three‐toed woodpecker was almost tenfold higher in a landscape of old‐growth forests (0.8 pairs/km2) than in managed forest landscape (0.1 pairs/km2) in north‐central Finland (Kainuu–Kuusamo; see Figure 4; Virkkala & Rajasärkkä, 2006).

The white‐backed woodpecker and the lesser spotted woodpecker breed in mature deciduous forests with plenty of dead wood as they both eat wood‐boring insects (Aulén, 1991; Lõhmus et al., 2010; Olsson et al., 1992; Virkkala et al., 1993; Wiktander et al., 1992). Thus, logically deciduous tree volume was the most important variable defining nesting‐habitat suitability of these species. The presence of marshlands in the surrounding landscape was an important variable for the lesser spotted woodpecker, probably because remaining mature deciduous stands are often restricted to marshland edges, which are typically left outside forestry practices. The predictive maps of suitable nesting sites for these species can provide very useful generic information for targeting new conservation actions, as both the white‐backed woodpecker and the lesser spotted woodpecker have been shown as indicators of high biodiversity, including many threatened beetles, in deciduous forests (Lõhmus et al., 2010; Martikainen et al., 1998; Roberge, Angelstam, et al., 2008; Roberge, Mikusiński, et al., 2008).

Advantages of methodology applied in the study

Among the numerous methods available for SDM, MaxEnt has been one of the most popular techniques globally to develop distribution maps for species (Bradie & Leung, 2017; Morales et al., 2017), with applications from many different species groups, including several studies on birds (e.g., Brambilla et al., 2015; Carroll, 2010; Moradi et al., 2019; Zhang et al., 2019). MaxEnt‐based assessment of variable importance has also been used to study population‐level responses (Searcy & Shaffer, 2016), such as colonizations and extinctions (Venne & Currie, 2021). However, our study differs from many of the earlier MaxEnt modeling works in two important ways. Firstly, most earlier MaxEnt studies involve combining species records with rather coarse‐scale land‐cover and climate variables. Here, we focused on local‐scale determinants of species distributions across large areas and combined species records with the fine‐grained information on forest characteristics, local land‐cover variables, and key climatic variables. By these means, we were able to develop detailed forest‐stand scale predictions of the suitable sites for the six indicator bird species systematically across the whole country. A key necessity in this is that the environmental data used as predictor variables is indeed systematically recorded with sound and comparable survey procedures across the whole fine‐grained lattice system (Kramer‐Schadt et al., 2013). When such environmental data are available, it has the potential to reveal suitable sites for the modeled species with much higher spatial accuracy than large‐scale SDMs, which are often based on more general environmental data such as broadly defined vegetation types (Bradie & Leung, 2017; Goetz et al., 2014; Seoane et al., 2004). Thus, we believe that our approach enables more robust local management and conservation planning.

Secondly, it is important to note that the habitat and environmental characteristics of the nest sites, where the studied bird species have successfully nested and produced offspring, can act as stronger biodiversity indicators than, for example, mere observations of adult birds. In this study, we used ringing data gathered from nest sites of boreal bird species that have different habitat requirements: the European honey buzzard, the northern goshawk, the common buzzard, and the Eurasian three‐toed woodpecker breed mainly in coniferous forests, whereas the white‐backed woodpecker and the lesser spotted woodpecker breed in deciduous forests. Thus, all the six species are indicative of different but likewise specialized forest conditions that are commonly in short supply in commercially managed forests, such as large stand volume and large coniferous and deciduous trees, high stand age, and accumulation of coarse woody debris (CWD), thereby providing useful indicators for conservation planning (Björklund et al., 2015; Burgas et al., 2016; Fayt, 1999; Pakkala et al., 2018). Importantly, because of the different habitat requirements of the six bird species, suitable nesting forests for these species can collectively cover a larger spectrum of forest biodiversity for conservation and management planning than each species could on its own.

The use of ALS and other RS data provided a highly important way of studying fine‐scale (resolution) habitat data at a large scale (extend). When ALS techniques are further developed for analyzing forest structure, for example, to separate dead wood important for many species, and ALS data are available for the whole country in Finland, this kind of fine‐resolution RS data related to nesting sites of bird species would make the data ecologically even more relevant (see, e.g., Klein et al., 2020).

Conservation issues

In Finland, old‐growth forests are focal habitats for the conservation of biodiversity, and specific protection programs for these habitats have been implemented (Auvinen et al., 2010). This is in line with the most recent Biodiversity Strategy in the European Union, where the aim is to protect all the remaining old‐growth forests by 2030 (European Commission, 2020). Currently, approximately 6% of the forest land (annual increment of at least 1 m3/ha) in Finland has been strictly protected and over half of all remaining old‐growth forests are included in PAs. Such PAs thus provide central sites for species conservation. However, great geographic differences exist in the cover of PAs, because more than 80% of the protected land is situated in northern Finland (National Resources Institute Finland, 2019; Virkkala et al., 2000), where the PA network is the most representative (Virkkala et al., 2000; Virkkala & Rajasärkkä, 2007) (Figure 1). At present, only about 3% of forest land is protected in the southern half of Finland and about 10% in the northern half of Finland (National Resources Institute Finland, 2019).

None of the studied species are concentrated in the northernmost regions with highest cover of PAs, as they have a more southerly distribution, and only the Eurasian three‐toed woodpecker is abundant in southern parts of northern Finland (Väisänen et al., 1998; Valkama et al., 2011). Nevertheless, PAs host considerable amounts of habitat for these species. In most parts of species’ ranges, less than 5% of (all) forests are protected, yet they host about 17%, on average, of the high suitability habitats for the species. Particularly species preferring old‐growth spruce forests, the northern goshawk (Björklund et al., 2020) and the Eurasian three‐toed woodpecker (Lõhmus et al., 2010; Roberge, Angelstam, et al., 2008) had a high percentage of their high‐suitability habitats (relative likelihood >0.8) inside currently existing PAs (32% and 26%, respectively).

Our analysis showed that there are also plenty of high‐quality habitats outside PAs, which should be taken account of in local and regional forest management planning by applying more biodiversity‐friendly forestry in these sites (Eyvindson et al., 2018). However, many of the high‐quality areas are fragmented, and therefore connectivity and contiguous cover of these areas should have a high priority in conservation planning. Unfortunately, because of recent accelerated pressures by the Finnish National Forest Strategy to increase intensive use of forest resources for bioeconomy and forest industry (Maa‐ ja metsätalousministeriö, 2015), many of these habitats are threatened by increased logging (Björklund et al., 2020; Byholm et al., 2020; Virkkala et al., 2020). Therefore, there is an urgent need to complement the existing PA network (METSO, 2020), as is also presented by the new Biodiversity Strategy in the European Union (European Commission, 2020) and for developing more biodiversity‐friendly logging procedures, such as, for example, continuous‐cover forestry, which benefits many species requiring mature and old‐growth forests (see Peura et al., 2018).

CONCLUSIONS

Identifying forest stands with high nesting‐habitat suitability for biodiversity indicator species, such as the six focal species in our study, can provide useful information for the future development of PA networks and targeting biodiversity‐friendly measures for forest management planning. The results of our analyses on habitat preferences of the study species were in line with the previous knowledge of species–habitat relations, which supported the reliability and validity of our modeling approach and the fine‐grained nest‐site suitability predictions it provided across the whole forested Finland. We conclude that predicting habitat suitability based on existing information on the nest sites of focal species provides a new tool of fine‐grained applicability for conservation planning, and this kind of approach would also be feasible and recommended elsewhere where similar data are available.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

Appendix S1

Click here for additional data file.