Effects of forest management on the spatial distribution of the willow tit ( Poecile montanus )

Modification, fragmentation and loss of boreal forest habitats have been intensive during the last century due to forestry practises and land use. This has been related to population declines of many forest species, yet the mechanisms affecting on the background are largely unknown. The willow tit, a primary cavity-nesting species that was once the 4 th most common bird species in Finland is nowadays endangered. Earlier findings suggest that the willow tit population is affected by the reduction of nesting sites, decaying snags in forests and the loss of mature forests which contain the food storages during the winter. In this study we are searching for the mechanisms how the forest management methods could explain the decline of the willow tit population. We used long-term breeding data of the willow tit nesting sites from 1990 to 2020 collected in a study area in northern Finland to analyse if forest management affected nearest neighbour distances and natal dispersal and breeding dispersal distances. We used Geographic Information System (GIS) methods to combine the ecological breeding data to accurate spatial forest management and habitat quality data. The data was analysed with linear mixed models. We found that clear-cuttings affected the willow tit dispersal and neighbouring nest distances more than thinnings. Both clear-cuttings and thinnings increased the nearest neighbour distances. The natal and breeding dispersal distances were lengthened by increasing proportions of clear-cuttings. The habitat loss caused by clear-cuttings and the decrease in habitat quality caused by thinnings has had a major role in the decline of the willow tit population. The forest management actions were estimated to explain about 65 % of the willow tit breeding density decrease in the study area. The effects of forest management were witnessed in a cumulative 0 – 30-year period meaning that forest management causes long-term habitat degradation and loss. Availability of deciduous snags in the forests can compensate the habitat loss to some extent by providing better breeding opportunities. As the effects of clear-cutting were more severe to the willow tit than thinning, we recommend using other forest management methods than clear-cutting as the main management method.


Introduction
Anthropogenic impact on environment (e.g., forestry, urbanization, land use and agriculture) has caused habitat modification, loss, and fragmentation worldwide (Sala et al., 2000;Newbold et al., 2015).According to IPBES (2019) report, 75 % of terrestrial areas on the globe have been significantly modified by human actions.The habitat loss is a major problem because it reduces the availability of e.g., breeding sites and wintering territories (Schmiegelow and Mönkkönen, 2002;Määttänen et al., 2022).Decreasing breeding success and survival (Lampila et al., 2006) may eventually lead to endangered populations, extinctions, and biodiversity loss (Newbold et al., 2015;Määttänen et al., 2022).The effects of habitat loss are rarely immediate, rather the change is usually seen as gradually declining populations in a long-term period (Hanski, 2011).
Forests cover approximately one-third of all land areas (FAO, 2012) and currently global forest area is estimated to be about 68 % of the size compared to pre-industrial times (IPBES, 2019).The modification of boreal forests due to forestry practises has been intensive during the last century (Edwards et al., 2022).About two-thirds of the boreal forests are under management (Gauthier et al., 2015).In Finland, the use of forests was intensified in the beginning of 20 th century, when even-aged forest management practices took place as the main forest management method (Leikola, 1987).The even-aged forestry method is based on forest cultivation, where mature forests are cut down, renewed by planting saplings or seeds, and managed by removing small-diameter trees i.e., thinning the forest from below regularly during the growing season (Kuuluvainen et al., 2012).
The even-aged forest management has reduced the number of oldgrowth forests (Shorohova et al., 2011) and decaying wood (Hanski, 2005) and increased the number of open clear-cuttings, sapling areas, and young forests (Gauthier et al., 2015).These actions have caused habitat modification, as well as habitat loss and fragmentation of oldgrowth forests (Määttänen et al., 2022).In addition, thinning the forest from below affects the biodiversity by reducing the number of niches, foraging and breeding opportunities and cover from predators e. g., for small passerine species (Eggers and Low, 2014).In recent years there has been an effort to develop forest management methods into a more sustainable and versatile direction in conserving biodiversity.By studying the effects of forest management on natural populations we can support this development towards more sustainable direction.
Intensive forest management methods greatly affect cavity-nesting bird species (Virkkala, 2004;Fraixedas et al., 2015).The forest management has been proposed to have caused population declines of forest birds (Hyvärinen et al., 2019).The lack of decaying wood reduces the breeding opportunities (Vatka et al., 2014a) and the lack of large mature trees harnessing great numbers of invertebrates (Schowalter, 1995) as food during the winter can increase the stress hormone levels and decrease the survival of the willow tits wintering in young forests (Cīrule et al., 2017).However, the exact mechanisms of the forest management effects on cavity-nesting birds are still unknown.Birds are important for biodiversity, because they manage various ecological functions at different trophic levels, e.g., being predators and prey, controlling pests, and dispersing seeds (Lindbladh et al., 2017).In addition, primary cavity-nesting birds provide breeding opportunities for other species while excavating breeding holes for themselves (Lampila et al., 2006).
We use the willow tit (Poecile montanus), an endangered (Hyvärinen et al., 2019) primary cavity-nesting forest species (Orell et al., 1999) as a model species in this research.The willow tit was ranked the 4 th most common bird species in Finland in the 1940 ′ s and in the 1950 ′ s ( Merikallio, 1958).The population has considerably decreased in numbers during the last few decades (Fraixedas et al., 2015;Virkkala et al., 2020).It has been previously proposed that the population decrease of the willow tit is highly affected by even-aged forestry practises, in particular the reduction of decaying snags in forests due to the management actions (Vatka et al., 2014a;Fraixedas et al., 2015;Hyvärinen et al., 2019).
In this study, we concentrate on the effects of clear-cutting and thinning actions on the willow tit, rather than on observing direct habitat quality variables.Studying spatial breeding distribution in the scope of environmental change gives us different kind of perspective compared with e.g., nest site selection (Lewis et al., 2007;Vatka et al., 2014a;van de Loock et al., 2020), breeding success (Vatka et al., 2011) and survival studies (Lampila et al., 2006).Dispersal and nearest neighbour distances are rarely studied in this context.Remote sensing methods used to map forest management areas combined with a highquality long-term breeding data gives us new opportunities to discover the effects of forestry practices (Pimm et al., 2015).Combining ecological breeding data with spatial forest management data enables the study to be spatially and temporally large-scale.
The effects of environmental change can vary in different temporal and spatial scales (Levin, 1992;Gustafson, 1998;González-Megías et al., 2007).For example, the spatial distribution can be affected by the changes in nesting site vicinity (e.g., the condition of the chicks and therefore the ability to disperse further away) or in landscape scale (e.g., the possibility to find their own territory).The forest habitat recovers slowly from forest management actions.The habitat degeneration in long-term and in a large-scale can therefore be more relevant than a recently occurred habitat loss, for example due to a small clear-cutting in the nesting site vicinity.
The aim of the study is to research the effects of environmental change, caused by forest management methods (clear-cutting and thinning) to the willow tit.We focus on three spatial breeding distribution factors: the nearest neighbour distances, the natal dispersal and the breeding dispersal distances.The distances to the nearest neighbour nests reveal if the breeding density has been affected due to changes in the environment.Correspondingly, natal and breeding dispersal distances provide us information on the effects of forestry practises on the possibilities to find and occupy suitable breeding territories.We expect that both studied forest management actions will have negative effects on the willow tit population, that is, forestry practises will cause loss of suitable breeding habitats of the willow tits forcing them to disperse further and nest further from each other.However, clear-cutting is expected to have a stronger effect than thinning.

Study area
The ca. 25 km 2 large study area is located on the north-east side of the city of Oulu (65 • N, 25 • 30 ′ E), in northern Finland.The area consists of managed commercial forests with varying fragments of clear-cuttings, sapling stands, young and mature forests.Peatland forests are heavily drained.There are also a few mires, two lakes, two conservation areas and a small river crossing the area.The surroundings of the river are almost in natural state.Most common tree species are Scots pine (Pinus sylvestris), Norwegian spruce (Picea abies) and Birch (Betula spp.).About 4000 natural and induced snags of decaying trees (mostly birches) for primary cavity-nesting species are found all over the area.The induced snags are replacing fallen natural snags, so they should not bias the study by attracting birds to certain areas (Vatka et al., 2014a).Intensive forest management actions have been ongoing in the area since 1980 ′ s.Similar forest landscape continues outside the study area.

Study species
The willow tit is a small (weight 10-12 g, (Rytkönen et al., 1996)), resident, primary cavity-nesting passerine species (Orell et al., 1999), which is dependent on suitable decaying snags for nesting in their territories (Vatka et al., 2014a).The species distribution range is in the northern parts of Eurasian continent (Cramp et al., 1993).The territory of a willow tit is formed after the first life year and in many cases the birds stay in the same territory for their lifetime (Orell et al., 1999).The willow tits can occupy several types of forests as their territories (Orell and Ojanen, 1983a).The study population is an open population where individuals can recruit to and from surrounding areas (Lampila et al., 2006).
The willow tits nest in cavities excavated by themselves in decaying snags (Orell and Ojanen, 1983a;Orell and Koivula, 1988) which in the study area are mostly birch (Orell and Ojanen, 1983a).The willow tits do not require old-growth forests as the breeding territory (Cramp et al., 1993), if the habitat requirements are fulfilled: decaying snags for excavating nest cavities (Vatka et al., 2014a), moist habitat type (enables suitable decaying snags to form (Lewis et al., 2007)) and forested habitat, which provides shelter from the predators (Siffczyk et al., 2003) and improves foraging possibilities (Rytkönen and Krams, 2003).

Collecting the long-term breeding data
The long-term breeding data of the willow tit had been collected since 1975 (Orell and Ojanen, 1983a).The mapping of breeding sites started in early spring.Information of previous nesting locations, known territories and the locations of decaying snags were used in searching for the nests.Found active nesting sites were monitored regularly during the breeding season until the nestlings fledged (Orell et al., 1994;Orell and Belda, 2002).The nestlings were individually marked with numbered aluminium rings at 13 days old and breeding adults with unique combinations of aluminium and coloured plastic rings during the nestling period (Finnish Ringing Centre Licence number 180), which enabled the recordings of individual dispersal distances.Standard methods were applied (Orell andOjanen, 1983a, 1983b) for determining the age and sex (Laaksonen and Lehikoinen, 1976;Svensson, 1992) of the parent birds and the weight of the nestlings while ringing.

Spatial data
Environmental changes (in forestry, building and other land use) have been documented in the study area since 1970 ′ s.Lacking information of some forest management actions were filled with utilizing open-source aerial images, National Forest Inventory data, information provided by the Forest Management Association and interpreting notifications of forest use actions provided by the Finnish Forest Centre.We formed a spatial vector data by digitizing the surveyed environmental changes with the program QGIS (QGIS.org,2022) from the years 1961 to 2020.A total number of 672 forest use actions were digitized, of which 316 were clear-cuttings and 356 thinnings.Clear-cuttings were on average (±SD) 3.02 ± 3.50 ha in size and thinnings (±SD) 3.04 ± 3.53 ha in size.88 % of the forest management actions were defined for the exact year the forest was cut.For the remaining 12 % we estimated a time frame of maximum of 10 years when the cutting had been executed.For these management areas, we marked the cutting year to be the earliest year of our estimation.Most of these were older forest cuttings, dating back to 1980 ′ s or before that.We could not expand our study to a longer period because of the lack of information of the older forest cuttings.
We created 57 m, 113 m, 329 m, and 500 m radius sized buffers over each nesting site in the study years 1990-2020 with the 'buffer' function.The first three buffer sizes were determined by following Vatka et al., (2014a) study, where 1 ha buffer (57 m radius) was considered as the habitat in the nesting site vicinity, 4 ha (113 m radius) the foraging area (Rytkönen and Krams, 2003) and 34 ha (329 m buffer) the maximum winter territory size (Siffczyk et al., 2003).The 500 m radius buffer (77 ha) describes the landscape scale.It was the largest buffer size which could be formed without the buffers crossing over the study area borders more than 20 % on average.
The area of the forest management actions (ha) inside each buffer size (57 m, 113 m, 329 m, and 500 m radius) were calculated with the function 'difference' in QGIS.For each buffer around a nesting site, we built a variable presenting either clear-cuttings or thinnings inside the buffer that had been executed 0-4 years, 0-10 years, and 0-30 years before the nesting occurrence.The cumulative forest management age periods include different forest growth periods; 0-4 years after cutting includes the open stage, 0-10 years after cutting includes the open stage and sapling stage and 0-30 years after cutting includes the open stage, sapling stage and young forest stage.The several spatial and temporal scales were considered because forest management can have different effects in a short-term or in a long-term cumulative period and in different spatial scales.The information of the original buffers is summarized in Appendix A; Table A1.For standardizing the different buffer area sizes, we converted the same information into proportions of managed forests (%) for each buffer size.
We formed the number of snags data of the known locations of decaying snags in the study area, from which we had accurate data from 2005 onwards.This was done with 'Heatmap (Kernel Density -estimator)' analysis according to the coordinate points of snags in the study area to find out the areas that had plenty of snags.The analysis was done with a 500 m radius and pixel size 1, Kernel shape was defined uniform.Only the snags that were more than 5 m away from each other were included in the Heatmap analysis, to assure that the same snags were not included twice in the data.We used the 500 m radius because it enabled us to form a heatmap of the decaying snags covering over the whole study area.The number of snags data was divided into three categories, describing different decades during the research period.Snag coordinates from 2005 to 2006 represented the 1990 ′ s and early 2000 ′ s years, as was done in Vatka et al. (2014a).Coordinates from 2007 to 2015 and 2015 to 2020 represented their current time.The estimated number of snags around a nesting site was calculated with 'Zonal statistics' analysis separately for the different categories.
The willow tit breeding density is naturally higher for example in moist habitat types (Lewis et al., 2007) and they avoid non-forested areas as breeding sites (Vatka et al., 2014a).To account for this we formed two explanatory variables, the forest type gradient representing the forest vegetation eutrophy (Cajander, 1926) and the proportion of non-forested areas.Both variables also cover the forest management areas.The forest type gradient was calculated using the open-source data of forest habitat type 2017 from the National Forest Inventory provided by the National Resources Institute Finland.This was a raster data where habitat type richness value was given for each raster pixel from the categories 1) herb-rich forest to 6) barren heath forest (see Appendix A; Table A2 for detailed information).The scale indicates decreasing forest type richness.The forest type gradient value for each nesting site within 57 m, 113 m, 329 m, and 500 m buffer areas was calculated from the forest habitat type raster by using 'Zonal Statistics' tool.
The non-forested areas data was formed by using the main forest habitat type 2017 raster data from National Forest Inventory provided by the National Resources Institute Finland.The data includes six main forest habitat type categories which are presented in Appendix A; Table A3.First, we calculated the number of forest areas around each nesting site within 57 m, 113 m, 329 m, and 500 m buffer areas by choosing categories 1) mineral soil, 2) spruce mire, and 3) pine mire with the 'Zonal statistic' tool.This information thus includes clearcuttings and thinnings.The proportion of non-forested areas was calculated by extracting the proportion of potential forest areas from the total buffer areas.The non-forested areas are typically roads, settlements, fields, mines, and open mires.

Distance measures
We calculated the distance to the nearest neighbour nest with the function 'NNjoin' in QGIS.The function recognizes the nearest points in a point cloud, joins them together and calculates a distance between them in meters.
The coordinates of nesting sites were marked according to Finland Uniform Coordinate System (KKJ), where one coordinate unit is one meter (Ollikainen and Ollikainen, 2004).This enabled us to use the twodimensional Euclidean distance formula (Anderson, 1971) for counting the natal and breeding dispersal distances.In the breeding dispersal distance data, we took account only breeding parents that had been nesting in two subsequent years.
Because similar forest habitat continues outside the study area, the nearest neighbour distances and the dispersal distances are affected by the location of the nesting in the study area.We did not look for nesting sites outside the study area so the birds breeding near the study area borders can have their actual nearest neighbour nests outside the area.This can lengthen the apparent nearest neighbour distances.When measuring dispersal distances inside the study area, the potential maximum dispersal distance for birds breeding near the study area centre is the radius of the study area.For birds breeding near the study area borders the maximum dispersal distance is two times the study area radius.In addition, the birds can recruit outside the study area thus these individuals are lacking from the analysed data.This causes a systematic artifact where the apparent nearest neighbour and dispersal distances can be longer from nesting sites near the study area borders than near the central point of the study area (Pakanen et al., 2016).We accounted for this artifact by including the distance to the study area centre as a covariate in the analysed models.The central point of the study area was determined with the function 'mean coordinate(s)' in QGIS based on all the nesting sites during the study period.The distance between a nesting site and the central point was calculated with the function 'NNjoin'.

Statistical analyses
We used linear mixed models (LMM) to analyse the data of distance to the nearest neighbour and natal and breeding dispersal distances.LMM models are suitable for analysing complex linear regression models, where we have repeated measures of the same individual breeding birds in different breeding years (Harrison et al., 2018).Sample size for distance to the nearest neighbour was 2481 nests, for natal dispersal distance 681 recruits and for breeding dispersal distance 2468 individuals.
We chose several explanatory variables in our global models based on previous knowledge of the habitat qualities and individual features of parents and offspring that are known to affect the spatial breeding distribution of the willow tits (Table 1).By considering the confounding effects of these variables, we can better reveal the effects of forest management on the nearest neighbour and dispersal distances.The variables forest type gradient and non-forested areas were used to eliminate spatial autocorrelation in the distance to nearest neighbour and natal dispersal distance analysis.The nestling weight was used to describe the nestling condition which can affect the natal dispersal distance for example, the probability for the heavier nestlings to recruit is greater than with the lighter ones (Verhulst et al., 1997).Age and sex of the birds were included in the analysis because they are known to affect the dispersal distances and can affect the nearest neighbour distances due to the social hierarchy (Koivula et al., 1993).For example, the willow tit females are known to disperse further than males (Orell et al., 1999) and the older birds can be able to maintain a larger distance to their nearest neighbours.The year was used as a fixed explanatory variable to overrule the possible temporal trend caused by the decline in population size (for other reasons besides the forest management effects).In addition, we used the year as a random block variable to bind the effects of annual fluctuating circumstances.Female and male bird ID (ring number) were used as random variables for taking into account pseudoreplication in case the same bird had been nesting many times during the study period.Territory was used as a random variable for correcting the breeding site dependencies and differences in territory qualities.

Model selection
To determine the best temporal (x, Table 1) and spatial (y, Table 1) attributes for the clear-cutting and thinning variables, we built several global models that were structured according to the final global models presented in the Table 2, except that the buffer sizes (y) and cumulative forest management age categories (x) were fluctuating.The fit of the temporal and spatial attribute of the variables were evaluated by comparing AIC-values (Burnham and Anderson, 2004) of the global models.The final global model (Table 2) was selected based on the lowest AIC-value.In Appendix B; Table B1. is a detailed description of how the AIC comparison was done.
The final model selection was done with the model averaging procedure in the MuMin package (Barton, 2020).Here we used the final global models (Table 2).First, we used the function 'dredge' to fit all possible variable combinations of the final global model.Then we used the function 'get.models' to subset models with ΔAIC lower than 2 units that were considered as the fittest models.If more than one model was within 2 AIC units, we used model averaging with the function 'model.avg'.For comparing the effects of the forest management actions, we scaled the explanatory variables and analysed the models again with the same procedure.Results of the scaled fixed effects analysis are found in Appendix C.
Models were fitted with the function 'lmer' from the package lmerTest (Kuznetsova et al., 2017) in R (R core team, 2021).We evaluated the normality of residual distributions with QQ-plots.The response variables in distance to the nearest neighbour and breeding dispersal distance analysis were square root transformed to normalize the distributions.The forest type gradient variable was found to have minor non-linear effect in relation to the nearest neighbour and natal dispersal distances.This was considered by fitting the forest type gradient variable as a polynomial and linear variable in the models.The multicollinearity of the models was examined with function 'vif' from the package car (Fox and Weisberg, 2019), which is used to inspect variance inflation factors in a model.Variance inflation factors were lower than two units within all variables.
We used the Moran's I test to detect if the final global models had significant spatial autocorrelation (Dormann et al., 2007).First, we used the function 'correlog' from package ncf (Bjornstad, 2020) for visual inspection of the autocorrelation.Global Moran I test was run with function 'moran.test'from package spdep (Bivand and Wong, 2018) for testing the significance of autocorrelation by using half of the maximum

Table 1
The descriptions and units of response, explanatory and random variables used in the analysis of the global models ( distance as a scale.Significant autocorrelation was not found in the final global models.

Distance to the nearest neighbour
The average (±SD) distance to the nearest neighbour was 323 ± 147 m.Clear-cuttings and thinnings significantly increased the distances between the nearest neighbours (Table 3).The more there were cumulative 0-30 years old clear-cuttings or thinnings within 500 m from the nesting site, the longer the neighbouring distances were (Fig. 1a,b).The scaled estimates were 0.866 for clear-cuttings and 0.377 for thinnings (see Appendix C; Table C1.for the scaled fixed effects table), indicating that the effect of clear-cuttings was stronger than of thinnings.On the other hand, the number of snags significantly decreased the distance between the neighbours (Table 3), meaning that in an area where there were plenty snags, the birds were breeding closer to each other (Fig. 1c).Also, the nearest neighbour distances were increased by low forest vegetation eutrophy levels (i.e., linear forest type gradient).The polynomial forest type gradient effect reduced the nearest neighbour distances, and the proportion of non-forested areas increased the distances though these effects were not significant.
The distance to the study area centre increased the distances between the nearest neighbours, implying that birds living near the study area borders had a longer distance to their nearest neighbour.This artefact was expected because similar forest habitats continue outside the study area.Birds breeding near the study area borders could have their actual nearest neighbour outside the study area.Also, the older females had longer nest distances than the juveniles.
To demonstrate the effect of the nearest neighbour distance, we estimated the roles of clear-cutting and thinning on the changes in actual breeding densities in the 25 km 2 study area based on the results given by the distance to the nearest neighbour analysis.The estimation procedure is described in detail in Appendix D. We calculated the average breeding density to have declined about 43 % during the 1990-2020.The effects of clear-cutting and thinning on breeding density were estimated from the original data to the relationship between the distance to the nearest neighbours (Table 3, Appendix D; Fig. D1) and the breeding density (Appendix D: Fig. D2).Clear-cutting effect corresponds to a 19 % and Table 3 The model averaged (conditional average) results of the linear mixed model analysis explaining the willow tits nearest neighbour distances (square roottransformed) in 1990-2020.Female and male ID (ring number) and the breeding year were assigned as random variables in the analysis.3, by keeping the other variable effects on their mean values.
S. Kumpula et al. thinning effect to an 9 % decline in the breeding densities.Together the forest management actions caused about 28 % of the decline in the focal willow tit population density, which explains about 65 % of the total 43 % decrease in the breeding density.

Natal dispersal distance
The mean natal dispersal distance (±SD) was 2043 ± 1241 m for females and 1843 ± 1066 for males.Of the examined variables describing forest management methods and habitat only cumulative 0-30 years old clear-cuttings within 500 m from the nesting site had nearly significant effect on the natal dispersal distance (Table 4, Fig. 2).The natal dispersal distances were lengthened with the increasing proportion of forests clear-cutted 0-30 years before the nesting occasion.
As expected, the distance from the study area centre increased the natal dispersal distance, meaning that recruits near the research area border seem to disperse further than recruits that have their natal nests near the study area centre.This artifact arises from the fact that natal dispersal distance for birds dispersing near the study area borders is two times the study area radius at maximum whereas for birds dispersing near the study area centre the dispersal distance is the radius of the study area at maximum.In addition, in an open population some of the birds may also disperse outside the study area, thus they are lacking from the analysed data.We found also that the females had longer natal dispersal distances than the males.

Breeding dispersal distance
The breeding dispersal distance (±SD) was on average 215 ± 224 m for females and 224 ± 228 m for males.The breeding dispersal of parent birds was significantly longer when the proportion of cumulative 0-30 years old clear-cuttings within 57 m of the nesting site was increased (Table 5, Fig. 3).Also, the juvenile birds dispersed significantly further than the older ones.In addition, the sex of the bird entered the averaged model, though the effect was not significant.

Discussion
The forest management affected the spatial breeding distribution of the willow tit.Both clear-cuttings and thinnings increased the nearest

Table 4
The fixed effects results of linear mixed model analysis with function 'lmer' of the fittest model, explaining the willow tits natal dispersal distances in 1990-2020.Year and territory (correcting the effect of the same natal nests) were assigned as random variables in the analysis.4. The other variable effects were kept on their mean values.

Table 5
The model averaged (conditional average) results of the linear mixed model analysis explaining the willow tits breeding dispersal distance (square roottransformed) in 1990-2020.ID (ring number) and territory (correcting the male and female breeding site dependency) were assigned as random variables in the analysis.neighbour distances in a cumulative 0-30-year period before the nesting though the effect of clear-cuttings was stronger.Distance to the nearest neighbour refers to breeding density which was estimated to have declined about 65 % because of clear-cuttings and thinnings based on the nearest neighbour distance results.Only the clear-cuttings in a cumulative 0-30-year period affected the dispersal distances by lengthening the breeding dispersal distance in a 57 m spatial scale.The natal dispersal distances were increased by clear-cuttings in a 500 m buffer; however result was only nearly significant.
The results indicate that clear-cuttings have more severe effect on the spatial breeding distribution than thinnings.The clear-cutting of forests explains almost half of the breeding density decrease in the study area, so their effect has been remarkable.It seems that clear-cutting of forests cause habitat loss for the willow tits and the effect is built in a cumulative perspective where at least in a 0-30-year period before the nesting are increasing the nearest neighbour (i.e., breeding density) and dispersal distances.The willow tits' habitat requirements (Siffczyk et al., 2003;Lewis et al., 2007;Vatka et al., 2014a) are lost for several years when all or most of the trees are removed in the clear-cutting.Though the thinning of forests does not cause direct habitat loss, it can reduce the breeding opportunities for the willow tits, if all the smalldiameter deciduous (decaying or alive) trees are removed from the forest during the thinning.Eggers & Low (2014) discovered that the reduction in small-diameter spruce density in the forest understory decreased the nesting success and survival of the willow tits.Our results are in line with Betts et al. (2022) study where they found that long-term habitat degradation and loss caused by forest management has resulted in forest bird population declines.
With a close-relative to the willow tit, the black-capped chickadee (Poecile atricapillus), the breeding territories were found to be larger in managed habitats compared with unmanaged ones (Fort and Otter, 2004).Also, Siffczyk et al. (2003) discovered that the wintering territories of the willow tits were larger in forests under management.The willow tits breed inside their wintering territories (Haftorn, 1999) meaning that there probably is a link between the sizes of wintering territories and breeding territories.For compensating the loss of highquality habitats, the breeding willow tit pairs may have to increase the size of their breeding territories.This can in theory cause a part of the witnessed increase in the nearest neighbour distances and explain some of the effect of forest management to the breeding density.However, distinguishing suitable habitats from non-suitable ones can be difficult because habitat quality is likely a continuous variable.All forested areas are not the same as a suitable habitat.
The number of decaying snags is an important factor for willow tit habitat selection and habitat quality, which was found by Vatka et al.   (2014a).Our study supported this finding because the increasing numbers of snags decreased the distances between the nearest neighbours.The number of snags seemed to function as a compensating factor over the effects of forest management.The forest management practises tend to reduce the number of snags in the forests (Hanski, 2005) which can limit the number of high-quality breeding territories for the willow tits and can therefore define the possible number of breeding pairs and breeding territories in a forest area.In our study area we have attached fallen deciduous decaying snags in the nearby trees to improve the willow tit nesting possibilities.Because of this procedure the habitat quality has likely not decreased as much inside the study area compared to areas of other commercial forests.This makes our results on nearest neighbour and dispersal distances conservative, meaning that the forestry effects are probably even stronger outside the study area.
The lengthening effect of cumulative 0-30 years old clear-cuttings in natal and breeding dispersal distances point to difficulties in finding a suitable nesting location as we expected.In breeding dispersal distance this effect was witnessed in the nesting site vicinity and in natal dispersal distance in the landscape scale, though the natal dispersal distance result was only nearly significant.The effect of clear-cuttings on the natal dispersal distance could have been stronger in a wider spatial scale because the average natal dispersal distance was about 2 km and the spatial scale that we used was 500 m.The lack of forest management information outside the study area restricted the size of the largest spatial scale.
About 35 % of the breeding density decrease was not explained by the forest management effects.There may still be some other factors lying under the decline of the willow tit population.The global warming can be one of them, though the warming climate has been found to

Table B1
AIC-comparison tables representing the model, AIC and ΔAIC of distance to the nearest neighbour, natal dispersal distance and breeding dispersal distance analysis.Model structures are from the final global models (Table 2).In each comparison both the spatial (buffer) and temporal (year) attributes in clear-cutting and thinning variables changes.The spatial attributes in variables forest type gradient and non-forested areas included in the distance to the nearest neighbour and natal dispersal distance analysis were built to match the spatial scales of the clear-cutting and thinnings variables.The other variables remained constant.A = 0-30 years old clearcuttings and thinnings within (A1) 57 m, (A2) 113 m, (A3) 329 m or (A4) 500 m of the nesting site, B = 0-10 years old clear-cuttings and thinnings within (B1) 57 m, (B2) 113 m, (B3) 329 m or (B4) 500 m of the nesting site and C = 0-4 years old clear-cuttings and thinnings within (C1) 57 m, (C2) 113 m, (C3) 329 m or (C4) 500 m from the nesting site.The chosen model (lowest AIC-value) is bolded.

Table C3
The model averaged (conditional average, scaled effects) results of the linear mixed model analysis explaining willow tits nesting dispersal distance (square root-transformed) in 1990-2020.ID (ring number) and territory (correcting the male and female breeding site dependency) were assigned as random variables in the analysis.reduce the mismatch in willow tit breeding synchrony when considering the timing of breeding and the highest caterpillar peak (Vatka et al., 2011).In addition, there might be increased interspecific competition in willow tit habitats as the temperate-origin species, the great tit (Parus major) and the blue tit (Cyanistes caeruleus) have spread to the northern Finland and increased in numbers in the last few decades (Valkama et al., 2011), also in our study area.Though the great and blue tits are secondary cavity-nesting species they might increase the competition of available resources (Siriwardena, 2004).The interspecific competition hypothesis has been studied in Britain by Siriwardena (2004) and Lewis et al. (2007) and they did not find support for it.However, the willow tits in northern Finland have no long-term history with the great and blue tits as they have in the middle Europe which can pronounce a competitive relationship, though this has not yet been studied.

Conclusions
Overall, our findings point out that the even-aged forest management has caused habitat loss for the willow tit which has affected the spatial breeding distribution measures and probably has a major role in the population decline of the species.It seems that clear-cuttings and thinnings have at least 0-30 years lasting cumulative effect on the breeding willow tits.The effect is probably caused by a long-term habitat degradation and loss due to the forest management.We have shown here that the effects caused by clear-cutting the forests are more severe than of thinning due to the direct loss of older forest habitats.Thinning of forests can also have a remarkable effect on the breeding willow tits if all  the snags suitable for being a willow tit's nest are removed from the forest.
For conserving this endangered species, we suggest considering other forest management methods than clear-cutting e.g., continuous-cover forestry as the main forest management practise whenever it is possible.While the increasing number of decaying snags improves the habitat quality for the willow tit, we recommend leaving enough smalldiameter deciduous decaying or soon decaying trees to areas that need to be cut.In the future, these will offer possible breeding opportunities for the willow tits.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.The distance to the nearest neighbour (sqrt-transf.) in relation to A) the proportion of 0-30 years old clear-cuttings within 500 m buffer (%), B) the proportion of 0-30 years old thinnings within 500 m buffer (%) and C) the number of snags within 500 m radius with the 95 % confidence intervals from the linear mixed model analysis of the fittest model.The selected variable effects are plotted from the model presented in Table3, by keeping the other variable effects on their mean values.

Fig. 2 .
Fig. 2. Natal dispersal distance (m) in relation to the proportion of 0-30 years old clear-cuttings within 500 m buffer (%) with 95 % confidence intervals from the linear mixed model analysis of the fittest model presented in Table4.The other variable effects were kept on their mean values.

Fig. 3 .
Fig. 3.The breeding dispersal distance (sqrt-transf.) in relation to the proportion of 0-30 years old clear-cuttings (%) within 57 m buffer with 95 % confidence intervals of the linear mixed model results of the averaged model, which includes all variables presented in the Table5.The other variable effects were kept on their mean values.

Fig. D1 .
Fig. D1.The relationship between distance to the nearest neighbour and proportion of A) 0-30 years old clear-cuttings and B) thinnings within 500 m buffer around the nesting site with the 95 % confidence intervals from the linear mixed model analysis.The dashed vertical and horizontal lines correspond approximately to the average proportion increase in clear-cuttings (Fig. D.1a) and thinnings (Fig. D.1b).

Fig. D2 .
Fig. D2.Empirical relationship between the distance to the nearest neighbour (m) and breeding density (pairs/km 2 ) in the willow tits 1990-2020.The equation of the best fitting curve is: Density = 84849*NNdist^-1.739(p < 0.001, R 2 = 0.9579).The corresponding clear-cutting proportions from Fig. D.1 are in solid lines and thinnings in dashed lines.

Table 2
Table 2).Model structures of the final global models and transformations of the response variables used in the analyses.

Table A1
Description of the forest management data: the buffer sizes, total buffer areas in hectares, age categories and average area of clear-cuttings and thinnings on the buffer in hectares.

Table A2
The forest type categories (i.e.forest vegetation eutrophy) of forest habitat type 2017 data from the National Forest Inventory provided by the National Resources Institute Finland.This data was used to calculate the forest type gradient data.

Table A3
Main forest habitat type categories of main forest habitat type 2017 data from National Forest Inventory provided by the National Resources Institute Finland.This data was used to calculate the proportion of non-forested areas.

Table C1
The results of scaled variables in linear mixed model analysis with function 'lmer' of the fittest model (fixed effects, scaled effects), explaining the willow tits nearest neighbour distances in 1990-2020.The response variable was square root-transformed.Year and female and male ID (ring number) were assigned as random variables in the analysis.