Remotely Sensed Single Tree Data Enable the Determination of Habitat Thresholds for the Three-Toed Woodpecker ( Picoides tridactylus )

: Forest biodiversity conservation requires precise, area-wide information on the abundance and distribution of key habitat structures at multiple spatial scales. We combined airborne laser scanning (ALS) data with color-infrared (CIR) aerial imagery for identifying individual tree characteristics and quantifying multi-scale habitat requirements using the example of the three-toed woodpecker (Picoides tridactylus) (TTW) in the Bavarian Forest National Park (Germany). This bird, a keystone species of boreal and mountainous forests, is highly reliant on bark beetles dwelling in dead or dying trees. While previous studies showed a positive relationship between the TTW presence and the amount of deadwood as a limiting resource, we hypothesized a unimodal response with a negative effect of very high deadwood amounts and tested for effects of substrate quality. Based on 104 woodpecker presence or absence locations, habitat selection was modelled at four spatial scales reﬂecting different woodpecker home range sizes. The abundance of standing dead trees was the most important predictor, with an increase in the probability of TTW occurrence up to a threshold of 44–50 dead trees per hectare, followed by a decrease in the probability of occurrence. A positive relationship with the deadwood crown size indicated the importance of fresh deadwood. Remote sensing data allowed both an area-wide prediction of species occurrence and the derivation of ecological threshold values for deadwood quality and quantity for more informed conservation management.


Introduction
Effective biodiversity conservation in managed forest landscapes requires knowledge about the distribution of key habitat features at relevant scales [1,2]. This knowledge is essential for assessing species' habitat selection, deriving threshold values for key features, and evaluating habitat quality across large spatial scales. Habitat suitability models (HSMs) [3] and their spatially explicit variant, Figure 1. For the purpose of this study, trees were classified as: living trees (LIVE) and standing deadwood objects (DEAD), which were further divided into dead trees (DEADTREE) and snags (SNAG) representing the stages of conifer tree decomposition after Thomas et al. [52]. Note the shrinking of the horizontal extension of the tree crown during the progress in decay.
Due mainly to harvesting and sanitation cutting, deadwood, to which TTW occurrence is closely linked, is an especially limited resource in managed forest ecosystems [21,54,64,65,73]. Deadwood thresholds of European forest-dwelling species range from 10 to 150 m 3 /ha with values of 20-50 m 3 /ha given for the majority of species as reviewed by Müller and Bütler [73]. This corresponds well with the 15-18 m 3 /ha to 30 m 3 /ha given for TTW occurrence [40,54,63,64]. Higher densities of deadwood are rare in Europe and occur only locally, mainly in protected areas [73,74], so that an upper deadwood limit could not be determined yet. However, the existence of a deadwood-optimum is likely, as a share of living trees would be necessary to allow for a continuous provision of dying and freshly dead trees.
In this study, we test the usability of remotely sensed single tree data for analyzing habitat selection and predicting area-wide occurrence of the TTW, identifying the most important variables explaining home range selection at multiple spatial scales, and deriving threshold values for conservation management. We test the hypothesis that extremely high amounts of deadwood lead to a decrease in the probability of TTW occurrence and assess the influence of deadwood quality on habitat selection. For the purpose of this study, trees were classified as: living trees (LIVE) and standing deadwood objects (DEAD), which were further divided into dead trees (DEADTREE) and snags (SNAG) representing the stages of conifer tree decomposition after Thomas et al. [52]. Note the shrinking of the horizontal extension of the tree crown during the progress in decay.
Due mainly to harvesting and sanitation cutting, deadwood, to which TTW occurrence is closely linked, is an especially limited resource in managed forest ecosystems [21,54,64,65,73]. Deadwood thresholds of European forest-dwelling species range from 10 to 150 m 3 /ha with values of 20-50 m 3 /ha given for the majority of species as reviewed by Müller and Bütler [73]. This corresponds well with the 15-18 m 3 /ha to 30 m 3 /ha given for TTW occurrence [40,54,63,64]. Higher densities of deadwood are rare in Europe and occur only locally, mainly in protected areas [73,74], so that an upper deadwood limit could not be determined yet. However, the existence of a deadwood-optimum is likely, as a share of living trees would be necessary to allow for a continuous provision of dying and freshly dead trees.
In this study, we test the usability of remotely sensed single tree data for analyzing habitat selection and predicting area-wide occurrence of the TTW, identifying the most important variables explaining home range selection at multiple spatial scales, and deriving threshold values for conservation management. We test the hypothesis that extremely high amounts of deadwood lead to a decrease in the probability of TTW occurrence and assess the influence of deadwood quality on habitat selection.

Study Area
The study was conducted in the Bavarian Forest National Park. Founded in 1970 as the first German National Park, it initially covered an area of 13,300 ha which was extended to 24,218 ha in 1998. The park is located in southeastern Bavaria (Germany) and borders the Šumava National Park (69,030 ha), Czech Republic to the East. The park covers a large part of the Bavarian Forest mountain chain with an elevational gradient ranging from 600 to 1453 m a.s.l. Depending on elevation, the mean annual temperature (1972-2001) varies from 3.5 to 7.0 • C, and the total annual precipitation varies from 1300 to 1900 mm [75]. The predominant vegetation is mountainous spruce and mixed forest with a share of Norway spruce (Picea abies) of 67.0%, European beech (Fagus sylvatica): 24.5%, Silver fir (Abies alba): 2.6%, and other tree species: 5.9% [76].
Following its non-intervention policy, the National Park authority allowed for natural forest dynamics in the core zone (currently encompassing 68% of the park area), with massive bark beetle outbreaks after severe storm and windthrow events in 1983 and 1984. This resulted in a dieback of spruce forests at an unprecedented rate in Central Europe in recent decades [77].

Remote Sensing Data
Habitat variables were extracted from a full, remote sensing-based tree inventory [78]. Full waveform ALS data was acquired on the 24th/26th/27th of July 2012 through the Milan Flug GmbH, using a Riegl LMS-Q 600i laser scanner of 350 KHz. A nominal point density of 30-40 points/m 2 was obtained from data recorded at a 0.32 m footprint. CIR aerial images were acquired in August 2012 using a DMC camera and a ground sampling distance of 20 cm. The images are composed of 3 spectral bands: near infrared, red, and green.
The preprocessing of the raw ALS data to the georeferenced three-dimensional (3D) point cloud, including the derivation of the intensity and the pulse width values using a sum of Gaussian functions, is described in Reitberger et al. [79] and Yao et al. [80]. Single tree detection and delineation was carried out with a 3D segmenting method solely based on ALS data and the geographical position and top height (H) were calculated for each segmented tree [81]. This resulted in a dataset containing 12,106,320 trees, consisting of two types of geographical data, point data for tree tops and polygons for crown delineation. In the next steps, for each tree, the tree type (conifer, broadleaf, or deadwood), projected crown area (C), crown base height and crown volume, were derived using both types of remote sensing data. Spectral information from CIR aerial imagery fused with segmented ALS point cloud data was used for tree species classification based on Reitberger et al. [79] and for the detection of snags and standing deadwood in line with Polewski et al. [30] and Polewski et al. [82]. The crown base height and the crown volume originated from a 3D Alpha-Shape-triangulation of the segmented ALS point cloud. Diameter at breast height (DBH), basal area (BA) and volume (VOL) were also calculated for live trees in an extra modelling step based on a calibration with an extensive ground reference database [83].
In addition to the automatically derived, tree-based information, we tested independent data of yearly visual assessment of deadwood areas based on the CIR aerial imagery [84] dated 2010-2015 (Table A1).

Species Data and Sampling Design
Presence locations of the TTW originated from the database of the Bavarian Forest National Park, including data from the biological monitoring, various research projects, and chance observations by trained park staff ( Figure 2). Locations of TTW, either observed directly or through sound identification, were recorded with a GPS. Data were collected in two time periods (2007-2008 and 2012-2014), however, to achieve a temporal synchronization with the remote sensing data from 2012, we only used the observations from the latter sampling period (N = 115).  [49,66,85,87].
The presence locations showed spatial clumping in some regions, indicating multiple observations originating from the same individual. To avoid pseudoreplication, we thinned the initial set of presence locations allowing a maximum 18% overlap of the sampling plots at the largest scale (r = 600 m), after Pechacek [85] who reported average territory overlaps of 17.6% (±3.9). Using R-package "Spatstat" package [88] to discard all presence locations that fell below the resulting minimum distance of 840 m, resulted in a final set of 52 presence locations.
In addition, we randomly created a similar number of pseudo-absence locations (in the following referred to as "absence") with the same minimum distance to any presence location and to each other.  [49,66,85,87].
The presence locations showed spatial clumping in some regions, indicating multiple observations originating from the same individual. To avoid pseudoreplication, we thinned the initial set of presence locations allowing a maximum 18% overlap of the sampling plots at the largest scale (r = 600 m), after Pechacek [85] who reported average territory overlaps of 17.6% (±3.9). Using R-package "Spatstat" package [88] to discard all presence locations that fell below the resulting minimum distance of 840 m, resulted in a final set of 52 presence locations.
In addition, we randomly created a similar number of pseudo-absence locations (in the following referred to as "absence") with the same minimum distance to any presence location and to each other.

Predictor Variables
Predictor variables were generated for each of the four sampling plot scales and encompassed three classes: general forest stand characteristics, specific tree features, as well as topographic information (Table 1).

Predictor Variables
Predictor variables were generated for each of the four sampling plot scales and encompassed three classes: general forest stand characteristics, specific tree features, as well as topographic information (Table 1). Based on a literature review, predictor variables describing general forest stand characteristics were selected according to their hypothesized ecological relevance for the species. Forest cover per plot was defined as the share of the horizontal plot area covered by all trees' crowns. The proportion of forest cover attributed to trees higher than 15 m was assessed as a proxy for mature forest, and the proportion the sampling plot covered by crowns of standing deadwood, as an indicator for areal deadwood availability. We also included a similar variable derived from the standard visual mapping (CIR): the proportion of area with standing deadwood originating from the period 2010-2015. This data was compared with the available ALS and CIR based data. In addition, the number and proportion of living trees (i.e., conifers, deciduous trees, and all trees) were calculated for each plot size, as well as the average height (H), diameter (DBH), and volume (VOL) and the variance thereof.
Specific tree features such as dead, "resource" or "veteran" trees have been observed to be important for the TTW. We defined resource trees (RESOURCE) as all living trees with a DBH > 30 cm, based on Pechacek and d'Oleire-Oltmanns [48], Kajtoch and Figarski [64] and Kajtoch et al. [54]. To approximate the different decay stages of deadwood in adherence to Thomas et al. [52] (Figure 1), and to distinguish between potential foraging trees (stages 2-4) and other dead trees (stages 5-7) as proposed by Bütler et al. [63], we subdivided the deadwood (category DEAD including all standing deadwood objects) into snags (SNAG) and dead trees (DEADTREE). The category SNAG encompassed all deadwood with either a crown area of C < 4 m 2 (1st Quartile of the DEAD crown area values) or deadwood with C ≥ 4 m 2 but a height of H < 15 m. DEADTREE objects were characterized as C ≥ 4 m 2 and H ≥ 15 m. The threshold of 15 m was chosen because in our study area, living spruce trees of that height had an approximate DBH of 20 cm ( Figure A1), allowing a direct comparison with DBH-based classifications of deadwood used by other authors (Bütler et al. [63]). Finally, we calculated the mean crown area of the trees in the respective classes (DEAD, DEADTREE, and SNAG) per plot to obtain a continuous metric of the crown status as a proxy of the stage of decay.
All structural predictors, except the variable derived from visual interpretation (DEADCIR_part), were thus calculated based on the single tree data originating from the ALS and CIR remote sensing dataset, either used directly to describe specific tree features or aggregated to describe forest stand related characteristics.
In addition, topographic variables were generated using a digital terrain model (DTM) with a 25 × 25 m resolution and included altitude, slope, eastern (sine of aspect), and northern (cosine of aspect) exposition as well as solar radiation, calculated using the Solar Analyst module in ArcGIS. We also included latitude and longitude to test for random spatial clustering of the woodpecker observations.
The preparation and calculation of variables with a horizontal dimension (i.e., referring to the proportion of the plot area) was carried out in ArcMap 10.4.1. [89]. The processing and calculation of the remaining variables was carried out in RStudio [90] using R [91] with the packages: "Raster" [92] and "Rgdal" [93].

Statistical Analysis
We modelled species occurrence as a function of the environmental variables using General Additive Models (GAM) facilitated in the R-package "mgcv" [94][95][96]. GAMs combine General Linear Models with smoothing splines [97], thereby allowing to fit the response curves "as closely as possible" to the data, within a permitted level of smoothing.
For each plot size we selected the best explaining variables following a hierarchical procedure. First, we ran univariate models for each variable, also testing their quadratic term. Predictor variables which significantly explained woodpecker occurrence (p ≤ 0.05) or showed a trend (p ≤ 0.1) and were significant in other studies (such as information on amount and volume of conifers and amount of live trees per plot) were retained. To avoid collinearity among variables in the multivariate models, we removed from any pair of correlated variables (Spearman's r ≥ |0.7|) [98] the "weaker" predictor based on Akaike's Information Criterion (AIC) [99].
In a first step, a common initial set of input variables (i.e., all variables retained at any of the four scales) was used for model calibration on all plot sizes. In addition to the environmental variables smoothed with a smooth term s(), a tensor smooth for the spatial location te(x, y) was added to account for effects of random spatial clustering of the TTW data ( Figure A2).
To avoid overfitting, as it was observed when running the model with automatic settings, we set the upper limit of the degrees of freedom associated with a smooth term to k = 3, as Guisan et al. [100] recommends after Hastie et al. [51] to use lower degrees of freedom (df < 4) for deductive species distribution and habitat modelling, while avoiding degrees of smoothing higher than 4 or 5 for predictive purposes. We used automatic variable selection (function "select=TRUE" in mgcv) which indicates variables that do not contribute to the model and can therefore be dismissed with p = 1. After the removal of these variables, the models were recalibrated and variables were again removed until no p = 1 occurred. In a second step, we used chi-square test statistics for assessing the significance of the smooth terms and removed variables with p < 1 but with chi-square equal 0 as not contributing to the model. At each step we compared the AIC of the resulting model with the previous step until no further reduction was achieved This way, we obtained one "best model" for each of the 4 plot sizes.
The models' fit was evaluated using 5-fold cross validation with 20% of the observations held back randomly with a condition of an equal proportion of presence and absence observations in folds. Multiple evaluation metrics, i.e., sensitivity, specificity, correctly classified rate, and Cohen's Kappa (all using the threshold 0.5), as well as the area under the receiver operating characteristics (ROC) curve (AUC) were calculated using the "caret" package in "R" [101] and evaluated according to Hosmer and Lemeshow [102]. The best model was then used to predict TTW occurrence probability across the entire National Park.
To analyze the model results, we plotted the single predictor variables against their smooths (function "gam.check" in the "mgcv" package) and against the target variable using the packages "mgcv" and "ggplot2" [103].
Finally, we calculated conditional inference trees (CTREEs), as implemented in the R-package "party" [104], to obtain thresholds for the most important variables for practical management recommendations. Trees based on maximally selected rank statistics were fitted using the Bonferroni correction for multiple testing and a minimum sum of weights in a node to be considered for splitting of 20 (minsplit = 20). All variables selected for the respective "best" GAM at each scale were included in the multivariate trees. In addition, univariate trees were fitted for variables with a significant split in the multivariate tree.

Habitat Selection
Univariate models revealed eight environmental variables that significantly contributed to explaining woodpecker presence at least one of the four scales and were retained for calibrating multivariate models (Table 2). With the exception of altitude all of the significant variables described stand and tree-related habitat characteristics. Three additional variables (DEAD_Nha, DEADTREE_Cmean, and DEADRSI_part) were significant but discarded as correlated with retained variables.
The final models consisted of 4 to 6 variables, depending on the plot scale ( Table 3). The amount of dead trees was the only variable with a significant contribution at all scales. Altitude and the mean crown diameter of all standing deadwood were also included in all models, but the former was only significant at the two smaller scales, while the latter was only significant at the two intermediate scales.
The proportional volume of conifer trees was included at 3 scales and the mean crown diameter of snags only on sampling plots with a 250 m radius, but neither of them had a significant contribution. The spatial location was also included in all models, suggesting a spatially clustered distribution of the woodpecker observations. The two preselected variables representing the living stand (the amount of living trees per ha and the amount of conifers per ha), although univariately significant at 3 scales, did not contribute to any of the multivariate GAMs.  The most meaningful variable at all scales was the amount of dead trees (DEADRTREE_Nha). The response plots ( Figure 3 and Table A2) indicate a unimodal response with adverse effects on woodpecker presence when the amount of dead trees increased beyond a threshold of 40-55 trees per hectare. However, these results need to be interpreted with caution due to only a few plots with extremely high numbers of dead tree driving this trend (i.e., two sample plots with DEADTREE_Nha > 90 on R = 100 m and R = 250 m and three sample plots with DEADTREE_Nha > 70 on R: 250, 450, and 600 m).
Another deadwood variable, the mean crown area of all standing deadwood (Dead_Cmean), showed a significant positive effect on TTW occurrence at the two intermediate scales (R = 250 and R = 450 m), suggesting a preference for deadwood with large crowns, i.e., in the early stages of decay ( Figure 1). This is in line with the species' opposite response trend to the mean crown area of snags (SNAG_Cmean) at the intermediate scale (R = 250 m). Finally, the share of conifers in the total volume of living trees (CONIF_VOLpart) showed a slightly unimodal, but non-significant relationship with TTW presence, with the highest occurrence probabilities in stands with about 50-80% conifers, depending on the sampling scale. Altitude was included in all models and significant at the two smallest scales, with higher probabilities of TTW presence at higher altitudes.

Model Performance
Model performance decreased with increasing sampling scale. This trend applied to both model fit and predictive performance over the 5-fold cross validation replicates and was consistent across all evaluation metrics (Table 4). Based on the AUC, our final models showed a good to excellent fit at the two small scales (R = 100 m and R = 250 m, AUC > 0.8) and an acceptable fit at the two larger scales (R = 450 and 600, 0.7 < AUC < 0.8) [102]. Five-fold cross validation confirmed an acceptable predictive performance of the models at the two smallest scales (AUC > 0.7), but less so at the two  Table 3). The blue line indicates the estimated smoothing parameter of a given variable, while keeping all other variables set on the median. Shadowed areas indicate the 95% confidence intervals conditional on the estimated smoothing parameter. Variable codes and descriptions are listed in Table 1.

Model Performance
Model performance decreased with increasing sampling scale. This trend applied to both model fit and predictive performance over the 5-fold cross validation replicates and was consistent across all evaluation metrics (Table 4). Based on the AUC, our final models showed a good to excellent fit at the two small scales (R = 100 m and R = 250 m, AUC > 0.8) and an acceptable fit at the two larger scales (R = 450 and 600, 0.7 < AUC < 0.8) [102]. Five-fold cross validation confirmed an acceptable predictive performance of the models at the two smallest scales (AUC > 0.7), but less so at the two larger scales (0.6 < AUC < 0.7). Similar trends were found for the other evaluation metrics R-Square (Adjusted), Correct Classification Rate (CCR), and Cohen's Kappa (Table 4. Sensitivity was generally higher than specificity at all scales indicating better detection and prediction of TTW presence than of absence locations. Complete results of the cross validation are presented in Appendix A, Table A2.

Model Prediction
The model calibrated at the smallest scale (R = 100 m) was employed to predict TTW occurrence throughout the National Park. The results are shown for a raster of 100 × 100 m (Figure 4), with the occurrence probability of each raster cell calculated based on the conditions within R = 100 m around the grid cell center. When using a presence probability of 0.5 as a threshold for occurrence, 36% of the park area was classified as potentially suitable TTW habitat. larger scales (0.6 < AUC < 0.7). Similar trends were found for the other evaluation metrics R-Square (Adjusted), Correct Classification Rate (CCR), and Cohen's Kappa (Table 4. Sensitivity was generally higher than specificity at all scales indicating better detection and prediction of TTW presence than of absence locations. Complete results of the cross validation are presented in Appendix, Table A2.

Model Prediction
The model calibrated at the smallest scale (R = 100 m) was employed to predict TTW occurrence throughout the National Park. The results are shown for a raster of 100 × 100 m (Figure 4), with the occurrence probability of each raster cell calculated based on the conditions within R = 100 m around the grid cell center. When using a presence probability of 0.5 as a threshold for occurrence, 36% of the park area was classified as potentially suitable TTW habitat.

Variable Thresholds
The results of the conditional inference trees ( Figure 5) revealed significant thresholds for two variables, the mean crown area and the amount of standing deadwood per plot. Multivariate trees were only found at the two smallest scales, with a first split indicating the highest TTW presence probability (>0.7) when deadwood with large crowns (>11-13 m 2 ) was available, intermediate probabilities (>0.5) when the abundance of dead trees per hectare was at least 4-5 (R = 100) or 3 (R = 250) respectively, and a low probability when none of the two variables exceeded these thresholds. These results indicate substrate selection with a first priority for fresh and then for other deadwood. No split was observed for variables measured on plots of 450 m radius or larger. Univariate models of the two variables with significant splits showed a TTW-presence probability of 0.7-0.8 when more than 8 dead trees per hectare were present in the surrounding of 100 and 250 m, respectively, or when the mean crown size DEAD_Cmean was larger than 11 m 2 (R = 100 m) to 13 m 2 (R = 250).

Variable Thresholds
The results of the conditional inference trees ( Figure 5) revealed significant thresholds for two variables, the mean crown area and the amount of standing deadwood per plot. Multivariate trees were only found at the two smallest scales, with a first split indicating the highest TTW presence probability (>0.7) when deadwood with large crowns (>11-13 m 2 ) was available, intermediate probabilities (>0.5) when the abundance of dead trees per hectare was at least 4-5 (R = 100) or 3 (R = 250) respectively, and a low probability when none of the two variables exceeded these thresholds. These results indicate substrate selection with a first priority for fresh and then for other deadwood. No split was observed for variables measured on plots of 450 m radius or larger. Univariate models of the two variables with significant splits showed a TTW-presence probability of 0.7-0.8 when more than 8 dead trees per hectare were present in the surrounding of 100 and 250 m, respectively, or when the mean crown size DEAD_Cmean was larger than 11 m 2 (R = 100 m) to 13 m 2 (R = 250).

Discussion
Our analysis shows the usability of area-wide, remote-sensing-based, single tree data for modelling the habitat selection of an endangered and highly specialized forest species such as the three-toed woodpecker. Combining remote sensing information from different sources with a comprehensive set of species observation data enabled finding species-relevant predictor variables and thresholds for practical forest management and species conservation.

Discussion
Our analysis shows the usability of area-wide, remote-sensing-based, single tree data for modelling the habitat selection of an endangered and highly specialized forest species such as the three-toed woodpecker. Combining remote sensing information from different sources with a comprehensive set of species observation data enabled finding species-relevant predictor variables and thresholds for practical forest management and species conservation.

Remote Sensing Data
The fusion of multiple sources of remote sensing data, especially ALS and aerial imagery technologies, has high potential for complementing traditional, field-based forest inventories [10]. Although original ALS point clouds deliver more detailed information on tree height and forest structure [105] and are widely used as input data for habitat modelling [106], the combination of ALS data with aerial imagery for deriving single tree-related information proved crucial for our purpose: While ALS data enabled accurate mapping of single trees and their projected crown areas with subsequent modelling of tree volume, multispectral data allowed deadwood detection and-in combination with the structural information-the detection of specific deadwood characteristics, such as fresh deadwood and snags. This approach offered two additional advantages: First, by summarizing structural information at the tree-level, our variables refer to a species-relevant ecological scale of habitat selection. Second, other than abstract point-cloud metrics, our data describe environmental features that can be directly translated into target values for conservation management.
Our deadwood variables at the tree scale outperformed the deadwood information (area-percentage of deadwood per plot, DEADCIR_part, Table 3) obtained from the yearly visual assessment of aerial imagery. Although univariately significant at the smallest plot size, it did not enter the final model. Only new deadwood areas were mapped each year [84,107], thereby neglecting forest dynamics such as ingrowth and regeneration in the dieback areas of previous years. Complete mapping of standing deadwood for a given year may therefore have improved the performance of this variable.
The corresponding variable based on remote sensing tree inventory data, the area percentage of deadwood per plot (DEADRSI_part), showed better explanatory power (lower AIC than DEADCIR) on 100 and 250 m plot sizes, but was correlated with the number of dead trees (DEADTREE_Nha), our most important predictor. It was therefore discarded. Nevertheless, the relationship of the two variables with TTW occurrence shows some potential for the planar mapping of standing deadwood areas when single tree crown delineation is not possible.
We show the usability of remotely sensed single tree data and derived variables using the example of the TTW, a keystone species of boreal and mountainous spruce dominated forests. However, these data could also be of high relevance for modelling the habitat of other species or species assemblages of that forest types [21,38]. Information on deadwood features and their quality (dead trees, snags and stumps) could be vital for species depending on deadwood in different decay stages either for food, shelter, or roosting such as saproxylic beetles [108], birds [46], or bats [109,110] and crown delineation, allowing the determination of canopy cover and forest gaps, could be used in studies deriving habitat thresholds for species responding to these structures e.g., capercaille [111] or hazel grouse [112].

Species Data
TTW presence data originated from three survey projects, including non-systematically collected chance observations of park staff. Despite thinning the original data according to the expected home range for one pair of birds, the models still indicated a clustering of the observations and a spatial correlation of the model performance with the locations of the observations. This may reflect a bias in sampling intensity, e.g., related to the road and walking paths network in the National Park, or be caused by a species-relevant environmental variable not included in the model. Our final models showed notably higher sensitivity than specificity, indicating a better classification of presence than absence data. This may be because of the random generation of pseudo-absence data outside the TTW presence areas, where false absences could not be ruled out.
We used species observations from three consecutive years starting in year one of the remote sensing data acquisition. At this time, the area of the National Park offered a large range of conditions, including optimal TTW habitat of mountainous, spruce dominated forests with a large amount of standing deadwood in different stages of decay. As both species and environmental data originated within a limited period of time and a unique environment, our models reflect only a snapshot of the species-habitat relationship [3]. The time lag of two years between the acquisition of remote sensing and species data we consider negligible, as also demonstrated by Vierling et al. [113], since no significant changes due to disturbance events were recorded in the respective period and the National Park is not subjected to regular harvesting. Moreover, our models showed a high predictive performance, with results largely conforming to those of other studies. This makes us confident that they captured TTW habitat requirements with a high level of generality.

Modelling Approach
GAMs are increasingly used in ecological modelling, especially when species-habitat relationships are complex and not easily fitted with the standard parametric functions of the predictors [100]. Using GAMs for the exploratory analysis of predictor variables is advantageous as GAMs fit the data in the most exact way possible [97]. However, being data-driven, they are prone to overfitting. We applied stronger smoothing to address this issue. The most important feature of GAMs for our study was the possibility of including a multidimensional smoother [100] for the spatial location (x,y) of TTW observations to account for spatial clumping of the data.
Although useful for identifying key variables and describing the species response, GAMs do not provide threshold values which are frequently required in ecology and forestry to define conservation targets [114]. We used conditional inference trees for this purpose as Müller and Bütler [73] found them particularly useful among a variety of methods [115]. The simplicity of the underlying model and the visualization of the results facilitate the development of applicable guidelines.

TTW Habitat Selection
From the initial broad set of environmental predictors (Table 1), only four structural variables indicating food resources, cavities, and altitude affected the occurrence of the TTW in our study area. Dead tree abundance was the most important variable. The species had a preference for deadwood in the early stages of decay when the abundance of insect food is highest [116,117]. Dying and dead spruce trees provide the major food sources of the TTW due to bark beetles (esp. Ips typographus) and wood-boring longhorn beetles inhabiting them. Müller and Bütler [73] showed the probability of TTW presence increasing from 0.1 to 0.9 when more than 0.81 (0.56-1.22, Switzerland) and 0.44 (0.25-0.62, Sweden) m 3 /ha basal area of standing deadwood corresponding to approx. seven and four dead trees with DBH ≥ 21 were present. Our results of 8 or more dead trees per hectare resulting in an 80% probability of TTW presence are in accordance with these findings.
In contrast with previous findings focusing on the minimum deadwood threshold, we show that very high amounts of deadwood, especially of late decay stages with little foraging value, negatively affect TTW occurrence probability. The detection of this tendency was made possible by very few observations at the extreme end of the gradient (i.e., sites with up to 120 trees per hectare), stemming from the large-scale area-wide bark beetle infestations. The lack of suitable research areas in Europe exhibiting the full possible gradient of deadwood abundance may be the reason that this effect has remained undetected, although Scherzinger [118] observed a recession in TTW occurrence, a few years after a significant increase following the bark beetle outbreak. This implies that a patchy distribution of bark-beetle infested trees and tree groups in the forest landscape is favorable compared to large-scale area-wide dieback, which is more likely in homogeneous, even-aged stands. Such heterogeneous deadwood distributions may be furthered by natural topographic complexity and increasing forest structural variability through active management or strict protection [119], as structural heterogeneity is expected to increase in unmanaged forests [120].
We also found a positive effect of the mean crown area of the dead trees per plot, indicating the availability of fresh deadwood with still complete tree crowns. This variable was selected into all models, although only significant at the two intermediate scales. Conditional inference trees indicated high probabilities (0.7-0.8) of woodpecker occurrence when the mean crown area per plot was larger than 11 m 2 (R = 100 m) or 13-13.5 m 2 (R = 250-450 m), respectively, corresponding to an average branch length of about 2 m. These findings are in line with Balasso [53], who found that the presence of TTW was related to abundance of fresh snags, and Scherzinger [84], who reported an initial increase in TTW occurrence shortly after bark beetle infestations with a subsequent decrease after some years. Nevertheless, the relationship between remotely sensed crown parameters and bark conditions needs further research.
In most field-based studies, deadwood is classified into standing dead trees, snags, and logs, representing different decay stages to account for TTW's prey diversity. In our study, the input data was limited to information that can be derived from the air. The first limitation was the omission of logs, the recognition of which, although theoretically possible by using ALS data from scanning in leaf-off conditions [121], was impossible with our data. Studies relying on field data often included this variable in HSMs [40,54], however it was rarely significant [51]. In addition, our remote sensing data could not provide information about the DBH, basal area (BA) and volume (Vol) of deadwood objects as often used in other studies [51,53,54,73,122]. This was due to the difficulty of modelling these values without reliable height measurements of the tree tops that are often broken in standing deadwood.
Similar to the findings of Braunisch et al. [21], our study suggested a positive, but non-significant correlation of TTW occurrence with the presence of conifers. In the Bavarian Forest National Park, conifer trees are predominantly Norway spruce, the primary host tree of Ips typographus which is the staple food of the TTW [123]. Scherzinger [118] concluded that not the pure amount of deadwood, but a permanent occurrence of dying and freshly dead trees originating from a continuous share of live spruce stands are crucial for the presence of TTW in the area. Mapping still alive, but degenerating spruce trees (the so called green attack stage) that were not detectable from our data and that remain a challenge for the remote sensing research [124][125][126] could potentially be of high explanatory value for TTW habitat selection. Further research using hyperspectral data could bring important progress here [127][128][129]. Resource trees that were an important variable in other studies [48,54,64] did not correlate with TTW occurrence in our study, due to a similar, very high resource supply in both presence and absence plots over the entire study area.
Decreasing model performance from the smallest to the largest sampling scale indicates habitat conditions, especially the amount and quality of deadwood, in the surrounding approximately 20 ha are most decisive for the TTW's habitat choice [17]. As species' area requirements depend on habitat quality, TTW home range sizes have been shown to vary considerably among regions and foraging conditions [47,57,62,85]. Bütler et al. [63] reports TTW ranges vary between 44 and 176 ha, depending on food availability and snag abundance. Kajtoch et al. [66] suggests at least 100 ha with optimal conditions and 200 ha in suboptimal stands are necessary, conforming to the results of other studies [43,58,[85][86][87]. Our study area, with its consistently high abundance of patchily distributed deadwood in different stages of decay therefore seems to represent an optimal habitat for the TTW.

Management Recommendations
Effective forest and biodiversity management requires habitat thresholds at a scale and resolution that are ecologically relevant to the species and can be practically implemented [105]. Bütler et al. [63] recommended a precautious 1.6 m 2 (basal area), corresponding to 5% of all standing trees or 14 standing dead trees with a DBH ≥ 21 cm per ha. We show the best response of TTW to habitat features within 100 to 250 m, i.e., related to a surrounding of up to 20 hectares. Within this area, at least eight dead trees per hectare should be retained, focusing on fresh deadwood in the early stages of decay, indicated by an average branch length of at least 2 m. Pechacek and Krištín [44] give similar management recommendations claiming that "dead trees should not be removed within a 250 m circle from nests".
To favor the coexistence of alternative prey for TTW and ensure a constant input of fresh deadwood, retaining and restoring dead coniferous trees in different stages of decay and a substantial portion of live spruce trees is required. At the landscape scale, Bütler et al. [40] showed an effect of the spatial arrangement and density of deadwood rich patches, and recommended a network of forest stands with high deadwood densities embedded in a forest landscape with lower deadwood densities.
As bark beetle spread was revealed to be strongly distance dependent with the most new infestations occurring within a 250 m radius of the previous year's infestation and 95% thereof in 500 m [130], safeguarding a wide enough inter-patch distance is crucial for preventing of bark beetle outbreak. Patches of declining and dead trees large enough to host bark beetle populations but disconnected from each other, would therefore aid forest managers in effectively controlling bark beetle dispersion [131], while at the same time promoting woodpecker habitat.
Our map showing current TTW habitat suitability allows distinguishing deadwood rich versus deadwood poor areas, so as to accurately target conservations measures.

Conclusions
Our study highlights the value of remote sensing, especially the fusion of ALS data with digital aerial imagery, for generating a full inventory of live and dead standing trees for large-scale, area-wide habitat analyses. Combining structural and spectral data enabled not only the identification of deadwood, but also of deadwood characteristics, which is indispensable for reliably modelling the habitat requirements of species highly specialized on particular types of standing deadwood. While our habitat analysis confirms the amount of standing dead trees as a key predictor of TTW occurrence, and the species' preference for fresh deadwood characterized by large and intact crowns, our study is the first showing a negative impact of very high deadwood amounts, with a tipping point at about 40-55 standing dead trees per ha. Moreover, we highlight the importance of resource diversity including also snags and live conifers. Based on tree-related remote sensing information, we were able to draw management recommendations. For example, keeping at least eight dead trees in the early stages of decay per hectare within 20 ha (corresponding to a small woodpecker's home range) leads to an increase in habitat suitability for the TTW.
Our models show a high predictive power, nevertheless, they may be improved by a more precise separation of fresh and old deadwood or even a further differentiation of decay stages or deadwood quality obtainable from field studies. Comparing decay stages from field assessments with time series of remote sensing data, and using hyperspectral imagery to detect tree decline in an early stage (e.g., the first stage of a bark beetle infestation), may further advance the set of predictors and aid foresters to better identify and carry out effective management measures to support biodiversity.

Conflicts of Interest:
The authors declare no conflict of interest.      Figure A2. Variable smooth effect plots for the predictor variables at all spatial scales produced using "gam.check". In brackets on y-axis: variables' edf (the estimated degrees of freedom of the smooth') from the GAM model. Variable codes and descriptions are listed in Table 1.