Population dynamics and regulation of the Middle Spotted Woodpecker in strictly protected and managed forests in Western Poland

Analysis of ecological time series allows to assess the relative contribution of density-dependent and density-independent factors influencing the regulation of the population of a species. In this study, we investigated whether the annual fluctuation in the population of Middle Spotted Woodpeckers Dendrocoptes medius can be attributed to the direct or delayed density-dependent regulation or to the density-independent factors such as weather conditions during the preceding winter and previous breeding season and/or disturbance in the habitat associated with sanitary felling and salvage logging of trees. This study was conducted in the Czeszewo Forest Reserve, which protects the remnants of semi-natural flood-plain forest, and in the Łówkowiec Forest, which is a managed oak-dominated stand, during 19 seasons (2001–2019). Time series diagnostic revealed that the growth rate of population of Middle Spotted Woodpecker in protected semi-natural riverine forest and in managed oak forest was determined by first-order negative feedback, suggesting an important role of direct density-dependent mortality in the regulation of an equilibrium in population size. The effect of density-dependent process varied between study sites, being lower in riverine forest, where population size has increased distinctly in recent years. In both study areas, we detected strong positive effects of wind chill temperature during the preceding winter on population growth rates. We did not find evidence for the assumption that the increase in population growth rate is affected by weather conditions in the previous breeding season and by the disturbance in the habitat associated with sanitary felling and salvage logging of trees. Our results suggest that Middle Spotted Woodpeckers can benefit from global climate warming by an enhanced rate of winter survival or an increase in the accessibility of food, e.g., insects and other arthropods. The lack of time-lag significant correlations between population fluctuations suggests local habitat conditions, e.g., forest management, can modify inter-annual variability in time series of the Middle Spotted Woodpecker.


Introduction
Long-term studies reveal the temporal fluctuations and trends affecting the size of the population in question (Royama 1992; White 2019). The primary goal of an ecological time series analysis is to assess the relative contribution of endogenous (density-dependent) and exogenous (densityindependent) factors affecting the regulation of the population of a particular species (Berryman and Turchin 2001). In spite of its potential, time series analysis has been seldom applied in the study of population dynamics of woodpeckers. A few studies have examined the long-term ecological time series by focusing only on density-independent factors, e.g., weather conditions during the winter and food supply during the spring, affecting population dynamics of woodpeckers (Nilsson et al. 1992;Saari and Mikusiński 1996;Wesołowski and Tomiałojć 1997;Steen et al. 2006;Selås et al. 2008). There is no previous study that describes the diagnosis and evaluation of density-dependent feedback mechanism. However, population fluctuations might be caused due to the interactions between exogenous influences and endogenous feedback mechanisms (Royama 1992;Saether et al. 2016).
In this study, we analyze a long-term time series of two populations of the Middle Spotted Woodpecker (Dendrocoptes medius) from two types of forests: fully protected and managed forests. The nonmigratory Middle Spotted Woodpecker, distributed over large parts of Europe from North-Western Spain to western Russia and from the coasts of the Mediterranean to the North and Baltic Sea, in Anatolia and the Zagros Mountains, and associated with mature, rough-barked deciduous forests (Pasinelli 2003), is listed in the Annex I of the European Union Birds Directive. This implies that the Middle Spotted Woodpecker is the subject of special measures of conservation in relation to their habitat to ensure their survival and reproduction in the area of distribution in European Union. Thus, understanding the mechanism regulating the long-term population dynamics may have practical implications for species conservation and management. Moreover, data describing the effect of weather conditions on population size of the Middle Spotted Woodpecker are scarce. To our knowledge, only one study has explored the effect of winter temperature and food supply to population fluctuations of this species (Wesołowski and Tomiałojć 1997). Earlier studies have found that low temperatures and high amount of rainfall during the nestling phase negatively affect their breeding performance (Pasinelli 2001;Kosiński and Ksit 2006;Kossenko and Kaygorodova 2007); however, the role of weather conditions during the previous breeding season in population dynamics is still unknown.
In recent decades, the weather conditions are modified by climatic changes. A lot of studies have documented the ecological effects of changes in weather patterns due to warm temperatures on living organisms (Visser et al. 2001;Ahas and Aasa 2006;Schaper et al. 2012). However, data about the effect of climate change on the population of woodpeckers are scarce (Schiegg et al. 2002;Koenig et al. 2017). Information about how weather conditions affect annual fluctuations in the population of Middle Spotted Woodpeckers can help to understand the ecological effects of global climate warming on their population.
Annual fluctuation in the population of Middle Spotted Woodpecker may also be caused due to the disturbance in the habitat which is associated with sanitary felling and salvage logging of trees, especially Oak Trees (Quercus sp.), a critical tree species for this species of woodpeckers (Pasinelli 2003). A considerable increase in the number of oak trees infested with wood-boring insects in recent decades, driven primarily by drought stress, is a widespread phenomenon in Europe (Thomas et al. 2002;Sallé et al. 2014). Sanitary felling is the most effective way to reduce bark and wood-boring beetles in managed forests, and salvage logging allows to recoup economic value of wood that would otherwise be lost (Thorn et al. 2018). However, the effect of these forestry practices on the population dynamics of woodpeckers, e.g., the Middle Spotted Woodpecker, is poorly known (Wiesner and Klaus 2018).
Therefore, in this study, we investigated whether the annual fluctuation in the population of Middle Spotted Woodpecker can be attributed to direct or delayed densitydependent regulation and to density-independent factors, such as weather conditions and habitat disturbance. Based on the previous literature, we propose the following hypotheses: (a) per capita rates of population change should increase with increasing temperature and (b) decrease with increasing amount of rainfall in the previous breeding season, (c) per capita rates of change should increase with decreasing severity of the preceding winter, and (d) per capita rate of population change should be inversely correlated with the volume of deadwood removed in the previous year.

Study area
This study was conducted in the Czeszewo Forest Reserve (CFR; 17°31ʹ E, 52°09ʹ N) and the Łówkowiec Forest (LF; 17°34ʹ E, 51°43ʹ N), Greater Poland voivodeship, Western Poland. CFR (222.6 ha in size) protects the remnants of semi-natural flood-plain forest, which encompasses 186.5 ha of forest area with ca. 5% (8.5 ha) area covered by coniferous stands. The vegetation consists of Quercus-Fraxinus-Ulmus (Fraxino-Ulmetum) woodland in the flooded parts and Quercus-Carpinus (Stellario-carpinetum) forest on the higher ground. The rest of the study area (39 ha, 17%) is formed by oxbow lakes and meadows. CFR differs in terms of age, species composition, and spatial structures because of previous human impact. The oldest, near-natural stands are distinguished from younger, human-altered stands with regard to the height of the trees (latest measurements revealed that the highest dead ashes and alive oaks reach up to 42 m), the multi-story profile of stands (4-5 layers), and the high diversity of tree species (up to 7 species). The proportion of broadleaved stands (mainly oakdominated) that are > 80 years (suitable for Middle Spotted Woodpeckers; Kosiński 2006) Fig. 1A). The growing stock volume of oak wood was assessed on 235 m 3 /ha, and broadleaved wood on 403 m 3 /ha in 2009 (based on data from Forest Data Bank; https ://www.bdl. lasy.gov.pl/porta l/mapy-en). However, in 2019 the stock volume of oak wood decreased to 218 m 3 /ha and broadleaved wood to 331 m 3 /ha, respectively (based on data from Forest Management Plan 2019-2028; Jarocin Forest District, unpubl. data). The decreased of growing stock volume was caused by the mass mortality of oaks and Ashes Fraxinus excelsior affected by lethal attacks from Oak Tree Borers Agrilus biguttatus and fungal disease caused by the ascomycete Hymenoscyphus fraxineus, respectively. The most of dead trees fell over during the last 10 years. Moreover, on 11 August 2017, some tree stands of CFR were damaged or destroyed due to the effect of the mesoscale convective system by gusty winds (Chojnacka-Ożga and Ożga 2018). CFR is a part of both the Natura 2000 Special Protection Area Middle Warta River Valley PLB300002, with the Middle Spotted Woodpecker as one of the target species (Wilk et al. 2010), and the Ostoja Nadwarciańska Special Area of Conservation PLH300009. Further characteristics of the forest stands and their management can be found elsewhere (Kosiński et al. 2018).
LF, formerly called as Stare Budy (Kosiński and Winiecki 2005;Kosiński and Hybsz 2006), is an isolated forest area of 643.5 ha (630.8 ha of forest area). It is a part of the largest concentration of old oak forests (mainly Quercetea roboripetreae, locally Calamagrostis arundinaceae-Quercetum petraeae) in Poland (Krahl-Urban 1941, 1943Kasprowicz 2010). The vegetation of LF consists of oak-dominated stands (617.5 ha; 98% of forest area). This forest is characterized by two-story stands and has low diversity of tree species (up to 2-3 species per stand, but in some forest area, there might be exceptionally more than 2-3) (Szmyt 2017 Fig. 1A). The growing stock volume of oak wood was assessed on 283 m 3 /ha, and volume of broadleaved wood on 375 m 3 /ha in 2019 (based on data from Forest Data Bank; https ://www.bdl.lasy.gov.pl/porta l/mapy-en). During the study period, the cumulative area of forest harvested (habitat loss) reached ca. 65 ha (Supplementary Fig. 2A). The mean volume of broadleaved deadwood harvested was 1.8 m 3 per hectare per year; maximum 6.2 m 3 /ha in 2018 (Supplementary Fig. 2A). LF is a part of Natura 2000 Special Protection Area Krotoszyn Oak Forest PLB300007, with the Middle Spotted Woodpecker as the target species, and of Uroczyska Płyty Krotoszyńskiej Special Area of Conservation PLH300002. Krotoszyn Oak Forest is the second most important Middle Spotted Woodpecker site in Poland, with an estimated population of ca. 480 pairs (Wilk et al. 2010).

Data on the woodpecker population
We counted Middle Spotted Woodpeckers from 2000 to 2019 in CFR and from 2001 to 2019 in LF; however, to compare the data, we used data from 2001 to 2019. To estimate the number and distribution of territories of Middle Spotted Woodpecker, the study areas were surveyed in each of the prebreeding seasons (from the second half of March till the end of April), which coincides with the peak of intraspecific aggression (Pasinelli et al. 2001). In CFR, we performed four to five complete censuses during 2000-2005 and two to three censuses during the rest of the study period. In LF, we performed two to three complete censuses in each year except in the years 2014 and 2017, when only one census was conducted in April due to adverse weather conditions for the most part of spring. There was an interval of minimum 1 week between censuses. All counts were based on call-playback method, which has been recommended for Middle Spotted Woodpecker censuses (Kosiński and Winiecki 2003;Kosiński et al. 2004;Weggler et al. 2013). At selected points, we played rattle-calls and advertising-calls of males to elicit vocal responses from birds. To reduce the probability of some individuals being attracted away from their territories by the playback technique, we applied the following parameters: a minimum distance of 150-200 m between points, a maximum time of 40 s for the stimulation in one bout, and an interval of 1-2 min for listening. After the first bird's response, the playback was stopped. We recorded and mapped the woodpecker's initial position, behavior, type of call, direction of movements, and sex, if possible. We paid special attention to register simultaneously active birds to separate neighboring pairs/territories. To locate the precise position of the bird, we walked along the forest roads, division lines, and shorelines of oxbow lakes in CFR. The observations began in the morning at 7-8 AM and continued until late afternoon. Our experience suggests that time of the day does not affect the detection of Middle Spotted Woodpeckers in the prebreeding season. This data helped us to estimate the number of birds holding territory in the study areas. Our long-term experience suggests that two registrations in the prebreeding period were sufficient to accept a territory (see also Weggler et al. 2013;Lõhmus et al. 2016). Moreover, even a single visit coinciding with the peak of territorial activity allows an experienced observer to find almost 100% of the territories in the prebreeding period (Kosiński et al. 2004;Kosiński unpublished data). In CFR, all Middle Spotted Woodpecker territories were checked to find a cavity, except during 2010-2012. We focused our field surveys mainly on the detection of calling nestlings in the second part of May and at the beginning of June. During 2001-2016, we found breeding cavities in 60-77% territories. However, during 2017-2019, the efficiency of cavity-searching method decreased to 40-56%, mainly due to the increase in fallen logs, which rendered the fieldwork much more time-consuming than in earlier years.

Environmental variables
Based on previous studies, we used environmental predictors that were biologically and ecologically relevant in our analysis. In winter, the heat loss to the surrounding environment (herein defined as the period from 1 December to 28 February) has substantially increased with decreasing temperatures and increasing wind speed (Bakken et al. 1991;Wolf and Walsberg 1996), which may affect winter survival (Saari and Mikusiński 1996). To include the effect of temperature and wind speed on changes in abundance over time, we converted them to the wind chill equivalent temperature in winter (WCETW; in °C) index (Osczevski and Bluestein 2005). Although WCETW has been previously found as highly correlated with mean temperature during the winter months (r = 0.98), we considered it to be a more appropriate parameter, because animals seem to respond to their perceived temperature rather than air temperature alone (Wiersma and Piersma 1994;Collop et al. 2016). Moreover, to test the effects of winter conditions on population dynamics, we used total rainfall and/or snowmelt (PW; in mm) and number of days of snowfall (DSF; excluding days when both rain and snow were recorded) during the winter months. To test the possible effects of weather conditions during the breeding season in the preceding year (t-1) on growth rates in the following season, we used mean temperature in May (TM; in °C) and total rainfall in May (PM; in mm). All values were assessed based on daily data obtained from the nearest meteorological station located in Kalisz (18º08´ E, 51º78´ N), ca. 36 km east from LF and 56 km southeast from CFR, provided by the https ://en.tutie mpo. net ( Supplementary Fig. 3A). To test the effect of intensity of forest management on population dynamics, we used the ln-volume of broadleaved deadwood removed (DWR) during 2004-2018 from Krotoszyn Forest Inspectorate (Supplementary Fig. 2A). The volume of the oak deadwood in most of the years constituted more than 90% of the broadleaved deadwood removed, except in 2005 when it reached 78%. Oak dieback and mortality was mainly affected by oak tree borers (Jabłoński et al. 2016).

Statistical analysis
We used the population analysis system (PAS) to explore the temporal changes in the population of Middle Spotted Woodpeckers using the ln-transformed data on abundance through time (Berryman 1999). We utilized autocorrelation function (ACF) to evaluate periodicity and partial rate correlation function (PRCF) to detect feedback dimensions in the time series (Berryman and Turchin 2001). The variability of time series was evaluated by the mean and the variance of return time (MRT and VRT, respectively). When time series was nonstationary, i.e., the mean value of the series changed gradually over time (P < 0.05, exceptionally P < 0.1), we detrended the data using the methods of Barryman (1994). First, we fitted a regression line Y t through the time series. Data were then adjusted to a constant mean using the following transformation equation: where N t is the new time series, X t is the untransformed time series, and ∼ N is the mean of untransformed series. We estimated the R-function (Berryman 1999) to identify the relative contribution of first-and second-order negative feedback (C1 and C2) by fitting the Gompertz function of the form: where R t is the per capita rate of increase (R t = X t -X t-1 ), A is the maximum per-capita rate of change, X t is the lntransformed data on abundance in year t, and e t represents the process noise. Subsequently, we inspected the PRCF to identify whether the population dynamics was dominated by first-(density-dependent) or second-order (density-independent) negative feedback (Berryman and Turchin 2001). The significance band of both ACF and PRCF was assessed by Bartlett's criterion (± 1.96/n 1/2 , where n is the number of observations). To assess relationship between both time series, we used cross-correlation function (CCF) in the "tseries" package. The magnitude of variability in the population was assessed using the coefficient of variation (CV) of the stationary distribution of population sizes around the equilibrium density K (Saether et al. 2016).
When population abundance is regressed against weather conditions, there is a risk of finding spurious climate-abundance relationships, because all variables change across years. Although weather variables did not show a significant trend over time (P > 0.05; Supplementary Fig. 3A, Table 1A), the time series were detrended before conducting the analyses using the first difference method (Royama 1992;Steen et al. 2006;Selås et al. 2008). Detrending by year decouples interannual variation in environmental predictors from long-term directional change (Anderson et al. 2012;Iler et al. 2017). Because the variable DWR significantly increased over time (slope ± SE = 0.146 ± 0.054, t = 2.69, df = 13, P = 0.019; Supplementary Fig. 2A), it was first detrended using the methods of Barryman (1994) and then expressed as growth rate using the first difference. Before analyses, we investigated the structure of correlation of detrended weather variables using Pearson correlation coefficient and the variance inflation factor (VIF). Our variables were generally not strongly collinear (|r| ≤ 0.64, VIF ≤ 2.02; Dormann et al. 2013) (Supplementary Table 2A).
To test the effects of density-dependent feedback, with lag-order determined by PRCF function (ln-abundance in year t-1) and weather factors on the population dynamics (R t ), we used linear model (LM) described by Eq. (1). We constructed a set of 13 candidate models: the interceptonly model (model 1), model with density-dependent feedback (model 2), and 11 models with different combinations of variables describing weather conditions in winter (models 3-9) and in spring (models 10-12), including the full model with all explanatory variables (model 13). Because the suitable habitats for Middle Spotted Woodpeckers in LF decreased during the study period, by cutting of old forest stands (Supplementary Fig. 2A), we used the log of the forest area as an offset in each model to control for this process. To perform a global validation Regression line (solid) and 95% confidence limit (grey area) are given of LM assumptions, we used the "gvlma" package (Peña and Slate 2006). In all models, the global test of model assumption was satisfied (P > 0.05). Finally, in the case of LF, we added the variable DWR to the set of predictors in candidate models. The global tests of the model's assumptions were not satisfied (P < 0.05); therefore, we inspected data to find outliers. Influential points (n = 3; data for 2008, 2015, and 2019) were removed from the analyses. However, the time series was short after removing the outliers and due to the missing data of prior to 2004 (n = 10). Consequently, we analyzed only a maximum of four covariates per model (i.e., N t-1 , DWR, one of the weather predictors, and log of the forest area as an offset); full model was omitted in this analysis. In all models, the global test of LM assumptions was satisfied (P > 0.05).
The relative support for each model was assessed using Akaike's Information Criterion corrected for small sample size (AICc). The best parsimonious models (delta AICc < 2 compared to the best model) were selected, with their parameter estimates and standard errors (SEs) averaged (Burnham and Anderson 2002), using the "MuMIn" package (Bartoń 2015). For all models, Akaike weights were calculated to provide a quantitative measure of relative support for each competing model, with higher weights indicating better explanatory power. We used 95% confidence interval (CI) to assess the magnitude of the effect.

Changes in woodpecker abundance
The number of pairs of Middle Spotted Woodpecker in CFR varied between 29 in 2003 and 2009 (1.6 pairs 10 ha −1 of the forest area) and 47 in 2018 (2.5 pairs 10 ha −1 ) and in LF between 48 in 2003 (0.8 pairs 10 ha −1 of the forest area) and 70 in 2008 and 2015 (1.2 pairs 10 ha −1 ). A linear regression analysis revealed that temporal trends in the abundance of Middle Spotted Woodpeckers in both study areas increased over time; however, only in the case of CFR, this trend was statistically significant (CFR: slope ± SE = 0.605 ± 0.161, t = 3.77, df = 17, P = 0.002; LF: slope ± SE = 0.410 ± 0.230, t = 1.78, df = 17, P = 0.093) (Fig. 1). In CFR, population time series was monotonic (MRT = 2.20, VRT = 6.76) and nonstationary. It was affected by an increase in abundance since 2016 (Fig. 1). Population dynamics of Middle Spotted Woodpeckers in LF was characterized by periodic fluctuations around zero with a short mean return time (MRT = 1.51, VRT = 0.64). The magnitude of population variability (CV) of the stationary distribution of ln-abundance around the equilibrium density K was 0.105 in the case of CFR and 0.087 in the case of LF. PRCF revealed that the population dynamics in CFR and LF (after detrending the time series) were dominated by the density-dependent firstorder negative feedback (ln N t-1 ) (Fig. 2). Using the multiple linear regression analysis, we detected that in both the study areas, first-order feedback process explained 29% (CFR) and 43% (LF) of the variation in the data ( Table 1). The analysis of relationships between both time series revealed lack of time-lag significant correlation between  Fig. 4A). Therefore, we modeled the effect of covariates on woodpecker dynamics separately for each studied site.

Modeling population abundance
In the case of CFR, two best models (delta AICc ≤ 2 compared to the best model) accounted for just over half of Akaike model weight (0.53), indicating that they were better supported by the data than that of the remaining models (Table 2). These models included three variables: ln-abundance in year t-1 (N t-1 ), wind chill temperature in winter (WCETW), and number of days with snowfall (DSF). The parameters based on model-averaging suggest that the per capita rate of increase (R t ) was positively affected by WCETW as this variable had effect sizes larger than their SE, and the 95% CI did not include zero (Table 3; Fig. 3). In the case of LF, the three best models included those included for CFR (sum of Akaike weight of models-0.74) ( Table 2). Averaged model coefficients suggest that the R t was positively affected by WCETW and negatively affected by N t-1 as both variables had effect sizes larger than their SEs, and their 95% CIs did not include zero (Table 3; Fig. 3). When the ln-volume of deadwood removed (DWR) was included in the analysis, then the two best parsimonious models contained first-order negative feedback, DWR, and WCETW, which accounted for almost all of the Akaike weight (0.88) ( Table 4). Effect sizes of variables, their SEs, and 95% CIs of estimates suggest that R t decreased with increasing N t-1 and increased with increasing WCETW (Table 5).

Discussion
We found that population dynamics of Middle Spotted Woodpeckers is affected both by density-dependent firstorder feedback process and stochastic density-independent mechanism, i.e., WCETW. The first-order feedback process is typical for populations constrained by intraspecific competition for food, mates, or territories (Berryman 1999), which is the case in Middle Spotted Woodpecker (Pasinelli 2003). The strength of this process varied between study sites, being lower in CFR, where population size has increased distinctly in recent years as compared to LF. This result is in-line with earlier suggestions that strength of density-dependence in birds typically declined as population size increased (Saether et al. 2016).
A recent meta-analysis that partitioned the effects of density-dependence and environmental stochasticity on the population growth rate of birds suggests that population dynamics are most strongly affected by environmentdrive variation in the number of new recruits entering the population (Saether et al. 2016). According to our results, WCETW positively affected per-capita rate of increase in the population of Middle Spotted Woodpecker in both study sites. This result is in-line with previous studies that examined the role of weather conditions during the winter (most frequently temperature) for the survival of different species of woodpecker in Europe (Nilsson et al. 1992;Saari and Mikusiński 1996;Steen et al. 2006;Selås et al. 2008). For example, in Białowieża National Park (Poland), the number of Middle Spotted Woodpeckers during breeding season was positively correlated with mean temperatures in the preceding winter months (December-March) (Wesołowski and Tomiałojć 1997). Earlier studies have revealed that energy expended on thermoregulation during winter is substantially increased with decreasing temperatures and increasing wind speed (Calder and King 1974;Bakken et al. 1991;Wolf and Walsberg 1996), whereas energy intake is often reduced due to lowered availability of food per se, by snow cover or ice glazing on vegetation (Pettersson 1983;Peterson et al. 1989;Saari and Mikusiński 1996;Newton 1998;Wesołowski et al. 2018). In consequence, birds spend at least part of the winter season at temperatures well below their thermoneutral zone (Eduardo et al. 2010). Both low wind chill temperature and low availability of live trunk-surface arthropod during winter probably increased the rate of mortality of Middle Spotted Woodpecker compared to warmer winter. Contrary to our expectations, we did not find the negative effect of number of days of snowfall and the total amount of rainfall and/or snowmelt during winter on population growth rate. Some authors suggest that snow cover and ice glazing on vegetation could be particularly destructive for birds that Table 1 Coefficients of variables explaining the effects of densitydependent feedback, i.e., ln-abundance in year n-1 (N t-1 ) and year n-2 (N t-2 ) on the population growth rate (R t ) of the Middle Spotted Woodpecker in Czeszewo Forest Reserve (CFR) and Łówkowiec Forest (LF) Significance code: a 0.001, b 0.01, c 0.05, d 0.1 are naive with respect to such conditions (Wesołowski et al. 2018), which was not likely the case in CFR and LF. However, our data showed that extreme amount of precipitation (mostly rainfall) in some winters, as it was observed in 2016/2017, might cause a sharp decline in the population abundance in the upcoming breeding season ( Supplementary Fig. 3A). We did not find evidence that an increase in the population growth rate is affected by weather conditions in the previous breeding season; we expected that higher reproductive rate due to an increase in the temperature and a decrease in the amount of rainfall in May could raise recruitment to breeding population, i.e., increase in growth rate. This is surprising, because according to previous studies, there were negative effects of low temperatures and high amount of rainfall during the nestling phase on breeding performance (Pasinelli 2001;Kosiński and Ksit 2006;Kossenko and Kaygorodova 2007). Pasinelli (2001) found a positive association between the average temperature during the nestling phase and nesting success and found a negative correlation between the mean amount of rainfall during the nestling period and the number of young fledged in northeastern Switzerland. Kossenko and Kaygorodova (2007) revealed negative correlation between fledgling production and amount of rainfall during nestling period in southwestern Russia; however, they obtained these results in food-poor forest fragments only. Moreover, Kosiński and Ksit (2006) found negative effect of rainfall during nestling period on the Table 2 Model-selection statistic for the analysis of relationship between ln-abundance in year n-1 (N t-1 ), weather conditions and per capita rate of increase (R t ) of the Middle Spotted Woodpecker in Czeszewo Forest Reserve (CFR) and Łówkowiec Forest (LF) The best models (separated by ΔAICc ≤ 2) are highlighted in bold WCETW-wind chill equivalent temperature in winter, PM-total rainfall and/or snowmelt, DSF-days of snow fall, TM-Temperature in May, PM-total rainfall in May a Degrees of freedom b Akaike's information criterion c Difference in AICc relative to the model with the lowest AICc d Models are ranked according to their Akaike weight; higher weights indicate more parsimonious models probability of nestling survival in CFR; however, the effect of year was not controlled in this study. On the one hand, the lack of association between population growth rate and amount of rainfall and temperature in the previous breeding season may be attributed to different methodological approaches. Earlier studies focused only on fledgling period, whereas we studied the effect of precipitation in broader, less accurately determined, time frame, including incubation (Kosiński and Ksit 2006), because accurate determination of nestling period and its calibration with weather factors in different years was not possible (Steen et al. 2006). It is also possible that variation in the amount of rainfall and temperature in May across years was insufficient to allow the detection of population growth rate-weather relationships. For example, mean amount of rainfall in our study area during May did not exceed 3.1 mm/day, exceptionally 4.8 mm/day in 2011, whereas in Switzerland the amount of rainfall during fledgling period, which negatively affected breeding success, exceeded 5.0 mm/day (Pasinelli 2001).
On the other hand, Saether et al. (2016) found that fecundity explains only a small proportion of variation in population variability in birds, whereas population growth is mainly affected by density-dependent mortality of both juveniles and adults. In general, our results are in-line with the hypothesis proposed by Lack (1966) that changes in population size in birds are determined by density-dependent mortality during the nonbreeding season and density-independent recruitment to the population (Saether et al. 2016). We did not confirm the hypothesis that rate of population change in LF is strongly and negatively affected by the volume of deadwood removed in previous year. In fact, the mean volume of deadwood harvested annually represents a small proportion of the growing stock volume of broadleaved wood (0.5%). Moreover, although the extent and intensity of sanitary felling and salvage logging increased distinctly in 2018 (Supplementary Fig. 2A), the breeding population of the Middle Spotted Woodpecker in 2019 has increased compared to previous years. Furthermore, recent study in forests around Jena (Germany) found that harvesting intensity did not affect temporal trends of the Middle Spotted Woodpecker; trends in plots with low and high timber harvest intensities changed in the same way in 2001(Wiesner and Klaus 2018. Taken together, this data suggests that sanitary felling and salvage logging might have limited negative impact on the temporal trend and growth rate of the Middle Spotted Woodpecker. However, long-term studies are needed to confirm these findings. We found that the abundance of Middle Spotted Woodpeckers in both study sites increased over time; however, the temporal trend in LF was less apparent. It was probably caused by reduce of forest area accessible for Middle Spotted Woodpeckers over study period. These results are in line with recent studies documenting increase in population size and range shift of the Middle Spotted Woodpecker in different parts of Europe (e.g., Berndt 2015;Lõhmus et al. 2016;Chylarecki et al. 2018). The most spectacular increase in number of territories was registered in Switzerland, where population size is now larger by factor 2 to 2.5 than in 2000 (Schuck et al. 2018). An increase in the population abundance is most frequently attributed to the Table 3 Model coefficients of the variables included in the best ranked models explaining the effects of ln-abundance in year n-1 (N t-1 ), wind chill equivalent temperature in winter (WCETW) and days of snow fall (DSF) on per capita rate of increase (R t ) of the Middle Spotted Woodpecker in Czeszewo Forest Reserve (CFR) and Łówkowiec Forest (LF) Estimation based on model-averaging. See Table 2  Although trends in population time series were generally similar in CFR and LF, we did not find the time-lag significant correlation between them. This suggests that local conditions of habitat, e.g., forest management, modified inter-annual variability in the time series of the Middle Spotted Woodpecker in LF. However, some synchronization in annual variation between plots with low and high timber harvest intensities was observed in Germany (Wiesner and Klaus 2018). The effect of global climate warming on fecundity in Middle Spotted Woodpeckers is not known. The study on red-cockaded woodpecker Picoides borealis revealed that birds were benefitted from the global warming as they started breeding earlier in warmer years. However, overall reproductive output decreased slightly over the study period (Schiegg et al. 2002). It has been suggested that forest thinning could increase population abundance of Middle Spotted Woodpeckers through providing a high density of oak trees with sunny crowns, which are highly preferred as foraging sites (Pasinelli and Hegelbach 1997;Schuck et al. 2018). As stated earlier, sanitary felling and salvage logging did not affect population growth rate in LF. Furthermore, we did not observe a temporal trend in the Middle Spotted Woodpecker abundance in CFR during 2000-2015, i.e., before and after the mass mortality of oaks and ashes (Kosiński et al. 2018). On the one hand, these results suggest that forest thinning and natural disturbances probably do not distinctly improve the accessibility of food for Middle Spotted Woodpeckers, mainly smaller bark surface-dwelling arthropods  Table 5 Model coefficients of the variables included in the best ranked models explaining the effects of ln-abundance in year n-1 (N t-1 ), wind chill equivalent temperature in winter (WCETW) and lnamount of dead wood removed (DWR) on per capita rate of increase (R t ) of the Middle Spotted Woodpecker in Łówkowiec Forest (LF) Estimation based on model-averaging. See Table 4 (Pasinelli 2003), or their effects are delayed. On the other hand, an increase in the accessibility of food can be more strongly affected by global climate warming. It was found that insects respond positively to increase in the temperature which favors reproductive rate and decreases winter mortality. Moreover, indirect positive effects on insects can occur through temperature and drought stress of host trees (Pureswaran 2018). The increase in population size of the Middle Spotted Woodpecker is also attributed to the availability of dead trees, which may improve their food base and play an important role in nest-site selection (Pasinelli 2000;Weggler et al. 2013;Schuck et al. 2018). However, as previously stated, we did not find any evidence that the Middle Spotted Woodpecker tended to increase in abundance after a pulse in deadwood availability in CFR (Kosiński et al. 2018). Moreover, we did not observe a negative effect of deadwood removal on population abundance in LF. We suggest that in old forest stands, as in CFR and LF, with a large number of oaks and other rough-barked tree species (in CFR) with characteristics suited for cavity excavation (dead limbs and branches, limbholes, open wounds/scars, fruiting body of tree fungus), the presence of dead trees does not limit the population size of Middle Spotted Woodpeckers (Pasinelli 2000;Lovaty 2002;Kosiński et al. 2018). However, recent studies showed that dead trees played a vital role as a resource for cavity excavation for this species (Kosiński et al. 2018). Excavation in dead substrate probably reduces the energy expenditure of Middle Spotted Woodpeckers, which are characterized by weak excavating abilities, and maybe adaptive by offering time and energy savings.

Conclusion
According to the results of this study, population dynamics of Middle Spotted Woodpeckers are driven both by density-dependent first-order feedback mechanism and stochastic density-independent factor, i.e., winter severity. This is probably the first report demonstrating the role of both mechanisms in the regulation of woodpecker populations. We did not find evidence that the population growth rate is affected by weather conditions in the previous breeding season and by the anthropogenic disturbance in the habitat. Our results suggest that increase in population abundance in recent years can be affected by global climate change. Middle spotted woodpeckers can benefit from warm temperatures, which might help to survive during the winters or by an increase in the accessibility of food such as insects and other arthropods. The lack of time-lag significant correlations between population fluctuations suggests local conditions of their habitat, e.g., forest management, which can modify inter-annual variability in time series of the Middle Spotted Woodpecker.