Our data originated from 29 villages and towns (hereafter, we use the term “towns” for simplicity) in northern coniferous forest biome in Finland (Jokimäki et al., 1996; Suhonen et al., 2009, 2010). Most of the study plots were about 30 ha (mean ± SD: 31.2 ± 7.7 ha; range: 13–51 ha). The study areas were situated relatively evenly along a 950-km north–south extent (60–68°N; Jokimäki et al., 1996). The study areas included urban settlements, ranging from large towns to villages; both habitat types were located evenly across Finland. The human populations in these communities ranged from 300 to 159,000 persons (mean ± SD: 21,694 ± 34,018 inhabitants) (Jokimäki et al., 1996). The variation in local habitat structure and edge effects were minimized by selecting each study plot from the most heavily urbanized area of each site. All study plots contained houses, roads, and scattered green areas (woods, parks, and fields) between houses (Jokimäki et al., 1996). The proportion (%) of block-of-flats, single-family house areas and wooded areas (mean ± SD) in the study plots has not changed during the study period [block-of-flats: 32.97 ± 39.40 (1990) and 32.86 ± 39.86 (2014); single-family houses: 33.31 ± 29.98 (1990) and 33.41 ± 30.05 (2014); green areas: 33.24 ± 27.01 (1990) and 33.24 ± 27.01 (2014); Mann–Whitney U-test; P > 0.05 in all cases].

The mean mid-winter temperature (January: mean ± SD) of the study sites differed between study winters, being lowest during 2009–2010 (−11.55°C ± 0.25), intermediate during 1999–2000 (−7.42°C ± 0.25), and highest during 1990–1991 (−4.7°C ± 0.25; F = 184.69; df = 2, 56; P < 0.001; Finnish Meteorological database). There was more snow (cm; mean ± SD) at the study sites during the winter of 1999–2000 (30.80 ± 1.86) than during 1991–1992 (20.15 ± 1.86) or 2009–2010 (19.69 ± 1.96; F = 11.38, df = 2, 56; P < 0.001; Finnish Meteorological database).

Winter birds were surveyed by the single-visit study plot method (Bibby et al., 1992). Most of these censuses were made in late December or early January. We counted the birds during three different winters within 8- and 10-years interval: 1991–1992, 1999–2000, and 2009–2010. The census work in each case was done by one person (mainly the same person during each period) in good (windless and no rain) weather conditions at mid-days, when the light availability allowed the surveys to be conducted (in our case between 10.00 and 14.00). Most of the surveys were done during working days, Monday through Friday. The surveys were not done along a single route through the plots but instead as a zigzag walk through the plots. Survey speed was ~30 ha/h. To try to find all the species and individuals from the study plot, we searched birds also from the backyards and gardens of houses during surveys. This kind of a transect count reduces many of the problems associated with counting birds in urban areas, e.g., varying noise level and visibility (Jokimäki and Suhonen, 1998). It is also important to note that wintertime territorial behavior by singing in northern latitudes is lacking. Thus, singing behavior did not affect species detectability in our study. During the winter season, at least in northern areas as in our case, wintering birds are mainly foraging on bird feeders, partly due to the low temperatures, a short daytime (4–6 h) and snow cover (Hurme, 1973). Because feeding sites are more abundant in urban than in more natural sites, most overwintering birds are concentrated in these sites within human settlements (Jokimäki and Suhonen, 1998; Tryjanowski et al., 2015). Also, a scare vegetation cover, due to the absence of leaves on deciduous trees and shrubs in winter, increases the visibility of birds. Therefore, we suppose that the accuracy of the surveys is high. We also assume that the probability of observing birds is about the same for each data that are compared. In general, a single-visit census during wintering enables the detection of ~90% of the species and 80% of the individuals in urban environments (Suhonen, 1984; Jokimäki and Suhonen, 1998). To avoid double counts of individuals, we used a relatively high census rate (10 ha/20 min). All birds were included in the survey, except for overflying individuals, which did not land and stay on the study plot. We excluded aquatic bird species such as the mallard (Anas platyrhynchos), the herring gull (Larus argentatus), and the common gull (Larus canus) from these analysis, since their occurrence is mainly restricted to open water. We also excluded one large flock (over 800 individuals) of the Bohemian waxwing (Bombycilla garrulus) from the Saarenkylä (Jokimäki et al., 1996; Suhonen et al., 2009) survey site during the first winter 1991–1992, because the flock was migratory.

Following the suggestion of McGeoch and Gaston (2002), bird species were defined as core species if they occupied ≥90% of the study sites and satellite species if they occupied ≤10% of the study sites; the rest of the species were defined as intermediate species. Migratory status of bird species was extracted from Väisänen et al. (1998). We included invasive bird species in the partial migratory species category, and, in our case, migratory species included only species migrating to Central Europe or the Mediterranean area since long-distance migratory species were not detected in our study sites during winters.

We tested for the modality of the occupancy frequency distribution in each winter of study. We used 10% occupancy classes and number or percentage of wintering bird species in each class to represent the geographical variation in the SOFD, as recommended by McGeoch and Gaston (2002). We applied the nonlinear regression analysis to the empirically ranked species occupancy curves (RSOCs; Jenkins, 2011). We made separate SOFD using by the pooled data, three winters separately, separate data of large and small towns as well as towns located in the south and north. For the regional SOFD pattern analysis, we grouped study sites for southern (<63°, 30′N) and northern towns (>63°, 30′N). For the town-size analysis, we grouped study sites for large-sized (>500 inhabitants/km²) and small-sized towns (<500 inhabitants/km²; see Jokimäki et al., 1996). For these analyses, in each of the data sets, species were listed in rows and towns as columns, and data were processed separately (Jenkins, 2011). All analyses described below are based on relative species occupancy (presence/absence) data at individual towns. First, we calculated the sum of all towns where a given species was observed. We then divided each species occupancy by the total number of towns to get the relative proportion where each species occupied (Jenkins, 2011). Second, we sorted the species by their relative occupancy values in decreasing order, where R_i is the rank value for species i. We plotted the relative occupancy of the species (O_i) as a function of R_i, resulting in the RSOC. Third, we compared which of the most common core–satellite species patterns (unimodal satellite dominant, bimodal symmetric, or bimodal asymmetric) gave the best fit for the assemblages (Jenkins, 2011).

We fitted each of the following three SOFD patterns:

Unimodal satellite dominant (exponential concave): Oi = y₀ + a^*exp(–bRi), where the initial parameters were y₀ = 0.01, a = 1.0, and b = 0.01.

where y₀, a, b, and c were the estimated parameters.

The Levenberg–Marquardt algorithm (999 iterations) was used in the nonlinear regressions analyses following suggestions by Jenkins (2011), and parameters were estimated by means of ordinary least squares (OLS). We graphically evaluated the assumptions of the regressions for normality of residuals, homogeneity of variance, independent error terms, as well as the tails and shoulders of the data and models.

We used the AICc (Akaike information criterion for small sample sizes) to compare the alternative models. The model with the smallest AICc is considered to be the best with respect to the expected Kullback–Leibler information (Burnham and Anderson, 2000). The approach is powerful to detect differences between models if ΔAICc (= AICcmin-AICci) values are higher than 4 (Anderson et al., 2000; Jenkins, 2011).

We then used linear regression analysis to examine the relationships between mean local abundance and regional occupancy frequency, with the former as the independent variable and the latter as the dependent variable. Both mean local abundance and regional occupancy were log₁₀ transformed to meet the assumptions of regression analysis. We used the general linear model to test the temporal variation in the strength and intensity in slope in species occupancy frequencies (log₁₀ transformed). The explanatory variables were study winter (winter), size of town (small and large), and location of town (south and north), and abundance (log₁₀ transformed) (continuous covariate) and interaction term winter × abundance. All the data analyses were performed using the IBM SPSS statistical package, version 23.

We observed a total of 35 wintering bird species in our study towns in Finland, which is ~30% of the total number of bird species wintering in Finland (Lehikoinen and Väisänen, 2014). These results indicate that residential areas are important areas for the many overwintering bird species (Jokimäki and Kaisanlahti-Jokimäki, 2012; Leveau and Leveau, 2012). The proportion of satellite species (43%) in the pooled data was much greater than that observed earlier in a breeding season study conducted in urban core zone areas in Europe (20%; Jokimäki et al., 2016). In addition, southern towns had slightly more satellite species than northern towns, probably due to a greater number of overwintering migratory species occupying in the south than in the north (Lehikoinen and Väisänen, 2014). These results indicate that there are some seasonal differences in bird community composition in urban areas.

We found that the same three sedentary core bird species dominated the winter bird assemblages during every study winters independently of the location and the size of the town. These species (great tit, magpie, and house sparrow) belong also to the species that could be found in most winter feeding sites in Finland (Väisänen, 2018), urban wintering bird assemblages in Poland (Tryjanowski et al., 2015), as well as in breeding bird communities of most European towns (Jokimäki et al., 2016). We did not detect any between-winter cases when a core species changed to satellite species, or vice versa. In practice, there is little empirical evidence that species move between the core and satellite classes in different areas, although there are some observations that a core species may be the same across a variety of habitat types (McGeoch and Gaston, 2002). In our case, two intermediate bird species (blue tit Cyanistes caeruelus and hooded crow Corvus corone cornix) were included in the core species list in some study winters, but no satellite species were changed into the intermediate species group. However, some bird species with intermediate occupancy frequencies showed changes in their relative occupancy rate between the study winters. For example, in two coniferous forest species, the willow tit and the crested tit, and a deciduous forest species, common redpoll, occupation frequencies decreased, whereas the relative occupation frequencies of the tree sparrow and blackbird increased during the study period. The occupancy rates of frugivorous birds such as the bohemian waxwing and fieldfare varied between years due to variation on fruit crop size of rowan berry trees (Sorbus aucuparia; Lehikoinen et al., 2010; Suhonen and Jokimäki, 2015). Our results in species-specific occupancy frequencies and changes in occupancy frequencies are in a good concordance with the long-term winter bird population changes observed both in the Finnish winter bird surveys 1987–2014 (Lehikoinen and Väisänen, 2014; Lehikoinen et al., 2016) and winter feeding site study during 1988–2017 (Väisänen, 2018). Both the satellite and intermediate species have several partial migratory or migratory bird species (Table 1). It is possible that these overwintering migratory bird species might increase species richness and total abundance birds in the winter bird assemblage, especially in the south, and also increase the number of satellite species in the SOFD pattern and decrease the slope in the SAOR pattern in the future.

Both SOFD and SAOR patterns may possibly depend on sampling artifacts as well as environmental and biological mechanisms such as (i) habitat disturbance, (ii) niche breadth, (iii) dispersal ability, and (iv) sampling site position within geographical range (see Gaston et al., 2000; McGeoch and Gaston, 2002). The sample unit should be kept constant when studying assemblages with a fewer than 20 species (we had at least 22 species), and surveys should be conducted using similar extents. To facilitate comparison, our data came from the core area of each town, and standard sampling procedures were used. The sampling artifact may cause a positive SAOR if locally uncommon species are not detected, although they actually occur at a site (Gaston et al., 1997b, 2000). It is, thus, likely that some locally uncommon bird species were not detected at every site where they actually occurred. The problem of undersampling rare species is to some extent unavoidable (McGeoch and Gaston, 2002), and the longer sampling continues, the more rare species are likely to be added to the satellite mode (McGeoch and Gaston, 2002). To avoid these problems, we standardize our bird counting time to 20 min per 10 ha. Also, samples from large areas are more heterogeneous than samples from smaller areas (McGeoch and Gaston, 2002; Heatherly et al., 2007), but our study sites were almost equal in size (about 30 ha), and we counted birds at each town center to decrease habitat heterogeneity between study sites. Moreover, in most of the cases, the same person surveyed exactly the same study sites in each study winters, which decreases possible interobserver bias related to bird detectability and searching intensity. Thus, the possible artefacts of study size and bird counting methods were probably marginal in this study.

We found a bimodal symmetric SOFD pattern that was constant during the study winters. A bimodal SOFD pattern has also been reported in many earlier studies conducted in more natural habitats (Mehranvar and Jackson, 2001). We did not find any differences in SOFD patterns between different sizes of towns and towns located in different latitudes. These results conflict with the breeding season data from the core zones areas of European towns, in which a bimodality SOFD pattern was observed only within large-sized and central European towns (Jokimäki et al., 2016). These differences might be related to both the different extent of these studies (national vs. European level) and season (winter vs. summer). Also, our winter season study did not include very large sized tows with millions of inhabitants. As the scale expand, it is possible that the sample will contain more satellite species. During winter, the regional species pool is smaller and many species are concentrated in urban habitats, at least in northern latitudes like in Finland.

Some earlier studies have found some variability in the SOFD pattern among regions (Heatherly et al., 2007; Jokimäki et al., 2016; Korkeamäki et al., 2018). Our results suggest that the SOFD pattern slightly changed with the latitude, southern towns having more satellite species than the northern towns, probably due to a higher number of migratory species in the south. However, in a breeding bird study, northern regions had more satellite species than southern regions (Jokimäki et al., 2016).

We found that all core species were sedentary generalist species, whereas most (about 70%) satellite and intermediate species were partial migratory or migratory species (Table 1). The NBM did not fit well with our data because it predict that most species are rare, and thus, the SOFD pattern should be a unimodal satellite dominant mode. An earlier study has indicated that local species richness either slightly increased or decreased between winters (Suhonen et al., 2009) as predicted by the MPDM. The temporal stability of the SOFD fits well with the predictions of the MPDM, which predict co-occurrence of many core (common) and satellite (rare) species in the assemblage. If the SOFD pattern is bimodal and species relative occupation rates vary temporally, then MPDM would be shown a better fit to the data than the NBM (Hanski, 1982; Collins and Glenn, 1991). We found that in some species, the relative species occupancy frequency either increased or decreased during the 18-year study period as predicted by the MPDM, although we did not find any case wherein a core species shifted to satellite species and vice versa.

Our data indicated that the temporal variability of the SAOR was rather stable over the years as has been observed in several previous studies conducted in more natural habitats (Blackburn et al., 1998; Webb et al., 2007; Heino, 2008). In addition, the SAOR pattern did not vary between different sizes of towns or towns located in different latitudes. Although both the NBM (e.g., Brown, 1984) and the MPDM (e.g., Hanski, 1982) predict a positive SAOR, our results give more support to the MPDM (Hanski, 1982), which bimodality should result from stochastic species colonization and extinctions. In other words, even common species will go extinct at several plots and become rare (Suhonen and Jokimäki, submitted). On the other hand, rare species can become common via new arrivals (as we observed, e.g., in the case of the great woodpecker and the blackbird), although such shifts should be primarily seen only in long-term studies (e.g., Hanski, 1982). An alternative NBM predicts temporally stable occupancy and abundance patterns for species (Brown, 1984), and our data did not support that.

The handling editor declared a past co-authorship with the authors.