This paper reports the long-term numerical trends of the thirty common forest bird species and explores changes in the community composition in the three main types of old-growth stands in the Białowieża National Park (E Poland, hereafter BNP) over 45 years (1975–2019). We present recent (2015–2019) data on abundance of birds for the seven study plots and pool them with the time series collected since 1975. The numbers of individual bird species strongly fluctuated, with most of the species showing alternating phases (the initial periods of population growth followed by the periods of population decline or stability). The numbers of 19 species increased; maximum growths by c. 3–5% per year included Columba palumbus, Dendrocopos major, Sylvia atricapilla and Regulus ignicapilla. Among a few declining species, Ficedula hypoleuca, Phylloscopus sibilatrix and F. parva experienced the strongest declines, respectively by 4.0%, 2.7% and 2.2% per year. Mostly the same species bred in the plots in the 1970s and in recent years, indicating a stable species pool. The total abundance peaked around 2005, declining thereafter in deciduous stands, but increasing further (along with the species richness) in the coniferous stands. The similarity index between the study plots (beta-diversity) changed little over 45 years; ash-alder and lime-hornbeam stands remained most similar, while coniferous forests stood more apart. The changes found in the old-growth stands of BNP (mostly coniferous fragments) could be partly explained by natural modification of the habitat structure and the processes acting in the large geographical scales, within or outside of the breeding grounds. The long-term studies such as the one in the BNP reported here, provide a basis for the rates of natural turnover in the bird communities in pristine habitats, directly unaffected by human impact.
Long-term research in ecology and evolution is fundamental for the understanding of the complex dynamics of natural systems. The results of long-term studies cannot be replaced by outcomes of short-term ones, since the former are necessary to enable detection of unpredictable, episodic, cumulative, slow-acting, non-linear phenomena, or processes subjected to thresholds of environmental change (Magurran et al. 2010, Silvertown et al. 2010). Long-term studies contribute a disproportionately large amount of new knowledge, helpful in understanding relationships and processes in wild populations and systems, which are helpful in advancing ecological and evolutionary theory. Despite their contributions becoming ever more relevant in the era of accelerating environmental changes, the current model of science funding makes running the long-term studies exceedingly difficult (Hayes & Schradin 2017, Hughes et al. 2017, Kuebbing et al. 2018 and references therein). Birds and bird communities are viewed as one of the important indicator groups (Fraxeidas et al. 2020), well suited to study environmental changes since decades (Furness & Greenwood 1993), and the same applies to their ecology (Fiedler 2009). However, observations of the bird community dynamics in forests least affected by direct human impact are rare, and long-term series of data collection in such conditions are even less common. Apart from the Białowieża National Park (BNP thereafter, this study) only a few other long-term study-programs focused on forest avian communities, i.e. research in Birdsong Valley in S Sweden (Svensson et al. 2010, Å. Lindström, pers. comm.), and in the Hubbard Brook Experimental Forest, in New Hampshire, USA (Holmes 2011, Holmes & Likens 2016).
Numerous anthropogenic changes of environmental conditions, both indirect (e.g. via climate change) and direct (deforestation, fragmentation, introduction of alien species, removal of old trees and dead wood) can markedly influence forest bird communities (reviews in Tucker & Evans 1997, Kampichler et al. 2014). However, measuring these anthropogenic effects remains difficult without knowing the rate of natural turnover (or ‘background dynamics’) in the communities free from such human impact. Such a reference information allows to differentiate between the effects of direct human actions (e.g. forestry operations) and variation due to other (‘natural’) reasons. The dynamics of avian communities in areas least affected by human influence are extraordinarily helpful in this respect (Tomiałojć et al. 1984, Wiens 1989, Wesołowski & Tomiałojć 1997, Gatter 2000, Wesołowski & Cholewa 2009, Magurran et al. 2010, Kampichler et al. 2014). This reasoning led to initiation of the breeding bird community studies in the strictly protected part of the BNP in 1975 (Tomiałojć et al. 1984).
The Białowieża Forest is one of a few places in temperate Europe where fragments of a primeval lowland forest have survived (Tomiałojć & Wesołowski 2005, Wesołowski 2005, 2007a). The pristine features, which include the enormous stature of trees, large amount of dead wood, an almost intact community of predators and large herbivores, and a great diversity of superabundant tree holes, have been well-documented (Faliński 1986, Tomiałojć 1991, Jędrzejewska & Jędrzejewski 1998, Wesołowski 2007b, Kuijper et al. 2013, Jaroszewicz et al. 2019, Wesołowski et al. 2021). Thus, the BNP data may provide a benchmark for bird community studies in modified woodlands (Tomiałojć et al. 1984, Tomiałojć & Wesołowski 1990, 2004, Wesołowski & Tomiałojć 1995, 1997, Wesołowski et al. 2006, 2010, 2015, Wesołowski 2007b, Kampichler et al. 2014, Czeszczewik et al. 2015).
The breeding bird community occurring in BNP exhibits a remarkable level of the compositional stability (Tomiałojć & Wesołowski 1994, Wesołowski et al. 2002, 2010, 2015, Kampichler et al. 2014). The local species richness (alpha-diversity) is very high, possibly one of the highest in any temperate forest (Tomiałojć & Wesołowski 2004). The majority of bird species can be found in a wide array of old-growth stands due to the fine-grain mosaic of the stands and the flexibility of breeding birds in their habitat requirements (Tomiałojć et al. 1984, Wesołowski et al. 2010, 2015). Consequently, the beta-diversity, defined as a turn-over of species among habitats, is relatively low. Even the formation of several hectares large windfalls did not result in marked exchange of species (Fuller 2000). However, some large-scale natural transformations of forest structure do affect the composition of the breeding bird community there. Deciduous trees have overgrown gaps in the coniferous forest (Walankiewicz 2002) that was later colonised by birds from ‘deciduous’ habitats, leading to an increase of bird species richness in mixed-coniferous stands (Tomiałojć & Wesołowski 2004, Wesołowski et al. 2015). Concurrent disappearance of sharp forest edges resulted in the disappearance or the decline of the open-landscape species (Wesołowski et al. 2010, 2015).
In this paper, we characterised the breeding bird community of BNP in 2015–2019 to complete time series on bird abundances, collected since 1975. All the analyses were done on 45 years long time series (1975–2019). To provide the longest trend assessment of forest birds in this part of Europe, we then determined the long-term changes in abundance of the 30 most numerous species. To answer our main research question — how quickly communities can change in pristine forests without direct human intervention (the natural turnover) — we assessed the degree of bird community modification by testing for any directional changes in species richness (alpha diversity) and similarity of bird communities among forest stands (beta-diversity). We expected that the most pronounced changes would occur in the conifer stands, which became increasingly similar to the lime-hornbeam ones due to overgrowing of windfall gaps over the study period, presumably leading to an increasing similarity of the bird breeding community between the coniferous and deciduous stands. We expected that the post-‘disturbance’ changes of bird community would occur faster in coniferous stands than in other deciduous fragments, reflecting the dynamics of natural modification of the pristine habitats. We discussed whether the observed numerical variation in the bird community may relate to the cumulative effects of the long-term habitat transformations, such as the disappearance of Norway Spruces Picea abies and European Ashes Fraxinus excelsior (in riverine forests), and the colonisation of canopy gaps by deciduous trees in the coniferous stands.
The study was carried out in 2015–2019 in the best-preserved, central part of the Białowieża Forest (52°29′–52°57′N, 23°31′–24°21′E), within the Białowieża National Park, strictly protected since 1921. Detailed descriptions of this area are given elsewhere (e.g. Wesołowski & Tomiałojć 1995, Tomiałojć & Wesołowski 2004, Okołów et al. 2009, Wesołowski et al. 2018). The census work was carried out within the seven permanent census plots established in 1975–1976 (total area of 187.5 ha, Fig. 1; Tomiałojć et al. 1984). All plots were covered by 50 × 50 m grid system of orientation marks. The detailed descriptions and photographs of the plots in the beginning of the study are available in Tomiałojć et al. (1984), Tomiałojć & Wesołowski (1990, 1994), Wesołowski & Tomiałojć (1995) and here we illustrate their current state (Figs. 2–4). Below, we describe the recent changes in the habitat structure.
Ash-alder riverine forest Fraxino-Alnetum — plot K (33 ha) situated at the forest edge (Fig. 1). It is a forest ‘peninsula’ adjoining, along c. 1000 m long border, former meadows in the Narewka river valley. Initially the boundary was sharp, as the meadows were regularly mowed. However, after cessation of mowing, young alders Alnus glutinosa colonised a 50–100 m wide belt of the meadows adjoining the old growth forest. Currently they form dense stands with trees c. 40 years old.
The plot is covered mostly by riverine stands on swampy organic soil, but there are also intrusions of c. 110–130 years old alder-birch Betula spp. regeneration and c. 6 ha of dry ‘islands‘of the lime-hornbeam type (see Tomiałojć et al. 1984). The riverine fragments were formerly covered mostly with alder (c. 40%), ash (c. 20%) and spruce. However, due to the continued ash dieback since c. 2000, caused by a novel fungal disease (Kowalski & Holdenrieder 2008), the share of ash strongly diminished in the canopy (Fig. 2C) and the density of shrubs (mainly currants Ribes sp.) increased. A number of canopy spruces were killed by bark beetles Ips typographus during the consecutive outbreaks (see coniferous plots below), and so spruces have become less numerous as well. Consequently, the already open-canopy forest became even more open, with the canopy cover reduced by up to 10–20% in some places and an increased number of uprooted trees. Uprooted trees were c. three times more frequent in this habitat than in the lime-hornbeam forest in the 1970s (Wesołowski 1983), while they became c. six time more frequent recently (Karpińska et al. 2022). The canopy gaps are overgrown by early-successional stages including rich, locally hardly passable undergrowth, composed of hazel Corylus avellana, young alders, Bird Cherry Padus avium, Rowan Sorbus aucuparia and Water Elder Viburnum opulus (Fig. 2D). The ground layer is luxuriant (locally up to 1.7 m tall in June) composed of nettles Urtica dioica, currants Ribes spp., reed Phragmites australis, ferns, sedges Carex spp., etc. (Fig. 2B). In comparison to the 1970s, the area has become much less swampy due to dramatically lowered ground water levels, which is probably due to more frequent snowless winters in recent years and an increased water outflow from the Narewka river.
Alder-swamp forest Ribeso nigri-Alnetum — plot L (25 ha). Situated c. 5 km from the forest edge (Fig. 1). The plot embraces the lower section of a small stream valley, flanked by lime-hornbeam stands on the higher ground. Along the valley sides, alder carrs gradually merge with ash-alder associations, both with admixture of spruce. It is a largely inundated area with a patchwork of almost permanent pools and tiny tree-covered dry islets formed by tree roots (Fig. 2A). These islets are overgrown by young canopy trees, as well as bushes (Alder Buckthorn Frangula alnus, willows Salix spp., Bird Cherry, hazel). Luxuriant ground layer (reaching over 2 m high in June) is composed of tall ferns, sedges, reeds and numerous other herb species (Fig. 2B). The canopy was open already in the 1970s, when the amount of uprooted and dead trees was high. During the study period this part of the forest changed relatively little, except that old ashes disappeared almost completely. Due to falling of canopy ashes and spruces the number of canopy gaps and amount of fallen logs on the ground increased in recent years. The area became also somewhat drier than in the 1970s.
Lime-hornbeam forest Tilio-Carpinetum — plots CM (24 ha), MS (30 ha), W (25.5 ha). The three plots in lime-hornbeam forest are located at progressively increasing distances from the forest margin. Plot W is situated at the former forest-edge, plot CM is c. 2 km from the edge, and plot MS is located 3 km inside the BNP (Fig. 1). We named it “oak-lime-hornbeam forest” in our earlier publications, but the “lime-hornbeam” reflects better both the actual tree composition of this forest habitat and the scientific name of this plant association. The stands developed on rich mineral soils that dominate in the BNP (c. 47% of its area). It is the most stratified habitat with the tallest and the thickest trees. The emergent spruces reach 45–50 m, and several deciduous tree species are over 40 m (Niechoda & Korbel 2011), while the circumference at breast height of Pedunculate Oaks Quercus robur can reach >600 cm and that of Small-leaved Limes Tilia cordata >500 cm (Niechoda & Korbel 2011, Fig. 3A, 3B). The upper canopy is composed of emergent spruces, oaks, limes and ashes, and the same tree species along with Norway Maples Acer platanoides, elms Ulmus spp. and locally birches and aspens Populus tremula form the main canopy. The under canopy mostly consists of hornbeams Carpinus betulus and younger individuals of other canopy trees. Except in gaps created by fallen trees, little light penetrates to the forest floor so that both shrub and herb layers are poorly developed and grow mostly in the gaps (Fig. 3C). The number of standing and fallen dead trees is lower than in riverine and conifer plots.
The plots are structurally similar, although the share of the main tree species varies locally. Hornbeams are the most abundant tree species in plot MS, and limes in plot CM (Wesołowski 1996, Maziarz & Broughton 2015). In 1980s', spruce formed c. 23% of tree stands in plots C and W and 12.5% in plot MS (Wesołowski 1996), but it dropped to 13–16% and 4–6% respectively in 2008 and 2015 (Maziarz & Broughton 2015; D. Czeszczewik and W. Walankiewicz unpubl. data). Elms remain most abundant (c. 7%) in plot W, and only single elm trees could be found in plots CM and MS in 2008 (Maziarz & Broughton 2015). The pattern of change in the share of some tree species was common in the whole strict reserve of the BNP (Zajączkowski et al. 2017).
Young spruces and spruce saplings, widespread in the 1970s and 1980s (Wesołowski 2011, Fig. 3D) virtually disappeared, and almost no spruce regeneration occurred in recent years. On the other hand, saplings of limes and hornbeams were present all the time and saplings of Norway Maples and Pedunculate Oaks became plentiful recently. Shrub layer remained poorly developed until the last decade, since when the saplings and young trees in gaps have begun to grow much faster than previously. In several old gaps, thickets of young trees 1–6 m tall have appeared, probably as a consequence of reduced grazing pressure by large herbivores; the numbers of Red Deer Cervus elaphus in the Białowieża Forest in the 2000s were reduced to a third of their 1990s level (Jędrzejewski et al. 2006).
The southern margin of plot W formed a sharp border with abutting fields (cf. photo in Tomiałojć et al. 1984) in the 1970s. However, after abandonment of farming, trees colonized the fields margins and a several dozen meters wide belt of secondary woods, composed of numerous species of deciduous trees and spruces, developed along the former edge. Currently, plot W is flanked mostly by c. 40 years old thickets and not by open areas as formerly.
Pine-bilberry coniferous forest Peucedano-Pinetum, Querco roboris-Pinetum — plots NW (25 ha) and NE (25 ha) situated c. 4 km from the forest edge. Mixed coniferous-deciduous stands are trophically one of the poorest local habitat types. Trees, mostly spruces and Scots Pines Pinus sylvestris with an admixture of birches, aspens and oaks, are of moderate to large size and grow very densely. Ground layer is low, composed mostly of mosses and berries Vaccinium spp. (Fig. 4A).
During the last 45 years the number of old pines and spruces decreased, and numerous canopy gaps appeared after windfalls, with the largest created by hurricanes in 1988 and 2002. In addition, bark beetle killed many spruces during outbreaks in 1983–1988, 1994–1997, 2001–2003, and 2012–2019 (Gutowski et al. 2021). Consequently, the number of living spruces that formed the canopy more than halved during the study (Fig. 4B). The amount of dead timber (mainly broken stumps or logs) was high already in 1998 (c. 100 snags/ha, Walankiewicz 2002), and it further increased in the last years. Ground layer remained poorly developed in the closed-canopy fragments, but young hornbeams, birches, limes, oaks, and some spruces rapidly colonised the windfall gaps (Fig. 4C, D) and lush herb vegetation (bracken Pteridum aquilinum, reed grass Calamagrostis spp.) developed. As a result, the formerly old growth, closed-canopy, conifer-dominated stands are now admixed with patches of deciduous trees typical for the lime-hornbeam stands up to 30 years old.
We applied an improved version of the mapping technique for censusing all breeding birds on the study plots to estimate densities of common species (Tomiałojć 1980, Verner 1985), used in all previous studies (Wesołowski et al. 2015 and references therein). The numbers of most species were registered with small errors, but in the case of some other, e.g. Turdus philomelos, Coccothraustes coccothraustes, Ficedula albicollis and probably Regulus spp., the breeding numbers could be underestimated by 20–33%, particularly when densities were high (Tomiałojć 1980, 2004, Tomiałojć & Lontkowski 1989, Walankiewicz et al. 1997). Since there are only seven plots studied annually since 1975, the results are representative for a group of common forest birds, mostly passerines, but inadequate for species holding large territories (like raptors). Within the commonest species, nearly all regularly breeding in BNP were represented in the plots.
Every year, at least ten visits (up to 11–12 to compensate for the adverse weather conditions or early onset of spring) were made between 5–10 April and 25 June to each plot. The observations began before sunrise and proceeded along marked transects 100 m apart on the plots (leaving the transect for minor detours when necessary). Each time a different route across the plot was chosen. Only plot K, the richest in bird species, was subdivided into two parts of c. 16 ha each, and either censused simultaneously by two observers, or by one observer on two consecutive mornings. Such a division was necessary to maintain a high quality of the censuses on large swampy area that was difficult to walk through. In order to achieve a high inter-plot comparability, all plots were visited on a rotational basis by all 6–7 observers.
Usually, one of the ten visits at the beginning of the breeding season was performed in the evening for mapping bird species that are active at dusk (Turdus spp., Erithacus rubecula, Scolopax rusticola and owls). To record all singing individuals during the short-lasting activity period of birds in the evening, each plot was censused simultaneously by two to four observers. All fieldwork was conducted by the experienced fieldworkers, and the new participants were admitted only after a period of apprenticeship (i.e. after making sure that they were able to gather field data of adequate quality).
All records from field maps were assembled on the individual species maps (scale 1:1000) for further evaluation. To assure the maximum consistency of the evaluation rules (see Morozov 1995), all estimates of cluster/territory numbers on the species maps were checked and later compared by a minimum of three experienced researchers to draw the final values. While drawing the ‘territories’ around the clusters of records, we relied mainly on the presence/absence of contemporary records. This helped to avoid an apparent tendency in mobile individuals/species to form double clusters instead of a single large territory (Tomiałojć 1980). As a rule, the records on the three different visits to the plots were required as a minimum to draw a cluster; a few exceptions included late arriving or inconspicuous species (e.g. Muscicapa striata, Locustella fluviatilis, Oriolus oriolus). In such cases, only two records indicating the presence of a territory were sufficient. Numerous located nests (especially for Sturnus vulgaris), and behavioural cues indicating the presence of active nests, such as carrying nest material or food, alarm calls, helped to establish the number of territories. The territories of polygynous (Wesołowski 1987) or bachelor males were treated as equivalent to those of monogamous pairs. In case of double-brooded species (e.g., Turdus philomelos), we assumed that repeated or second broods fell within the territories established in the beginning of spring.
The estimates of species richness may be somewhat higher due to inclusive ‘territories’ of large or rare bird species, partially overlapping with the study plots (relates to all species marked with ‘+’ in the tables).
To determine the long-term changes in the number of the individual bird species, and assess the variation in the composition of the breeding bird community, we pooled the data from 1975–2014 with the data from 2015–2019, all collected in the same way on the seven permanent study plots (Tomiałojć et al. 1984, Tomiałojć & Wesołowski 1994, 1996, Wesołowski et al. 2002, 2006, 2010, 2015).
The trends and yearly indices of the 30 most numerous bird species were estimated using ‘rtrim’ package (Bogaart et al. 2020) in R 3.6.3 (R Core Team 2019). Trim uses generalised estimating equations (GEE) to compute trends and indices, and, at the same time, accounts for overdispersion, serial correlation, and allows for missing data in the time series. The abundance index reflects changes in abundance in relation to the base year (1975 in our case), making interpretation easy and intuitive. For each species, we have chosen among the three models to balance complexity between the simplest model structure (model 1 hereafter; assumes constant trend for the whole period, usually unrealistic for longer time series) and the most complex one (model 3 hereafter; estimates separate trend parameters for each pair of the adjacent years). Model 1 usually suffices with short time series only and becomes too simple as the time series get longer. In contrast, model 3 will frequently be unnecessarily complex, since the trends in adjacent years are often so similar that they can be described with a single parameter (instead of two or more, depending on the number of years). Model 2 can be viewed as a consensus model between 1 and 3 as it represents a piecewise approach. Being a simplified version of model 3, it identifies so-called changepoints. At changepoints, a trend becomes different enough from the previous one to warrant a separate parameter estimate; all the remaining periods in the data are characterized with single parameters. Therefore, model 2 is more parsimonious than model 3, while reflecting trend changes sufficiently well. We have fitted models 1–3 to the data of all 30 species and the Akaike Information Criterion (AIC) was used for model selection. As supposed, model 2 was most often the top-supported one, while models 1 and 3 appeared too simple or too complex, respectively and much less supported. Model fit assessment has been done with the goodness-of-fit test based on Pearson chi-square statistics (Pannekoek & van Strien 2005), while overdispersion and serial correlation (both accounted for with GEE) were given for each of the top-supported models (see Appendix 1). Illustrated trends come from the top-supported models; yearly indices are point estimates with their standard errors. For each species, an overall trend (a multiplicative slope), expressed as an average annual growth rate can also be computed and is presented; it allows a classification of the species into one of the five categories originally implemented in TRIM software (strong or moderate decline, stability, strong or moderate increase, along with unclassified; Pannekoek & van Strien 2005). While packaging complex, alternating abundance changes into a single value is an obvious simplification, it is necessary when one wants to concisely characterize changes in the long-term. Despite these values can be easily interpreted, one should keep in mind, that very few species showed truly simple patterns, missing irregular or alternate changes in abundance.
We have also estimated the dynamics in species richness and trends in the total abundance over 45 years using generalized additive mixed models (GAMMs), for their flexibility in capturing non-linear patterns (Zuur et al. 2014) in the ‘mgcv’ package (Wood 2017). In these models, the response represented either the number of species or the total number of all the breeding pairs in a study plot in a given year. We assumed that response variables followed a Poisson distribution. The differences between forest types were accounted for with the fixed factor (three levels: coniferous, lime-hornbeam and riverine), while trends were estimated with separate smoothers (thin plate splines with an automated selection of k, the number of knots defining the wiggliness of the fitted curves) for each forest type (via an interaction term with forest type factor). To address the differences in plot sizes, we included log(area) as an offset term, making the model effectively estimating the number of species or the number of breeding pairs per 10 ha. Potential autocorrelation due to the same plots surveyed in subsequent years was addressed with a random ‘Plot’ intercept in both cases. Species marked with ‘+’ (breeding but < 0.5 territory within the plot) in the Tables 1–7 were also included.
Additionally, widely used community ecology indices (alpha- and beta-diversity) were applied to illustrate (a) within-plot changes in species composition and (b) between-plot similarities across 45 years. We opted for simple versions of these indices. To assess the alpha diversity we used the Gini-Simpson index (the probability that two randomly chosen individuals/breeding pairs belong to different species). To define beta diversity, we used the Sørensen index, measuring the similarity in species composition between plots, and the Bray-Curtis index, which used the abundance data and also considered the differences in species/community numbers between plots (Tuomisto 2010, Jost et al. 2011). All indices were computed in the ‘vegan’ package version 2.5-7 in R (Oksanen et al. 2020).
To estimate the direction and the rate of change in the community composition we used Time Lag Analysis (TLA) with Hellinger distances (Collins et al. 2000). So-called Hellinger distances represent Euclidean distances computed from abundance data transformed to proportions (Legendre & Gallagher 2001) at increasingly larger time-lags (e.g., year 1 to year 2, year 1 to year 3 and so on). In a linear regression, Hellinger distances were then regressed on square roots of time-lags (to minimize bias induced by smaller number of data points at larger time lags). Slopes resulting from this regression are informative in respect to the direction and rate of change a community undergoes (see Fig. 2 in Collins et al. 2000). If the community is stable, or the species vary randomly without temporal autocorrelation, distances between samples do not increase with time lags and the slope is either not different from zero or not significant. A significant positive slope (meaning distances increase with time lags) indicates instability: with time, the community veers away from its initial state (a negative slope, with distances decreasing over time means community converges to one of the early sample periods). With 45 year-long data series, there exist 990 ((n2 – n)/2, where n is the number of years) pairwise distance measures. We used Monte Carlo procedure — nonparametric bootstrapping in which columns representing years are resampled 10,000 times (see Thibault et al. 2004) — to determine variance and error probability (‘significance’ of the slope; the proportion of times among all resamples the slope was smaller than the observed slope). Year-to-year change was characterized with the Hellinger distances at lag 1 (with means and 95% CIs estimated with the Monte Carlo approach).
Breeding bird community in 2015–2019
We recorded 69 bird species breeding at least once within the study plots in 2015–2019 (Tables 1–7), including the first confirmed case of Accipiter gentilis breeding in the NW plot (Table 7). Accipiter gentilis regularly breeds in other parts of the BNP, but it has not been found nesting within the study plots previously (Wesołowski et al. 2003).
Species richness was highest in the forest-edge riverine plot K (mean 44.6, total 57 species, Tables 1, 8), followed by the riverine interior plot L (Table 2), the edge lime-hornbeam plot W (Table 3), and the coniferous plots NW and NE (Tables 6–7). The lime-hornbeam interior plots CM and MS (Tables 4–5) were least species rich (Table 8). Fringilla coelebs, Erithacus rubecula and Sylvia atricapilla were the most numerous species in all plots, almost in all years (Tables 1–7). Ficedula albicollis and Coccothraustes coccothraustes remained regular dominants (≥ 5% share) in the lime-hornbeam forest (Tables 3–5), while Phylloscopus collybita was dominant in the riverine plots K and L (Tables 1–2) and in the coniferous areas NE and NW (Tables 6–7). Sturnus vulgaris continued to be dominant in the forest-edge plots K and W (Tables 1, 3). Parus major and Turdus philomelos were common species in the lime-hornbeam (Tables 3–5) and coniferous habitats, while Phylloscopus sibilatrix remained a dominant species only in some years in the plots MS, NE and NW (Tables 5–7). The cumulative share of dominant species ranged between 43% (plot K) and 61% (plot MS) on average (Table 8).
The overall densities in the coniferous habitats (mean c. 61 p/10 ha) were by 24–25% lower than in the interior deciduous forest, and by 40–49% lower than in the forest edge plots W and K (Table 8).
Long-term changes in the number of individual bird species
Population trends of individual species fluctuated to a large degree, with many species exhibiting alternating patterns: the initial periods of population increase, followed by periods of the decline or stability (Fig. 5). Overall, 19 of 30 species were classified as increasing in numbers over 45 years (Table 9). The most spectacular increases concerned Regulus ignicapilla (by c. 5%), Dendrocopos major (3.5%), Columba palumbus (3.3%) and Sylvia atricapilla (3.1%). Few other species increased by 2% per year on average: Parus major, Phylloscopus collybita, Ficedula albicollis, Coccothraustes coccothraustes and Turdus merula.
Only six species have declined significantly, with Ficedula hypoleuca and Anthus trivialis showing the strongest declines (4–5% per year). The declines of Phylloscopus sibilatrix and Ficedula parva (2.0–2.5% per year), or Certhia familiaris and Regulus regulus (0.5–1.0% per year) were less severe.
The long-term trends of five species were classified as stable (Table 9); these included Oriolus oriolus, Lophophanes cristatus, Muscicapa striata, Prunella modularis and Fringilla coelebs.
Long-term variation in species richness and the total abundance
Species richness was stable both in riverine (K, L) and lime-hornbeam (CM, MS and W) habitats, but it clearly increased in the coniferous stands (NE, NW) over the 45 years of study (Fig. 6, Appendix 2). The largest number of species (annual average ∼43–44) bred in riverine plots, a bit smaller (∼40) in lime-hornbeam plots and initially the lowest, but steadily increasing (from ∼31 to ∼44) in coniferous plots. The alpha-diversity, expressed with the Gini-Simpson index, was approximately stable in the riverine (K, L) and lime-hornbeam stands (CM, MS, W) over the whole study period, but increased abruptly since c. 1985 in coniferous stands (NE, NW), almost levelling the values recorded in the riverine plots in the last 15 years (Fig. 7).
The breeding bird assemblage of the ash-alder forest (plot K, 33 ha). ‘+’ indicates breeding, but less than 0.5 territory; ‘–‘indicates non-breeding; bold type indicates a dominant species (comprising ≥ 5% of the community). In the species for which number ranges are given, the means were used for all further calculations.
The breeding bird assemblage of the alder swamp forest (plot L, 25 ha). For explanations see Table 1.
The total number of breeding birds (the whole community) showed relatively similar dynamics in both riverine and lime-hornbeam plots, where it clearly increased from 1985 to 2005, and then declined (Fig. 8, Appendix 3). In coniferous plots, the total abundance increased steadily but at a slower rate than in the deciduous stands (Fig. 8).
The breeding bird assemblage of the lime-hornbeam forest (plot W, 25.5 ha). For explanations see Table 1.
Similarity in the composition of bird community between the plots (beta diversity) was greatest within the same forest types (Fig. 9). Also, lime-hornbeam and ash-alder communities were more similar to each other (no matter which pair to be compared) than to the coniferous plots.
The overall resemblance of community structure across the habitats was due to low species turn-over among the areas: most species bred in all three types of forests and high repeatability of species proportions across years was evident. The dominant five species made up 45–70% of the whole community and most of them remained in this group over the whole period (Table 10). Fringilla coelebs and Erithacus rubecula belonged to dominants in all habitats, in almost every season (Table 10). Several other species, though breeding in all habitats, belonged regularly to the dominants only in some forest types: e.g. Ficedula albicollis —in both deciduous forest types, Coccothraustes coccothraustes in lime-hornbeam forest or Regulus regulus in coniferous forest (Table 10).
The breeding bird assemblage of the lime-hornbeam forest (plot CM, 24 ha). For explanations see Table 1.
Time lag analysis (TLA) indicated that despite high consistency in annual similarity of species composition within the plots (Fig. 9), the community composition changed directionally over 45 years in all plots. Species composition became less and less similar to the initial state in 1975 as the time interval increased (Fig. 10). The rates of community change (regression slopes) over 45 years were lower in the deciduous habitats (c. 0.04, Table 11) than in the coniferous plots (c. 0.07). Mean Hellinger distance at time lag 1 amounted to 0.36–0.37 in coniferous and 0.24–0.31 in the deciduous habitats implying quicker change in the former (Table 11). Though significant, these changes were rather modest in comparison with the year-to-year variability — Hellinger distances between the most different adjacent years could be of similar magnitude as the distances measured over the whole period (Fig. 10).
The breeding bird assemblage of the lime-hornbeam forest (plot MS, 30 ha). For explanations see Table 1.
The study constitutes the five-year extension of the long-term research on the dynamics of breeding avifauna in the primeval habitats of BNP (Tomialojć et al. 1984). The work provides novel information on population trends for the 30 commonest forest bird species and the directional changes in the composition of bird community over 45 years. Despite some discrepancies, the observations from 2015–2019 resemble the patterns recorded in 2010–2014 and in the last 40 years in general (Wesołowski et al. 2015); the structure of the breeding bird community remained similar and most bird species bred in BNP in all types of forests. The overall density of bird species was alike as in the previous five-year period (Wesołowski et al. 2015), but it was by c. 30% lower than around the peak year of 2005 (Wesołowski et al. 2006), and higher than in the 1970s (Tomiałojć et al. 1984). As previously, the overall density of bird species varied from the highest in the riverine plot at the forest edge to the lowest in the coniferous plots in the forest interior (Wesołowski et al. 2015).
Long-term changes in numbers of individual bird species
The numbers of individual bird species varied extensively over 45 years, although they usually stayed within the less than fivefold range (except for Regulus ignicapilla). For the majority of species, the numbers showed alternating phases, with the initial periods of growth followed by decline or stability. Steady decline or increase of bird numbers occurred rarely.
The breeding bird assemblage of the pine-bilberry coniferous forest (plot NE, 25 ha). For explanations see Table 1.
The changes of species numbers could be partially explained by natural transformation of the habitats within BNP. For example, the initial increase and then the decline of Phylloscopus collybita numbers in coniferous forest followed the appearance of numerous canopy gaps in the early 2000s and their later ‘closure’ due to natural regeneration by young, mostly deciduous trees. Also, the occurrence of additional ‘deciduous forest’ bird species, such as Ficedula albicollis, Poecile palustris or Coccothraustes coccothraustes, followed the replacement of coniferous trees by deciduous trees in the formerly conifer-dominated habitats.
However, the numbers of other bird species did not vary as predicted by habitat change. For example, the decline of the ‘spruce-dependent’ Regulus regulus appeared evident around 2000, whereas the Norway Spruce decline started in mid-20th century and continued in BNP since then (Brzeziecki et al. 2016). The recent ash dieback in ash-dominated stands might result in the reduction of tree-hole availability. The decline of Sturnus vulgaris would be therefore expected due to the limited abundance of the potential nesting sites in ashes (Wesołowski 2021), but this did not seem to happen, possibly because holes suitable for Sturnus vulgaris were plentiful in other tree species as well.
The breeding bird assemblage of the pine-bilberry coniferous forest (plot NW, 25 ha). For explanations see Table 1.
Main structural parameters of bird communities in different BNP stands in 2015–2019. Mean values ± standard deviations for individual plots are given; densities in pairs/10 ha.
As such, the numbers of breeding birds in the BNP were most likely affected also by factors operating at a wider geographical scale (see e.g. Portaccio et al. 2021). These factors may include habitat transformation or climate change along the migratory routes and/or on wintering grounds (Jones & Cresswell 2010), and explain similar population trends in BNP and nationally since 2000 (Chodkiewicz et al. 2019, Wardecki et al. 2021). For example, the increasing numbers of Columba palumbus or Dendrocopos major, the decline of Anthus trivialis or similar fluctuations of Erithacus rubecula in BNP and nationally are striking.
Long-term variation in the composition of bird breeding community
The composition of the breeding avifauna of the BNP has not changed substantially since 1975 (Tomiałojć & Wesołowski 2004, Wesołowski et al. 2003, 2015); the same bird species breeding in 1970s (Tomiałojć et al. 1984) were present recently. Also, the core of the breeding bird community in the BNP was stable, with almost 40% of species found breeding in the plots each year, and 57% of the species breeding in > 35 years. Similar pattern of the dominant species forming a core of the breeding bird community is widespread in old-growth deciduous forests across temperate Europe (review in Wesołowski et al. 2018).
Multiplicative slopes (rates of change) from trend analysis for 30 forest bird species in BNP, 1975–2019. All trends were significant and categorized as either moderate growth (↑), moderate decline (↓) or stability (↔) over 45 years. ↓
Species richness in the riverine and lime-hornbeam habitats did not show a clear long-term change, but the number of species increased in both coniferous plots (Tomiałojć et al. 1984), since large canopy gaps appeared after hurricanes. The incremental change in the composition of community structure was rather slow in deciduous plots (regression slope of c. 0.04) and greater in coniferous stands (c. 0.07). Such slow incremental changes may be typical for natural forests, but when natural disturbances are large, one would expect relatively high slope values. Kampichler et al. (2014) suggested a rate of change between 0.04 and 0.06 for disturbed or managed stands, which was similar to what we estimated for the BNP's stands. Thus, although natural, the ‘disturbances’ in the BNP’s coniferous stands, such as an appearance of large canopy gaps and their colonisation by deciduous tree species, appeared substantial as manifested by the quickest rate of change and year-to-year variability of the avian community.
Dominant (constituting > 5% of assemblage) breeding bird species in relation to habitat in BNP, 1975–2019. Only species belonging to dominants in > 25% of all season-plot combinations (sample size, n, given in parentheses) in any habitat type are shown. Riverine — plots K & L; Lime-hornbeam — plots CM, MS & W; Coniferous — plots NE & NW.
Whereas species richness (alpha diversity) of breeding bird community in the primeval old growth stands was very high, the change of the community structure between habitats (beta diversity) remained rather low. The current study confirms the conclusions by Tomiałojć and colleagues, who stated in 1984 that the same bird community inhabits all the BNP plots, and only small local deviations may occur (Tomiałojć et al. 1984). Most bird species breeding in the Białowieża NP have rather wide habitat niches and can occupy an array of structurally different places. Thus, even large disturbances, like formation of windfall gaps, do not lead to substantial species turnover, provided the forest tree composition remains similar (Fuller 2000). Such broad niches seem to be characteristic for European forest birds in general (Blondel 2018, review in Wesołowski et al. 2018), but differ from the North American birds, which show much higher habitat specialisation (Mönkkonen 1994, Blondel 2018).
Changes in the breeding assemblage composition with time in relation to forest type in BNP. Time Lag Analysis results: slopes and Hellinger distances at time lag 1 for the seven studied plots, 1975–2019. Slope measures the rate of change (the higher the value, the faster the change), p is the error probability while year-to-year variability is measured by the mean Hellinger distance at time lag 1 (see Methods).
The general compositional and structural stability of the breeding bird community observed in the primeval old-growth forests of the BNP is remarkable, as each spring most of the community is formed by different individuals, which are probably new to the area (most of species are not philopatric; review of evidence in Wesołowski et al. 2015). The breeding community largely ‘disassembles’ in autumn (less than 20% of the breeding species and less than 10% of individuals spend the winter here; Wesołowski et al. 2015) and is ‘re-created’ the following spring by newly-arriving’colonisers’. Thus, the arriving birds have to assess habitat suitability and decide on where to settle using mostly the cues available to them at the arrival time. Therefore, the high constancy of species composition of the breeding bird community in the BNP seem to result largely from the interplay of two factors: 1) long-term stability of the forest habitats (particularly the deciduous stands), causing places suitable in one year to remain so over many seasons, and 2) cross-generational reproducibility of the selection criteria used by the birds in their settlement decisions (Wesołowski et al. 2015). As long as the area retains its attractiveness, it would continue to be chosen by the same species in consecutive years. However, where the major transformations in habitat structure occur (as happened in the coniferous habitat), this may make the area unacceptable for some species and more attractive for the others. This would then result in changes of the breeding community composition. The processes of community formation by birds in BNP are, thus, driven both by the population demography of species (no birds may settle when none survives), and active acceptance of the same spots by mostly different individuals in consecutive seasons (Wesołowski et al. 2015). Similar situation might apply to other European forests that undergo little anthropogenic or natural transformations, but further studies would be needed to test this idea.
Gathering reference information on the functioning and natural dynamics of forest ecosystems, which are free of direct human impact, and studying ecological and behavioural adaptations of birds in their evolutionary (primeval) settings, is fundamental for understanding the ecology and evolution of forest birds. This served as a motivation to begin the study on the breeding bird community in the strictly protected part of the BNP in 1975 (Tomiałojć et al. 1984). After 45 years of continued study, we can also see how these pristine communities react to the large scale environmental changes occurring during this time. So far, they appeared highly adaptable and resilient to varying environmental challenges. However, no plasticity is limitless and the increase of external pressures may lead, at some point, to changes in the BNP avian communities becoming stronger and increasingly quicker as well. Provided that the current study is continued, we would be able to detect such break points. Understanding how species, communities and entire ecosystems respond to environmental change is critical to our efforts to arrest the current diversity crisis, but it can be done only with truly long-term studies (White 2019). Sadly, long-term studies are in decline often due to difficulties in ensuring stable funding (Hughes et al. 2017, Kuebbing et al. 2018). Thus, more endeavour is needed to assure the continuation of the long-term research which constitutes an invaluable source of crucial information.
Our work would not be possible to accomplish without the support of numerous people. First of all, we express our gratitude to M. Czuchra for her help in elaborating census data. We thank also L. Tomiałojć for his long lasting participation in the verification of species maps evaluations. The kind co-operation of the Białowieża National Park administration is acknowledged as well. Fieldwork and analyses were, to a large extent, supported by the internal grants from the Faculty of Biological Sciences, University of Wrocław to TW, GN and MM. Field work of DC was supported by Siedlce University.
[Długoterminowe zmiany w ugrupowaniach ptaków lęgowych w lesie pierwotnym: 45 lat cenzusów w Białowieskim Parku Narodowym (Polska)]
W pracy omówiono długoterminowe trendy zmian liczebności 30 leśnych, pospolitych ptaków lęgowych na 7 stałych powierzchniach próbnych w Białowieskim Parku Narodowym w latach 1975– 2019 (wschodnia Polska, dalej BPN, Fig. 1). Przeanalizowano zmiany w składzie ugrupowań ptaków lęgowych w trzech najważniejszych typach drzewostanów pierwotnych w BPN: łęgach (Fig. 2), grądach (Fig. 3) i borach (Fig. 4).
Cenzusy kontynuowane w latach 2015–2019 z wykorzystaniem kombinowanej metody kartograficznej na siedmiu stałych powierzchniach próbnych w BPN, monitorowanych od 1975 r., w tym dwóch na skraju lasu i pięciu w jego wnętrzu wykazały gniazdowanie 69 gatunków ptaków. Najwięcej gatunków gniazdowało w łęgach (57 gatunków na skraju lasu i 45 we wnętrzu), mniej w grądach (45 gatunków na skraju lasu, 38 i 39 gatunków we wnętrzu) i borach (44 i 50 gatunków, obie powierzchnie we wnętrzu lasu) (Tab. 1–8). Liczebność poszczególnych gatunków ptaków była mocno zmienna — większość gatunków charakteryzowała się naprzemiennie występującymi okresami istotnych wzrostów, spadków i stabilizacji (Fig. 5). Większość gatunków (19 z 30) w 45-letnim okresie badań (1975–2019) wykazywała wzrost liczebności. Najsilniejsze wzrosty (średnio 3–5% rocznie) odnotowano u grzywacza, dzięcioła dużego, kapturki i zniczka (Tab. 9). Najsilniejszy spadek liczebności dotyczył świergotka drzewnego (–4,7% rocznie) — związany był on z zarastaniem luk pohuraganowych i prawie całkowitym wycofaniem się gatunku z powierzchni próbnych. Ponadto 5 dalszych gatunków wykazywało istotne spadki: muchołówka żałobna (–4,1% rocznie), świstunka leśna (–2,7% rocznie), muchołówka mała (–2,2% rocznie), mysikrólik (–1,1% rocznie) i pełzacz leśny (–0,4% rocznie). Zięba, muchołówka szara, wilga, czubatka i pokrzywnica fluktuowały liczebnie, ale ich populacje w omawianym okresie 45 lat były stabilne. Zięba, rudzik i kapturka były najliczniejsze na większości powierzchni w prawie wszystkich latach (Tab. 10). W większości te same gatunki ptaków gniazdowały na powierzchniach w latach 1970-tych i obecnie. Bogactwo gatunkowe ugrupowań było stabilne w lasach liściastych, a borach rosło (Fig. 6), podobny wzorzec zmian dotyczył alfa-różnorodności (Fig. 7). Całkowita liczebność wszystkich gatunków w lasach liściastych rosła od 1975 i była najwyższa około 2005 roku, po czym spadła, a w borach wzrastała przez cały okres badań (Fig. 8). Wskaźnik podobieństwa ugrupowań ptaków lęgowych między powierzchniami (beta-różnorodność) nie zmienił się znacząco w okresie badań: łęgi i grądy pozostawały podobne do siebie, a bory były bardziej odmienne od lasów liściastych (Fig. 9). Ugrupowania ptaków lęgowych w borach wykazywały najszybciej postępujące zmiany składu i stawały się coraz bardziej odmienne od stanu początkowego z 1975 r., mimo że analogiczne zmiany były też widoczne w lasach liściastych, ale zachodziły wolniej (Fig. 10, Tab. 11). Szybkie zmiany ugrupowań zachodzące w borach można wiązać z zarastaniem luk pohuraganowych przez drzewa liściaste typowe dla grądów, i kolonizacją powierzchni borowych przez „grądowe” gatunki ptaków. Wskazuje to, że nawet w lasach gdzie brak jest bezpośredniej ingerencji człowieka, ugrupowania ptaków lęgowych znacząco zmieniają skład w dłuższych okresach czasu, a dynamika ta jest powodowana przez czynniki pochodzenia naturalnego.
Długoterminowe badania, jak powyższe w BPN, dostarczają cennych informacji dotyczących zmian liczebności poszczególnych gatunków oraz tempa naturalnej wymiany gatunków w ugrupowaniach ptaków w lasach pierwotnych, pozbawionych bezpośredniej ingerencji człowieka.
GEE models used for trend estimation of forest bird species in BNP, E Poland, 1975–2019. Model 1 assumes a constant trend for the whole study period (usually unrealistic for long time series), while model 3 estimates separate trend for each pair of adjacent years — in most cases it has too many parameters. Model 2 is a consensus-type model and reduces the number of parameters of model 3 by a piecewise approach: it identifies periods over which a trend can be summarized into a single paramater and only at identified changepoints between periods trend becomes sufficiently different to warrant a separate estimate. Non-convergence of a model due to insufficient number of observations is marked with ‘–‘. The top-supported model, used for inference is given in bold. Abbreviations: AIC — Akaike Information Criterion, Corr — serial correlation, OD — overdisperson (both estimated under the top supported model), LR — likelihood ratio statistic, df — degrees of freedom, p — probability.
Coefficients from a GAMM explaining long-term variation in species richness (the total number of breeding species) in relation to forest type in BNP, E Poland, 1975–2019. For the fixed terms, estimates for habitats are given as contrasts to the intercept; annual change in species richness was modelled with separate smooth terms for each habitat type.
Coefficients from a GAMM explaining long-term variation in the total abundance (summed abundances of all individual breeding species) in relation to forest type in BNP, E Poland, 1975–2019. For the fixed terms, estimates for habitats are given as contrasts to the intercept; annual change in abundance was modelled with separate smooth terms for each habitat type.