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

**Raimo Virkkala 1,\*, Juha Aalto 2,3, Risto K. Heikkinen 1, Ari Rajasärkkä 4, Saija Kuusela 1, Niko Leikola <sup>1</sup> and Miska Luoto <sup>3</sup>**


Received: 19 December 2019; Accepted: 31 January 2020; Published: 1 February 2020

**Abstract:** 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 buffering mechanism against climate change may be provided by the heterogeneity in topography and consequent local climate conditions. However, the buffering capacity of this topoclimate has so far been insufficiently studied based on empirical survey data across multiple sites and species. Here, we studied whether the fine-grained air temperature variation of protected areas (PAs) affects 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 effects of temperature buffering, 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 significant effect 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 buffering effect stemming from the local temperature variation of PAs in the population trends of northern forest birds.

**Keywords:** boreal; buffering; climate change; forest bird; macroclimate; population decline; protected areas; topographic heterogeneity

#### **1. Introduction**

Climate change is a major threat to biodiversity [1,2], increasingly affecting present-day species populations and communities [3–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–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

notable fine-grain heterogeneity in local climates caused by variation in topography (i.e., topoclimatic heterogeneity) [16–19]. The conservation of landscapes with high topoclimatic variation allows species to move shorter distances while tracking suitable conditions. Thereby, it may facilitate the extended persistence of trailing-edge (the contracting or retreating edge of range) populations by providing refuges and holdouts with favorable local climate conditions [20–23]. Thus, variation in topoclimates provides a potentially important buffering mechanism against the impacts of climate change on biodiversity [17].

However, although theimportance of topoclimatic variation for biodiversity preservationisintuitively credible, support for it largely comes from theoretical and model-based studies, or climatological assessments of variability in local climatic conditions [21,24–26]. In essence, the evidence on the benefits of topoclimatic variation emerging from species surveys conducted across multiple sites or protected areas (PAs) is still scarce (see, however, [19]). Monitoring studies of species trends in topographically diverse landscapes both in PAs and outside them have revealed contrasting trends, including upslope as well as downslope expansion of populations, increases in abundance as well as declines of species, and even local extinctions [8,22,24]. Repeated surveys of only one or a few PAs provide useful insights but do not enable systematic comparisons of biodiversity changes across topoclimatic gradients.

Surprisingly, very few studies have specifically investigated the potential of topoclimatic buffering to support species persistence based on repeated surveys of multiple populations. Using a large set of climate-threatened species from the UK, Suggitt et al. [27] showed that fine-grained climatic variation caused by topographic heterogeneity may indeed buffer species against regional extirpations. These results are encouraging, but it should be noted that studying changes in species occurrences over coarse scale grid cells (10 km × 10 km) brings two challenges: separating potential land use-induced impacts on species populations (cf. [25,28]) from climate change impacts, and the potential omission of declines in species abundances which often become apparent earlier than changes in occupancies [29]. In another study, Maclean, et al. [30] examined the impacts of fine-resolution microclimatic buffering on grassland plant species responses and showed that cooler slopes can support the prolonged persistence of populations. What remains open to debate is how well these results apply to other species groups in other environments showing more rapid responses to climate change than plants which have many means to survive during unfavorable periods [31].

Clearly, we need more studies on changes in the abundance of species showing potentially rapid responses to climate change. It would be especially valuable to examine changes in species abundances systematically across PAs, with different potentials for topoclimatic buffering where human-induced land use changes are not allowed. Among different species groups, birds can be particularly useful in studying the impacts of topoclimatic variation. This is because bird populations have the potential to respond rapidly to the changes in climate, as shown by the recent changes in their distributions and in the composition of local bird communities [26,32]. Moreover, there are sufficient data on bird species abundances available from systematic bird surveys [29,33,34].

Here, we studied the temporal changes in bird species densities in 129 individual PAs with varying levels of topoclimatic heterogeneity, using systematically recorded data from a national PA network which may effectively alleviate the negative impacts of climate change on species [35–38]. We investigated the variation within the PAs in fine-grained air temperature patterns to assess whether larger temperature variation can slow down the negative effects of climate warming on bird species populations. We focused specifically on population size changes in northerly distributed bird species in a region where these species have contracting range-margin populations [39,40]. Alongside the temperature variability, we addressed two other potentially important drivers of bird population sizes: changes in the observed broad-scale air temperature (i.e., macroclimate) [41,42] and the size of the PAs [43,44], and also took the variation in forest proportion in the PAs into account.

Specifically, we addressed the following two questions: (1) Is the density of bird species populations higher in PAs showing a larger variation in local temperature (conditions)? (2) Is the density of northern bird species declining less in PAs with larger variation in local temperature (conditions)?

#### **2. Material and Methods**

#### *2.1. 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.

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 km2, with a mean of 33.7 km<sup>2</sup> (size range 1.6–300.3 km2, 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].

#### *2.2. 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/km2. 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 km in 1981–1999, and 4585.9 km in 2000–2017, and the mean total length of the transects in a PA was 36.3 km in 1981–1999 and 35.5 km in 2000–2017 (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).



#### *2.3. 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 (*TAMJbroad*) 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 (σ*TAMJlocal*, 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.

#### *2.4. 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 *TAMJbroad* and σ*TAMJlocal*. 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 σ*TAMJlocal* 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, σ*TAMJlocal* and forest proportion were intercorrelated (Spearman's rank correlation rs = 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, 2 = 2000–2017). The model structure in R syntax:

$$\begin{array}{c} \text{Bird density} = \, TAMf\_{\text{Ironal}} + \sigma TAMf\_{\text{Ironal}} + \text{forest proportion} + \sigma TAMf\_{\text{Ironal}} \times \text{forest} \\ \text{proportion} + \log(\text{PA size}) + \text{time} + (1 \mid \text{species}) \end{array} \tag{1}$$

Secondly, we modelled the change in bird species density from 1981–1999 to 2000–2017 (i.e., our second response variable) as a function of change in broad-scale air temperature (Δ*TAMJbroad* change in April–June air temperatures between 2000–2017 and 1981–1999) between the two time periods, local temperature variability and size of the PA. In this model (GLMM), forest proportion in the PA was regarded as an offset variable because of the intercorrelation with σ*TAMJlocal*. The offset term is a structural predictor, the values of which are not estimated by the model but are assumed to have the value of one. The values of the offset are added to the linear predictor of the target. Here, we assumed Gaussian error distribution and included species as a random factor:

$$\text{Change in bird density} = \,\,\Delta TAMf\_{\text{broad}} + \sigma TAMf\_{\text{local}} + \log(\text{PA size}) + (1 \,\, \text{[species]} \,\text{m}^{-1}) \tag{2}$$

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 Δ*TAMJbroad* and σ*TAMJlocal* (rs = 0.388), σ*TAMJlocal* and size of the PA (rs = 0.312), and Δ*TAMJbroad* and size of the PA (rs = 0.208).

When comparing species-specific density changes in the PAs between the time periods, we used a non-parametric Wilcoxon signed rank test.

#### **3. 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 *TAMJbroad* 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 σ*TAMJlocal*, 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 (*TAMJbroad*), local air temperature variability (σ*TAMJlocal*), forest proportion, interaction between σ*TAMJlocal* 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 (Δ*TAMJbroad*) 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 (Δ*TAMJbroad*), local air temperature variability (σ*TAMJlocal*) 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.


#### **4. 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 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]).

#### **5. 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.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/1424-2818/12/2/56/s1, Table S1. Protected areas studied. Figure S1. The agreement between observed and interpolated monthly average air temperatures (1981–2010).

**Author Contributions:** Conceptualization, J.A., M.L., R.K.H and R.V.; methodology, A.R., J.A., M.L., N.L., R.K.H. and R.V.; data analysis, J.A., R.H.K., M.L., N.L., R.V.; data curation, A.R., J.A., N.L., R.V.; writing—original draft preparation, A.R., J.A., M.L., N.L., R.K.H, R.V. and S.K.; project administration and funding acquisition, R.K.H., R.V. and S.K. All authors have read and agreed to the published version of the manuscript.

**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).

**Acknowledgments:** Bird censuses were carried out by numerous field ornithologists, whom we gratefully acknowledge. Most of the census data have been collected in Metsähallitus, National Parks Finland.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*
