Can Topographic Variation in Climate Bu ﬀ er against Climate Change-Induced Population Declines in Northern Forest Birds?

: Increased attention is being paid to the ecological drivers and conservation measures which could mitigate climate change-induced pressures for species survival, potentially helping populations to remain in their present-day locations longer. One important bu ﬀ ering mechanism against climate change may be provided by the heterogeneity in topography and consequent local climate conditions. However, the bu ﬀ ering capacity of this topoclimate has so far been insu ﬃ ciently studied based on empirical survey data across multiple sites and species. Here, we studied whether the ﬁne-grained air temperature variation of protected areas (PAs) a ﬀ ects the population changes of declining northern forest bird species. Importantly to our study, in PAs harmful land use, such as logging, is not allowed, enabling the detection of the e ﬀ ects of temperature bu ﬀ ering, even at relatively moderate levels of topographic variation. Our survey data from 129 PAs located in the boreal zone in Finland show that the density of northern forest species was higher in topographically heterogeneous PAs than in topographically more homogeneous PAs. Moreover, local temperature variation had a signiﬁcant e ﬀ ect on the density change of northern forest birds from 1981–1999 to 2000–2017, indicating that change in bird density was generally smaller in PAs with higher topographic variation. Thus, we found a clear bu ﬀ ering e ﬀ ect stemming from the local temperature variation of PAs in the population trends of northern forest birds.


Introduction
Climate change is a major threat to biodiversity [1,2], increasingly affecting present-day species populations and communities [3][4][5][6]. Upslope shifts in the distribution of species have already been observed with rising temperatures both in temperate [7,8] and tropical regions [9,10], although the elevational shifts may vary considerably in terms of direction and magnitude among species [8]. Changes in climate are projected to cause further poleward and upward species range shifts in the future [2,11]; but, in many areas rapidly changing conditions will challenge the ability of species to track the spatial reorganization of suitable locations [12][13][14]. Thus, climate-smart conservation strategies are required to identify landscape characteristics which best support the range-shifting and persistence of species populations [15,16]. One promising strategy is the prioritization of areas with

Study Areas
Finland stretches 1100 km across the boreal biome of northern Europe (Figure 1). In this biome northern species show increasing and southern species decreasing trends towards the north. The boreal biome can be divided into southern, middle and northern boreal zones. Here, we focused on bird population trends in the middle boreal zone, where the ranges of southern and northern species largely overlap [45]. The PAs occurring in the middle boreal zone are mostly covered with coniferous (dominated by Scots pine Pinus sylvestris or Norway spruce Picea abies), mixed and deciduous (dominated mainly by birch Betula spp.) forests and pine and open, treeless mires.
Diversity 2019, 11, x FOR PEER REVIEW 3 of 12 density of northern bird species declining less in PAs with larger variation in local temperature (conditions)?

Study Areas
Finland stretches 1100 km across the boreal biome of northern Europe ( Figure 1). In this biome northern species show increasing and southern species decreasing trends towards the north. The boreal biome can be divided into southern, middle and northern boreal zones. Here, we focused on bird population trends in the middle boreal zone, where the ranges of southern and northern species largely overlap [45]. The PAs occurring in the middle boreal zone are mostly covered with coniferous (dominated by Scots pine Pinus sylvestris or Norway spruce Picea abies), mixed and deciduous (dominated mainly by birch Betula spp.) forests and pine and open, treeless mires. Figure 1. Location of the study areas (protected areas (PAs) in black lines) in Finland with the mean air temperature for April-May-June at 50 m × 50 m grid size. As an example, the topographic variation (meters above sea level; DEM = digital elevation model) and mean April-June air temperatures are shown for a specific region with a high variation in temperature caused by variation in topographic heterogeneity. Location of the study areas (protected areas (PAs) in black lines) in Finland with the mean air temperature for April-May-June at 50 m × 50 m grid size. As an example, the topographic variation (meters above sea level; DEM = digital elevation model) and mean April-June air temperatures are shown for a specific region with a high variation in temperature caused by variation in topographic heterogeneity. In this study, birds were counted in 129 PAs situated in the middle boreal zone both in 1981-1999 and in 2000-2017 (including the southernmost part of northern boreal zone; uniform coordinate system [epsg: 2393], grids 70-74, Figure 1). In these PAs the total land area was 4347 km 2 , with a mean of 33.7 km 2 (size range 1.6-300.3 km 2 , Figure 1, Table S1). The forests (including wooded mires) covered 63% of the land area of the studied PAs, the remaining being open, treeless mires. This information on the proportion of forests in the studied PAs is based on the data of habitat classification in National Parks Finland, which governs state-owned protected areas. More than two thirds of the protected forest stands are over 100 years old [46].

Bird Censuses
Land birds in the studied 129 PAs were counted using the Finnish line transect census method, which is suitable for counting birds over large areas [45,47]. The line transect method applies a one-visit census in which birds are counted during the breeding season along a transect with an average length of 5-6 km.
The census is carried out in late May and in June in the early morning, when the singing activity of the birds is highest. In the line transect method, a 50 m (25 m on each side)-wide main belt along the walking line and a supplementary belt outside the main belt are separated. Birds are counted both from the main belt and the supplementary belt which together constitute a survey belt. The supplementary belt consists of all birds observed outside the main belt (for details, see [45,47,48]).
In the Finnish line transect method, the densities of species based on the observations in the censuses are calculated in standard units of pairs/km 2 . A pair is inferred (following the instructions for the Finnish line transect census) by: a male heard singing, an otherwise observed male or female, or a group of observed fledglings. Densities of bird species were calculated on the basis of observations in the whole survey belt. In calculating the density, a species-specific correction coefficient (based on the proportion of main belt observations in the whole survey belt) was used to take into account the differing audibility and other detectability of different species (see, e.g., [47,48]).
When counting the birds, in both time periods the total length of the transects in each of the 129 PAs was at least 1.0 km (see Figure 1, Table S1). The total length of the line transect censuses in the studied PAs amounted to 4677. 8 (Table S1). The median number of years when the censuses were carried out was two in each protected area, both in 1981-1999 and in 2000-2017. The median census year was 1994 in the first period and 2009 in the second period. Precisely the same transects were not surveyed, but the censuses in each protected area included the same proportion of habitats during the two periods.
The population changes of northern forest bird species revealed by the census data from the two time periods were related to the topoclimatic variation (50 m × 50 m) within the PA, the changes in the broad-scale (10 km × 10 km) climate, and the size and forest proportion of the PA. We concentrated only on forest species because open mire species-such as many waders-occur on flat peatlands and thus are inherently associated with breeding sites with a low variation in local climatic conditions. We included all bird species breeding in forests with a northern distribution in Finland. Based on these criteria, the following 17 species [49,50] were selected for our study: rough-legged buzzard Buteo lagopus, golden eagle Aquila chrysaetos, merlin Falco columbarius, northern hawk owl Surnia ulula, great grey owl Strix nebulosa, three-toed woodpecker Picoides tridactylus, Bohemian waxwing Bombycilla garrulus, red-flanked bluetail Tarsiger cyanurus, Arctic warbler Phylloscopus borealis, Siberian tit Poecile cinctus, great grey shrike Lanius excubitor, Siberian jay Perisoreus infaustus, brambling Fringilla montifringilla, common redpoll Carduelis flammea, two-barred crossbill Loxia leucoptera, pine grosbeak Pinicola enucleator, and rustic bunting Emberiza rustica. Over 80% of recorded bird pairs of the 17 species studied belong to the three most abundant species, i.e., the brambling, common redpoll and rustic bunting (see Table 1).

Climate Data
Our climate data consisted of data recorded at two distinctly different spatial resolutions. First, as a measure to examine the effect of the broad-scale air temperature in the two time periods (i.e., 1981-1999 and 2000-2017), as well as the changes between the two periods, on the bird species densities, we used the mean air temperature values of April-June (TAMJ broad ) interpolated for each of the 10 × 10 km grid squares where a protected area was located. We focus on April-June air temperatures, because this is the time of arrival of migratory birds and the breeding period for all the studied northern bird species. The temperature values were obtained from the daily gridded climate dataset by the Finnish Meteorological Institute [51] both for the period of 1981-1999 and 2000-2017.
Second, in order to capture the variation in the local climatic conditions in the studied PAs, we modelled the monthly average air temperatures (1981-2010) across the study domain based on daily data from 318 meteorological stations derived from the European Climate Assessment and Dataset database (ECA&D; [52]). For these models, we employed generalized additive models (GAM), as implemented in R-package mgcv version 1.8-7 [53]. Following Aalto et al. [54], we used predictors of the geographical location, topography (elevation, potential radiation, relative elevation) and water cover [55] to produce monthly average air temperature surfaces at a high spatial resolution (50 m × 50 m) for all months. The modelled monthly average air temperatures agreed well with the observations from the meteorological stations, with a root mean squared error (RMSE; leave-one-out cross-validation) ranging from 0.56 (June; Figure S1) to 1.49 • C (January). We used the fine-grained (50 m) variation within the PAs (σTAMJ local , quantified as the standard deviation of April-June mean temperatures across PAs) in the subsequent analyses, as an indicator of local temperature variation driven by topographic heterogeneity, and thereby an indicator of potential buffering effect.
We could not account for the effect of forest canopy cover in our modeling of local temperatures. This is due to the placing of meteorological stations on the open landscape [56]. Although this is likely to increase modelling uncertainty in forested areas, the microclimatic effect of forest canopy (i.e., buffering temperature variability) is indirectly accounted for in the models by using the forest proportion predictor.

Statistical Analyses
We modeled the spatial variation in the population density of northern forest birds over 1981-2017 (i.e., our first response variable; combined density in 1981-1999 and in 2000-2017 in Table S1) as a function of the TAMJ broad and σTAMJ local . In addition, we used forest proportion and the size of the PA as predictors in our models. Large PAs are likely to have higher bird densities compared to small PAs [e.g., 44]. The size of the PA was log-transformed for the analyses to approximate normality.
In studying the bird density variation between the PAs, we used generalized linear mixed modelling (GLMM [57]; assuming negative-binomial errors with a log-link function), that accounts for zero-inflation in the response variable. In addition to the main effects of σTAMJ local and forest proportion, we included their interaction with the model. This interaction was included because we expected that forest cover can mediate the effect of local temperature variability on bird densities. Indeed, σTAMJ local and forest proportion were intercorrelated (Spearman's rank correlation r s = 0.706). Species were considered as a random factor in the model, and the two time-periods were included as an ordered factor (1 = 1981-1999 The modeling was carried out in R (version 3.6.1 [58]) and the GLMMs were fitted using the package glmmTMB [57]. There was a slight positive correlation between ∆TAMJ broad and σTAMJ local (r s = 0.388), σTAMJ local and size of the PA (r s = 0.312), and ∆TAMJ broad and size of the PA (r s = 0.208).
When comparing species-specific density changes in the PAs between the time periods, we used a non-parametric Wilcoxon signed rank test.

Results
We found that the combined density of the 17 studied northern forest bird species declined significantly (p < 0.001) by ca. 38% between 1981-1999 and 2000-2017 (Table 1). Five species-including all the three most abundant northern forest bird species (brambling, redpoll, rustic bunting)-and two less-abundant species (Arctic warbler, two-barred crossbill) declined significantly while three species increased (three-toed woodpecker, Bohemian waxwing, red-flanked bluetail).
Our modelling showed that broad-scale air temperature had a negative effect, and local temperature variability a positive effect, on the densities of northern forest bird species (Table 2). On average, densities of northern forest bird species increase northwards (Table S1). Therefore, the variation of TAMJ broad was negatively related to the bird abundances as values of April-June average air temperatures generally decline northwards. In addition, we found a significant (p < 0.05) negative interaction between forest and σTAMJ local , indicating that an increase in forest proportion inside PAs tends to decrease the effect of local temperature variability on bird densities. Table 2. Results of generalized linear mixed modelling (GLMM; assuming negative-binomial errors with a log-link function) of bird population densities in PAs. Bird densities were modelled using broad-scale air temperature variation (TAMJ broad ), local air temperature variability (σTAMJ local ), forest proportion, interaction between σTAMJ local and forest proportion and PA size as predictors. Species were considered in the models as random factors. Predictor Time shows the division of data into two distinct parts. The protected area size was log-transformed prior to the analysis. Local temperature variability had a significant (p < 0.001) negative effect on the density change of northern forest birds from 1981-1999 to 2000-2017 (Table 3). Importantly, the decline in bird density was smaller in PAs associated with large temperature variation. In contrast, change in broad-scale air temperatures (∆TAMJ broad ) did not have any significant effect on the density change of northern forest birds. The size of the PAs had a positive effect on the change of population densities of northern forest birds, so that densities declined more in large PAs. Table 3. Results of generalized linear mixed modelling (GLMM; assuming Gaussian errors) of bird density change in PAs. Bird density was modelled using change in broad-scale air temperatures (∆TAMJ broad ), local air temperature variability (σTAMJ local ) and PA size as predictors. Forest proportion was an offset variable in the analysis. The protected area size was log-transformed prior to the analysis.

Discussion
Our analyses show that local temperature variability had a positive effect on the population density of northern forest species in the time periods of 1981-1999 and 2000-2017. High temperature variation also seems to buffer the negative decline of northern forest birds. Our results are thus in line with earlier studies investigating the impacts of topoclimatic variation on the persistence of species populations that are based, as our study is, on data from repeated surveys [27,30,59].
Interestingly, broad-scale temperature affects the density patterns of northern forest birds, but change in temperature does not explain the decrease of these species in the model. This can be due to the fact that in the model local temperature variability outweighs the negative effects of climate warming on the northern declining bird species. Moreover, the range of change in broad-scale temperatures is much less than variation in local temperatures, so broad-scale temperatures probably affect more similarly everywhere. The negative effects of climate warming on the northern bird species have been observed in earlier studies with their ranges shifting northwards [50,60], population densities in PAs declining [61] and mean weighted densities moving northwards [47].
Thus, bird densities tend to be higher in topographically more heterogeneous PAs, and the local temperature variability seems to provide a buffering effect against the declining population trends of northern forest birds in our study area. High-latitude regions, such as our study region, have experienced more severe warming than other areas-on average twice the global rate of warming [62]. In our study region, the mean April-June temperature has increased by an average of 1.0 • C, and the mean annual temperature has increased by an average of 1.1 • C between the periods 1981-1999 and 2000-2017 [40]. Mirroring our results against these trends in climate, detection of workable mechanisms for retarding the negative effects of climate warming can be considered as of special importance.
Controlling for the potential land use-induced impacts on populations from climate change impacts in coarse-scale grids can be difficult, particularly as human interventions can cause both negative and positive impacts on species viability (e.g., loss vs. restoration of habitat). Previous studies have not been carried out in protected areas and, therefore, land-use change may bias the observed results, because in a flatter landscape land use is much more intensive than in a rough terrain or in mountains. In our PAs, human-caused habitat alteration is not allowed; therefore, our results are more explicitly related to the topographic heterogeneity and climate-induced changes than in earlier studies.
As regards the extend of vertical variation, our results are highly interesting. This is because topographic variation of our PAs is relatively moderate, with a maximum vertical within-PA difference of about 250 meters. Yet, this variation, together with differences in incoming solar radiation between North-and South-facing slopes [54], is sufficient to counteract the recent macroclimatic forcing impacts on bird populations. Thus, the buffering impacts of local air temperatures may not require very large elevational gradients but can be effective already in moderately varying landscape. A study by Gaüzère, Prince and Devictor [59] investigated the climatic debt in bird communities in France using the species community temperature index (CTI), and showed that an increasing range in elevation can affect bird species populations positively by reducing the rate of their change. However, the study area in France also included mountainous landscapes, where the range of elevational variation was measured in several hundreds of meters compared to our study PAs with the variation commonly measured in tens of meters.
A central difference between the present results and a study carried out in the UK [27] is that the latter was based on presence records recorded in Atlas data at a 10-km resolution, while our work examined species abundances and the buffering capacity of fine-grain temperature patterns at a 50 × 50 m resolution. In the UK study, the extirpation risks stemming from climate warming were reduced by 22% in plants and 9% in insects due to topographic variation creating potential for microrefugia [27]. However, assessing changes in species trends based on occurrences from such coarse-resolution Atlas data may only partly conceal the changes in species populations which would otherwise be visible in abundance trends [33,47,63]. Thus, by using abundance analyses, the effect of topographic heterogeneity might be even more clearly observable.
The benefits of topographically more heterogeneous PAs may also manifest themselves in single species ecological and behavioral phenomena; for example, by helping a species to avoid misusing phenological triggers and the associated temporal mismatches between trophic levels [64]. In the boreal region, in topographically heterogeneous southern Norway, the breeding success of black grouse Tetrao tetrix and the capercaillie Tetrao urogallus has been observed to increase under the warming climate [65]. In contrast, in the flatter boreal forests of Finland, populations of black grouse have been declining-which has been attributed to the mismatch between the advanced time of mating and chicks hatching too early, resulting in declining breeding success [66].
Our findings suggest that topoclimatic variation may provide an effective conservation measure against climate change in moderately rugged landscapes. Such areas can be common in many places in the western boreal Palaearctic, where landscape is covered by a peneplain with a moderately rough land surface with small-scale variability in topography produced by a long period of erosion and other physical processes [67]. In these vast regions in the boreal biome, topographical heterogeneity may provide a significant buffer against rapid climate-induced changes in many areas [19].
The results reported here also have important consequences for climate-smart conservation planning for boreal regions. In addition to support for protecting particularly topographically heterogeneous large areas for declining northern species, another tentative consideration and a potentially useful tactic in climate-smart conservation planning might be targeting for a sufficiently Diversity 2020, 12, 56 9 of 12 connected network of topoclimatically heterogeneous PAs. This kind of approach could provide climatically suitable stepping stones and holdouts for a population to persist for a limited period of time, thereby facilitating the range shifts of contracting species under deteriorating climatic conditions (cf. [20]).

Conclusions
Our study was carried out in protected areas where harmful land use, such as logging, is not allowed; therefore, our results of the positive effects of topographical heterogeneity on the populations of declining northern birds are rather straightforward. In contrast, in comparisons of lower topographical variability (e.g., flatlands) and higher topographic variability (rough terrain, mountains), land-use patterns may differ considerably and, thus, the lower land-use intensity may affect the buffering capacity of rough terrains against climate change. Further studies are also needed to demonstrate the significance of topographical heterogeneity in preserving biodiversity in a rapidly warming boreal climate outside protected areas with varying intensities of land use.

Funding:
The study was part of a project evaluating the protected area network in the changing climate (SUMI) which was funded by the Ministry of the Environment in Finland and was financially supported also by the Strategic Research Council (SRC) at the Academy of Finland (IBC-Carbon, no 312559).