Multilevel landscape utilization of the Siberian flying squirrel: Scale effects on species habitat use

Abstract Animals use and select habitat at multiple hierarchical levels and at different spatial scales within each level. Still, there is little knowledge on the scale effects at different spatial levels of species occupancy patterns. The objective of this study was to examine nonlinear effects and optimal‐scale landscape characteristics that affect occupancy of the Siberian flying squirrel, Pteromys volans, in South‐ and Mid‐Finland. We used presence–absence data (n = 10,032 plots of 9 ha) and novel approach to separate the effects on site‐, landscape‐, and regional‐level occupancy patterns. Our main results were: landscape variables predicted the placement of population patches at least twice as well as they predicted the occupancy of particular sites; the clear optimal value of preferred habitat cover for species landscape‐level abundance is a surprisingly low value (10% within a 4 km buffer); landscape metrics exert different effects on species occupancy and abundance in high versus low population density regions of our study area. We conclude that knowledge of regional variation in landscape utilization will be essential for successful conservation of the species. The results also support the view that large‐scale landscape variables have high predictive power in explaining species abundance. Our study demonstrates the complex response of species occurrence at different levels of population configuration on landscape structure. The study also highlights the need for data in large spatial scale to increase the precision of biodiversity mapping and prediction of future trends.


| INTRODUCTION
Animals use and select habitat at multiple hierarchical levels and at different spatial scales within each level (McGarigal et al., 2016;Wiens, 1989). At the local level, site occupancy is determined by local habitat and community characteristics (Dunning et al., 1992;Van Buskirk, 2005). At the broader level, processes such as movement of individuals between local populations, regional patterns in landscape structure, or species interactions also affect occupancy patterns (Boscolo & Metzger, 2011;Dunning et al., 1992). Consequently, the spatial scale at which habitats affect the species differ at different levels that result from individuals behaviors in space, such as home ranges, population patches, or ranges (McGarigal et al., 2016).
The hierarchical nature and scale dependency of habitat use complicates evaluation of the habitat characteristics required to support populations. For example, the commonly used approach to determine conservation or management units from the habitat requirements of individuals may not be appropriate tool for securing population persistence (Jackson & Fahrig, 2015). Rather, knowledge of the spatial distribution of habitats, populations, population subunits, and individuals is a key element for the success of conservation and management programs. However, there is currently little understanding of how processes at different scales simultaneously influence species occupancy patterns (Du Toit, 2010;McGarigal et al., 2016;Wheatley et al., 2005).
Arboreal mammals are an animal group that have very explicit habitat and dispersal limitations and thus are interesting species to study scale-dependent habitat utilization. They depend on forest habitat, which is often heavily managed. Many of these species are reluctant to cross large forest gaps (Bakker & Van Vuren, 2004;van der Ree et al., 2010) and may, therefore, be unable to colonize suitable and empty habitat patches (but see Selonen & Hanski, 2004;Fey et al., 2016). Consequently, the distribution and abundance patterns of arboreal mammals are frequently influenced by habitat composition and configuration at both the patch and landscape level (Hurme et al., 2008;Mortelliti et al., 2011;Nupp & Swihart, 2000;Pardini et al., 2005). One such species is the Siberian flying squirrel (Pteromys volans), a forest-dwelling rodent with a declining population trend (Hanski, 2006;Jokinen et al., 2015;Santangeli et al., 2013a;Selonen et al., 2010). Habitat selection of the flying squirrel is well studied at the site level (i.e., selection of nesting and feeding sites; Hanski, 1998;Reunanen et al., 2002), but also at the landscape level (Hurme et al., 2008;Reunanen et al., 2000;Santangeli et al., 2013b). However, as is the case for other arboreal rodents, we lack knowledge about scaledependent landscape effects on population configuration. The species requires nest cavities that are usually found in old aspens (Populus tremula) in Norway spruce (Picea abies)-dominated forests (Hanski, 1998;Reunanen et al., 2002). The distribution of aspen and spruce is in turn correlated with the presence of fertile soils, which are frequently used for agricultural cultivation. Thus, and somewhat surprisingly, an earlier study found a positive correlation between flying squirrel occurrence and the presence of fields in the landscape (Santangeli et al., 2013b). In the long term, however, agricultural land use creates large areas of unsuitable open habitat in landscapes, which, depending on the dispersal potential of the species, negatively affects the population distribution of arboreal species (Van Apeldoorn et al., 1994). In other words, the response of flying squirrels to certain landscape variables, such as agricultural areas, should be nonlinear.
The objective of this study was to examine scale dependency and nonlinear effects of the landscape characteristics that affect occupancy patterns of flying squirrels. We expect habitat variables to affect at different spatial scales, when studied at different levels of population configuration (compare to multilevel and multiscale definition by McGarigal et al., 2016). In our case, the lowest level is (1) the occupancy of study site (hereafter site occupancy). At this level, we expect landscape variables to dominate at spatial scales that describe the home-range utilization of individuals. At the broader level, (2) the abundance of occupied sites in the landscape (hereafter landscapelevel abundance) may also be related to processes that are affected by landscape variables that operate at larger spatial scale than home range (like movement of individuals). We define the landscape-level abundance to depend on the configuration of functional subpopulation patches. Subpopulation patches are the areas where a species' local average abundance is higher than in surroundings. Finally, we predict that large-level (3) regional patterns in population density within the studied distribution range (hereafter regional density) may affect landscape utilization at site, and landscape levels (levels 1 and 2 above), that is, the responses may vary in areas of high and low population density (see e.g., Mazerolle & Villard, 1999;Sundell et al., 2012).

| Study species and study area
The flying squirrel is a small arboreal rodent ( Fig. 1) distributed from eastern Siberia and Japan to Europe, where the species is found only in Finland, Estonia, and Russia (Ognev, 1966;Wilson & Reeder, 2005).
In Finland, flying squirrels occur from southern Lapland to the southernmost coast of the country. The species is mainly associated with mature and old spruce-dominated mixed forests containing deciduous trees for foraging and with cavities for nesting and roosting (Hanski et al., 2000a;Santangeli et al., 2013b). Adults are site-tenacious, with home ranges averaging 8 ha for females and 60 ha for males (Hanski et al., 2000a;. Mean natal dispersal distances are between 1 and 2 km, with a maximum of almost 9 km (Hanski & Selonen, 2009

| Data collection
We used presence-absence data for flying squirrels obtained from 9 ha study site (300 × 300 m; n = 10,032) surveyed in [2003][2004][2005]. The survey was originally organized by Finnish Ministry of Environment (Hanski, 2006) to get knowledge on population status of the flying squirrel in Finland, because the species is protected by EU habitats directive. Size of the study site was selected to match average female home range (8.3 ha, Hanski et al., 2000b). The study sites were placed randomly in forest habitats with an average density of 10 sites per 200 km² (10 plots in every second 10 × 10 km square; for more information, see Santangeli et al., 2013b). The sites were located on mineral soil and were at least 1 km apart. They were inspected for signs of flying squirrel activity once by searching extensively for the presence of fecal droppings. This is a commonly adopted protocol for surveying flying squirrels (e.g., Hurme et al., 2005Hurme et al., , 2008Mönkkönen et al., 1997;Reunanen et al., 2002). All surveys were carried out between April and June by trained surveyors. Pellets are relatively easy to find due to their yellow color and deposition location: usually at the bottom of large aspen and spruce trees. Therefore, detection probably can be assumed to be close to one (Santangeli et al., 2013b) and the risk of false absences to be low in our data (Hurme et al., 2008). Thus, we have no reason to expect detection efficiency to be uneven within the surveyed area and do not anticipate spatial bias in the results.

| Analysis
We distinguished three levels in the information of species presenceabsence: site occupancy-present or absent (0 or 1), landscape-level abundance-proportion of occupied sites (0-1) within the range of spatial dependence (i.e., autocorrelation), and regional densityglobal spatial trend of occupancy proportion (0-1; Fig. 2). We build two models: In the first model response variable was the site occupancy. In the second model, response variable was the landscape-level abundance. The study site was, however, the unit of analysis in both models, because for landscape-level abundance, we used site-specific values based on components of universal kriging (local weighted mean). When analyzing site occupancy, we included landscape-level abundance and regional density as covariates in the model alongside other landscape variables. For the landscape-level abundance model, regional density was included as a covariate. Regional density based on trend surface component of the universal kriging. It was not modeled separately, but was used to analyze interactive effects in site occupancy and landscape-level abundance models.
The step-by-step details of the data management and statistical analysis are explained subsequently.

| Estimation of regional density and landscapelevel abundance with the geostatistical method universal kriging
That accounts spatial dependence of a variable as two components: large-scale trend surface interpolation (regional effects), and local structures of spatial autocorrelation (landscape patchiness) (Lam, 1983;Legendre & Fortin, 1989). The two-scale procedure is about 40% more effective for describing spatial patterns than ordinary F I G U R E 2 Principle of decomposing species presence-absence information: (a) population spatial pattern and sampling (arrangement simplified on the chart); (b) three additive components of species presence (i, ii, iii) projected on the diagram using one dimension of the study area. Note that two-dimensional models were used on the real data. Population patches and gaps are defined as areas of higher (+) and lower (−) local abundance in relation to the global trend of regional density. In the analysis, gross values of site occupancy (i + ii + iii) and landscape-level abundance (ii + iii) were explained by landscape metrics, while the partial components of landscape-level abundance (ii) and regional density (iii) were used as covariates accounting for the effects of lower level patterns kriging, which does not consider the large-scale trend (Ghiasi & Nafisi, 2015).
First, we estimated a two-dimensional quadratic trend surface model to explain species distribution density as a global gradient (twoway binomial regression). This was the regional density, that is, the largest level (level 3 in introduction), we used to describe the population configuration of flying squirrels. Second, using the residual values from the trend surface, a spherical semivariogram was modeled to describe the spatially autocorrelated portion of the population density structure in landscape level (Fig. S1). The subsequent weights of spatial dependence between sampling sites were based on the semivariogram model. Third, we estimated distance-weighted average occupancy of surrounding sites within the semivariogram range for each site. This was the landscape-level abundance of the species, that is the medium level (level 2 in introduction), used to describe population configuration of the flying squirrel. In other words, the landscape-level abundance described average occupancy rate at the scale of population patches estimated by the kriging method.

| Selection of explanatory landscape variables
We used five variables that are known or expected to correlate well with habitat and landscape utilization of flying squirrels.
Among the radii, the most parsimonious one was used in the final models, selected in the third step of the analysis (2.3.3; same for the rest of the variables). Habitat data from the same year as sampling of flying squirrels were based on the multisource national forest inventory (MS-NFI) in Finland (www.paikkatietoikkuna.fi), which employs satellite images, digital maps, and field measurements for producing forest estimates in the form of thematic maps (Tomppo et al., 2008). The thematic maps include, among others, volumes of main tree species: pine, spruce, birches, and other broadleaved trees.

2.
Aspen and alders, key tree species related to nest site presence and foraging conditions (Hanski et al., 2000a;. Data for these trees are not presented in MS-NFI maps, but we estimated the growing stock volume (m³/ha) of these species within a 1 km² that contains each study site by combining NFI field measurement plots and MS-NFI thematic forest maps. Based on volumes of aspen, grey alder, and black alder calculated for each NFI field plot (Tomppo et al., 2008), we created thematic map layers presenting the volumes of aspen, grey alder, and black alder. These map layers were produced by dividing the pixel values of the map layer "volume of broadleaved trees" into strata representing aspen, grey alder, and black alder, based on geostatistical interpolation of the volume proportions of the respective tree species in NFI plots.
Inverse distance weighting of the geographically nearest 12 NFI field plots was used in the geostatistical interpolation. The outcome of the interpolation was volume maps with similar resolution as the original MS-NFI thematic maps. Based on these maps, the volumes of the three tree species were then calculated for 1 km² grid elements, forming a uniform map grid covering the entire study region.

3.
Soil fertility (rank of 1, high to 8, low fertility), representing the longterm suitability of land for deciduous tree species and spruce, as well as for agriculture. We measured the average site fertility index at six radii around the centers of study sites: 100 m, 250 m, 500 m, 1 km, 2 km, and 4 km. The fertility index based on ranks using habitat type as a proxy for fertility: 1. Herb-rich sites, 2. Herb-rich heath forests, 3. Mesic, 4. Subxeric forests, 5. Xeric forests, 6. Barren forests, 7. Rocky and sandy soils, and 8. Summit and fjeld forest. This is commonly used proxy to describe soil fertility in Finland (Tomppo et al., 2008).

4.
Forest cover represents landscape connectivity for the species that prefer to not cross open habitats. We measured the proportion of the area covered by forest at eight radii around centers of study sites: 100 m, 250 m, 500 m, 1 km, 2 km, 4 km, 10 km, and 20 km.
Forest cover within radii of 10 km and 20 km was measured from CORINE land cover map data (European Environment Agency; proportion measurements for other radii were taken from MS-NFI data), as they might reflect very large-scale effects on the dispersal ability of the species, although these radii were not available for landscape variables (i), (ii), and (iii).

5.
Agricultural fields represent an expected negative coeffect at otherwise favorable fertile soils. We measured the proportion of fields at eight radii around centers of study sites: 100 m, 250 m, 500 m, 1 km, 2 km, 4 km, 10 km, and 20 km. Agricultural field cover within the 10 km and 20 km radii was calculated from CORINE land cover map data.

| Univariate models for every landscape variable for site occupancy and landscapelevel abundance
We used nonparametric, nonlinear generalized additive models (GAM; Hastie & Tibshirani, 1986;Li & Wang, 2013). Among several statistical methods used for studying species habitat selection, GAMs have a good optimum of flexibility for input data types, computational complexity, and popularity in the scientific community (Li & Wang, 2013).
We compared the effects of landscape variables at all measured radii in models for site occupancy and landscape-level abundance. The measurements of landscape variables within different radii are intercorrelated (Table S1). Thus, for further analysis (multivariate GAM), we selected the radius at which the correlation of each landscape variable was strongest with site occupancy and the landscape-level abundance of the species (DeCesare et al., 2012;Wheatley, 2010). Selection was based on the Akaike Information Criterion (AIC, Johnson & Omland, 2004;Anderson, 2008), which measures model parsimony. As a result, we maximized the ecologically relevant information that was entered into the final step of the analysis (multivariate GAM).

| Multivariate GAM
We used the landscape variables at the radii selected based on the univariate models described above. In the multivariate model (GAM) selection procedure, all possible combinations of the explanatory variables were compared along with their interactions with regional density of flying squirrels. The interactive terms were used to test for a difference in habitat use between sites with low and high regional density and to account geographically variable landscape effects.
Selection of the best set of explanatory variables was based on the AIC.
All explanatory variables and covariates were continuous variables.
The only categorical variable in the models was the response variable site occupancy, measured as presence-absence (binomial variable).
As variation in the explained variables was limited between 0 and 1, a logit link-function was used to generate normal residual distributions, and to ensure that all fitted values lay between 0-1. In the site occupancy models, variogram-based landscape-level abundance was used as a covariate covering the effect of spatial autocorrelation.
In the landscape-level abundance model, we incorporated spherical spatial correlation structure in the GAM to account the effect of spatial autocorrelation. To avoid over-fitted models, the degrees of freedom (df) for smoothing splines was limited to df ≤ 2. R packages ape (Paradis et al., 2004), gstat (Pebesma, 2004), mgcv , and sp (Bivand et al., 2013;Pebesma & Bivand, 2005) were used.

| Site occupancy model
In total, 10,032 plots in forest habitats were studied for flying squirrel occupancy. About 10.3% of the plots were occupied (1,030 plots, Fig. 3). According to the selection of optimal spatial scale from univariate models, the strongest correlations between site occupancy and landscape variables occurred with the proportion of preferred mature spruce-deciduous forest habitat within the 250 m radius, followed by three metrics within the 500 m radius: proportion of forest cover, average soil fertility index, and proportion of fields (see the Δ AIC 1 and 2 in the Table 1). That is, the relationships between 9 ha site occupancy and landscape variables were strongest within radii that represent a several times larger area around a study site than the plot itself (250 m radius = 19.63 ha, 500 m radius = 78.54 ha).
In the multivariate model selection, two of the highest ranked candidate models gained similar explanatory power (Δ AIC = 1.9, Table 2a). These two models included the four variables listed above along with stock volumes of grey and black alder, but the importance of the latter was low (  Table 2a and Fig. 4a): The effect of average soil fertility index increased toward higher fertility that is low index values; the proportion of fields had a positive effect up to 30% of landscape composition; the proportion of forest cover had a negative effect at higher than 50% of landscape composition; and a positive effect of preferred habitat was apparent up to 35% of landscape composition (Fig. 4a).

| Landscape-level abundance model
The semivariogram analysis revealed evident patchy landscape-level abundance of flying squirrel occurrence (R 2 = 0.17; Fig. 3). The average size of patchy population structures was approximately 32.5 km (semivariogram range), and according to the spherical semivariogram model, spatial autocorrelation of neighboring sites was halved at a distance of 11.3 km (Fig. S1). Comparison of univariate models of different spatial scales showed that landscape-level abundance was most strongly correlated with the proportion of preferred habitat, followed by soil fertility (see the Δ AIC 1 and 2 in Table 1). Both of these variables correlated most strongly with landscape-level abundance at the largest radii studied (4 km radius).
This suggests that the true optima probably lie beyond the range of studied scales. The next most important variables were forest cover within the 4 km radius and proportion of fields within the 1 km radius (Table 1).
The best multivariate model for landscape-level abundance (R 2 = 0.44, Tables 2b and 3b) included all variables except grey alder, which was included in the second highest ranked model (full model; Table 2b). While the importance of other variables over all candidate models was high, grey alder had relatively low importance (Table 2b).
In particular, the proportion of preferred habitat and the average soil fertility index had strong effects on landscape-level abundance (Table 3b). The effects of the four most important landscape metrics according to the best multivariate model are summarized as follows ( Fig. 4b): Average soil fertility index had a sharp optimum between average rank of 3.0-3.5; the proportion of forest cover showed an optimum between 50%-60% landscape composition; the proportion of agricultural fields showed a slight negative effect below 30% of landscape composition; and the proportion of preferred mature sprucedeciduous stands showed a relatively narrow optimum between 10%-15% of landscape composition.

| Differences at high and low regional density
Based on two-dimensional trend surface analyses, the regional density of flying squirrels was higher in western Finland and decreased Rank of the model Importance  (Fig. 3). The regional density, that is the global trend surface value, formed significant interactions with proportion of fields, average soil fertility index, and stock volume of grey alder in the site occupancy model (Table 3a) and with all variables except the proportion of forest cover in the landscape-level abundance model (Table 3b). That is, the positive effect of agricultural fields on site occupancy appeared strong in the high-density region (western Finland, up to 35% of landscape composition), but was absent in the low-density region (eastern and northern Finland; Fig. 5a). The site occupancy optimum occurred in slightly more fertile sites in the lowdensity region (around 2.7) compared with the high-density region (around 3.0). The effect of grey alder stock volume on site occupancy was always linear, but was strongly positive in the high-density region and slightly negative in the low-density region. However, this pattern was driven by a very small proportion of samples, as alder stock volume remained very low in the great majority of sites.
In the landscape-level abundance model, the effect of preferred habitat formed by far the strongest interaction with regional density (see the β-values in the Table 3b). The preference for preferred habitat was significantly higher in the low-density than in the high-density region (Fig. 5b). Despite this, there was only small difference in the optimum value for the preferred habitat: around 10% of landscape composition in high-density area and around 13% in the low-density area (Fig. 5b). For other variables, despite the significant effects of interactions with regional density were accompanied by low-effect sizes.

| DISCUSSION
In this study, we distinguish between scale-dependent effects of habitat characteristics on site occupancy and landscape-level abundance of flying squirrels, as well as explained geographic differences in the landscape effects. We are unaware of similar attempts to explain the partial components of presence of a species using landscape metrics.
The main results of our study are as follows: Landscape variables predicted the placement of population patches (landscape-level model) at least twice as well as they predicted the occupancy of particular sites (site occupancy model); the clear optimal value of preferred habitat cover for species landscape-level abundance is a surprisingly low value (10% within a 4 km buffer); landscape metrics exert different effects on species occupancy and abundance in high versus low population density regions of our study area.

| Site and landscape components of species presence
Our models clearly predicted placement of population patches better than occupancy of particular sites, that is, the model for landscape-level abundance produced a considerably better fit than that considering site occupancy. At the site level, extinction-recolonization dynamics of territories (sites) generate noise in the data because even high-quality sites can be temporarily unoccupied. The spatial scale that the landscape T A B L E 3 The best multivariate models (n = 10,032) of (a) site occupancy (model: df = 10,015.71, R 2 = 0.21) and (b) landscape-level abundance of flying squirrels (model: df = 10,013.35, R 2 = 0.44). See Fig. 4 for visual presentation of the modeled effects. The variable df values reflect whether an effect is linear (df = 1.0) or nonlinear (df > 1.0). To avoid over fitting of the models, limit was set at df ≤ 2 Main effect Interaction with regional density Stock of black alder within 1 km² 1.0 F = 16.91 <.001 −0.77 ± 0.14 t = −5.49 <.001 Species regional density 1.0 t = −1.62 .106 variables operated on, differed between site and landscape models in the manner we predicted (see also Fletcher et al., 2016;Jackson & Fahrig, 2015;Michael et al., 2016), probably because population-level processes should affect landscape-level abundance more than site occupancy. Landscape-level abundance was predicted well by the landscape variables measured at large spatial scales, in some cases at scales larger than we could measure in this study. Landscape-scale features are similarly more important than stand-scale variables in explaining species diversity in birds (Mitchell et al., 2001(Mitchell et al., , 2006. While many earlier mammal studies concentrated on species abundance at small spatial scales, for example in relation to habitat patch characteristics (Bowers & Matter, 1997;Mortelliti et al., 2011), our study supports the view that such studies would benefit from considering large-scale landscape features when aiming to predict the abundance of species (see also DeCesare et al., 2012;Lindman et al., 2015;Wheatley et al., 2005).
Conclusions about species responses to site and landscape characteristics clearly depend on the habitat measures considered (e.g., Schindler et al., 2013). For example, Michael et al. (2016) concluded that F I G U R E 4 Partial effects of four landscape variables in multivariate GAM on (a) occupancy of 9 ha sites and (b) landscape-level abundance (patches of 32.5 km radius). Gray areas represents 95% confidence intervals; line width represents sampling density. See details on model parameters in Table 3 local vegetation structure was more important than forest cover in the surrounding landscape in determining occupancy of 14 reptile species.
Unfortunately, we were not able to evaluate site-level variation in the nest cavity or food resource availability for flying squirrel. Earlier studies indicate that both large aspens providing cavities and alder trees providing food resources are important determinants of occupancy and fitness of flying squirrels .
However, aspen stock volume did not influence flying squirrel site occupancy in the current study, and the effect of alder stock volume was also small. This was an unexpected result, but suggests that the effect of these trees on site occupancy can be detected only at finer scales than we were able to measure in the current study (tree species at 1 km² correlated with flying squirrel presence at 9 ha scale). Instead, both alder and aspen had effects, albeit minor, in the landscape-level abundance model.

| Importance of agricultural fields and movement ability of flying squirrels
Unsurprisingly, the previously observed positive relationship between agricultural areas and flying squirrel occupancy (Santangeli et al., 2013b) weakened (high population density region) or disappeared (low population density region) when soil fertility was accounted for in this analysis. It is clear that fields as such are not a required resource for an arboreal animal like the flying squirrel, although the extent of agricultural land is correlated with soil fertility and the presence of edge habitat (Santangeli et al., 2013b). Nevertheless, it is interesting that even after accounting for soil fertility, fields appeared to have a positive effect on site occupancy in the high-density region (see discussion below). Furthermore, the landscape-level abundance of flying squirrels did not decline even if there were 30%-80% of fields in the landscape. This supports earlier studies indicating that the occupancy patterns and movement ability of flying squirrels may not be limited by a lack of structural connectedness even in highly fragmented landscapes (Selonen & Hanski, 2004;. Nevertheless, the population dynamics of the species are highly dependent on migration patterns (Brommer et al., 2017), which may also be related to the observed patchy density pattern that is characteristic to metapopulations (Hurme et al., 2008).

| Habitat threshold
In his review of bird and mammal occupancy patterns, Andrén (1994) concluded that the threshold where occupancy of a species starts to decline sharply lies somewhere between 10% and 30% landscape composition of the habitat (see also Swift & Hannon, 2010). For flying squirrel, landscape-level abundance declined sharply when the proportion of preferred habitat was below 10%-15% within a 4 km buffer (see also Reunanen et al., 2004). However, the abundance of flying squirrel also declined when the proportion of preferred habitat increased above the optimal level. One explanation for this unexpected pattern could be flexible habitat use and a preference for edge habitats (see below), but also good dispersal capacity that allows the species to inhabit fragmented landscapes (Brommer et al., 2017;Selonen & Hanski, 2004. Critical habitat thresholds are likely to be species-specific and our understanding of the mechanism underlying them remains poor, especially for mammals at landscape scales, as studies have concentrated on other F I G U R E 5 The two strongest interactive effects of regional density on the effect of (a) agricultural fields on site occupancy and (b) preferred habitat (mature spruce and deciduous stands) on landscape-level abundance of flying squirrel. The partial effects in the best multivariate GAM are presented for the lowest third of population regional density (blue), and the highest third of population density (red) according to the global trend surface model. The line width represents sampling density, and the gray area represents the 95% confidence interval. The darker gray area in part (a) represents overlap of the 95% confidence intervals of high and low population density models. See the model details in the Table 3 (a) (b) Landscape-level abundance species groups such as birds (Swift & Hannon, 2010; but see e.g., Pardini et al., 2010). Thresholds clearly depend on study scale (Swift & Hannon, 2010), and in our case, the site occupancy peaked when the proportion of preferred habitat (250 m buffer) was around 40%.

| Habitat use in low versus high-density regions
We found that landscape variables had different effects in highversus low-density regions of our study area. For landscape-level abundance, preferred habitat had higher explanatory power in the low than in the high-density region in western Finland. This was a surprising result and indicates that further study on flying squirrel habitat preferences in different regions may be needed. In the west (the high-density region), the species preferred sites near fields and with more grey alder, but this was not the case in the low-density region.
Instead, soil fertility had a higher optimum, and landscapes with black alder were more preferred in the low compared with high-density region. Some of these differences were very modest, but they nonetheless indicate differences in habitat use between regions. For example, clay soil type is more common in west than in northeast Finland (edge habitat between fields and forest in clay soil is suitable habitat for grey alder), whereas shoreland forests near lakes (suitable habitat for black alder) may be more typical habitat for flying squirrels in particular areas in eastern Finland. Alternatively, our results could indicate density-dependent habitat use, such that some individuals are forced to use low-quality habitat when population density is high (Gill et al., 2001;Rozenzweig, 1991;Sundell et al., 2012). However, this hypothesis is not supported by the fact that the response to preferred habitat was not density dependent in the site occupancy model.
In this study, we demonstrated that landscape variables have different effects when occupancy patterns are described at different population levels (site occupancy vs. landscape-level abundance). One factor that we did not control for was the effect of predators on site occupancy. However, earlier analysis of these (F. G. Blanchet et al. unpublished) and other data (Selonen et al., 2010) indicate that predators do not play a major role in explaining flying squirrel density in the landscape, although it is clear that predators affect individuals locally.
We are also unaware of any diseases or parasites that might affect flying squirrel densities. We did not control for historic land use, but the whole study area is under intensive forest management, and the age structure of forests used in the current analysis describes changes in forest structure during the last 50 years. Land use change from forest to agricultural areas has been minor during recent decades in Finland. Thus, in recent history, there have not been major changes in agricultural land area, although the total cropland area has slightly increased in the country (Greenhouse Gas Emission, 2016). Land use in the form of forest cutting during historic times (about >150 years ago) was very intense in western and southwestern Finland (Gyldén, 1850), and historic records from that time indicate that flying squirrels were more abundant in eastern than in western Finland (Mela, 1882).
Interestingly, the current situation is the opposite, that is, population density is highest in western Finland, indicating no negative effects from historic forest use on current flying squirrel density.

| Conservation implications
A large majority of the study sites occupied in this study are found in landscapes where the proportion of preferred habitat is below the observed optimal value for site occupancy and landscape-level abundance. This highlights a potential threat to the future survival of the population, in particular in the low-density region in the northeastern distribution range of flying squirrels in Finland. As the species' population dynamics considerably depend on immigration (Brommer et al., 2017), it is possible that there is a threshold level for the landscape proportion of preferred habitat at which immigration is no longer able to maintain population stability. In this kind of situation, sudden extinctions are possible even over large areas, as has been observed for the Greater glider (Petauroides volans) in Australia (Lindenmayer et al., 2011). Based on our results, the optimal landscape for site occupancy by the flying squirrel is a mosaic containing around 35%-40% of mature spruce and deciduous forests situated in semi-forested landscapes on fertile soils. The optimal size of landscape grain for site occupancy appeared to be in a radius of 250-500 m. This is a larger area than the average female homerange size (the 100% minimum convex polygon is around 8 ha, i.e., a 160 m radius, Selonen et al., 2001) or current management units used for flying squirrel conservation (radius 10-30 m; Jokinen et al., 2015). The optimal area of preferred habitat within 250 m buffers was around 2 ha, but our analysis cannot be used to evaluate the minimum forest area needed at forest cutting sites containing nest sites of flying squirrels (EU Habitats directive). Instead, we suggest that efficient species conservation follows ecological patterns and should cover, in addition to protection of high-quality habitat in occupied sites, also its near surroundings, landscape connectivity, protection of temporarily unoccupied sites, and other population-level effects. For this purpose, conservation might benefit if management considered the scale of functional subpopulations, which based on our results could be an area of at least 1,000 km²-the approximate size of a population patch.

| CONCLUSION
Our result supports the view that large-scale landscape variables have high power to predict species abundance (Mitchell et al., 2001(Mitchell et al., , 2006.
We also conclude that the optimal management plan for species should take into account spatial scales that capture population-level processes and not only scales that describe the site preferences of individuals. In addition, knowledge of possible regional variation in landscape utilization or density-dependent habitat use will be essential for the success of a species management or conservation program. This is important even when a species appears to specialize on a certain habitat, as in our case does the flying squirrel, a gliding mammal that requires forest habitat with specific key elements providing food and nesting cavities. We suggest that inclusion of landscape-level population patterns can significantly increase the precision of biodiversity mapping and prediction of future population trends.