Disentangling direct and indirect effects of local temperature on abundance of mountain birds and implications for understanding global change impacts

Data collection

Data were collected as previously described in Ceresa et al. (2020a). Specifically, fieldwork took place in the central-eastern Italian Alps (Wipptal, South Tyrol), immediately south of the main Alpine watershed, within a 3,400 ha-wide area (46.96° N, 11.50° E). This area spans from approximately 1,300 m asl (just above the bottom valley) to more than 2,700 m asl at the highest peaks. At lower elevation, mountainsides are mostly covered by woodlands, dominated by spruce Picea abies and European larch Larix decidua. Above the timberline (approximately 2,000 m asl, but highly variable), wide areas are covered by bushes (mainly Rhododendron spp.) and scattered larches, whereas alpine grasslands, rocks and scree slopes characterize the upper elevation belt.

Birds were surveyed during the breeding season of most local breeding species (May–July) in 2018, by means of point counts (N = 109, Fig. 1) carried out by expert ornithologists. The rugged topography, and the consequent low accessibility of some areas, made a random selection of the sampling points unfeasible. For this reason, we carried out the point counts along accessible transects, often along footpaths; we chose the location of each point based on a minimum distance of 200 m to the previous one, in order to avoid potential replicated counts of the same individuals (Chamberlain et al., 2016).

At each sampling point, the surveyor recorded all birds observed or heard during 10 min within a 100 m radius. All points were visited three times (survey sessions: 30 May–6 June, 19–23 June and 13–18 July). The survey timing allowed excluding the main migration periods of locally breeding species. Counts were carried out in the morning, between 5:30 and 11:30, and avoiding poor weather conditions (i.e., rain or strong wind).

Vegetation and other land cover characteristics within a 100 m radius were recorded at each sampling point by the surveyor, with the aid of detailed aerial photos (as described in Ceresa et al. (2020a)). These variables included the percentage cover of tree canopy (i.e., vegetation higher than 2 m), bushes (woody vegetation lower than 2 m, mostly represented by Rhododendron shrublands), grassland (areas with no canopy and covered by grassland vegetation) and rocks/scree (emerging bedrock or scree-covered patches). The percentage cover of each variable was thus visually estimated in the field to the nearest 5%. Variables covering less than 5% of the plot surface (as defined by the 100 m radius) were assigned a 1% cover value. We also recorded information about woodland composition, by visually estimating the proportion of tree canopy cover occupied by each arboreal species (excluding rare species, i.e., representing less than 20% of the tree canopy cover).

In order to obtain climatic information at the site level, at each sampling point we placed an iButton data logger (models DS1921G and DS1922L; Maxim Integrated, San Jose, CA, USA), which recorded temperatures hourly during the entire study period. Loggers were placed at the ground level and were protected from the direct solar radiation by means of a white plastic panel. Further details about loggers’ placement and use, as well as the temperature trends recorded at each point during the study period, are provided in Ceresa et al. (2020a). As we collected such fine-scale temperature values, we did not consider topographic variables in our study. In fact, the effects of topography on bird abundance are very likely to be largely indirect, through a strong influence on microclimate (especially temperature) and, consequently, also on vegetation. Our fine-scale temperatures were indeed clearly related to both elevation and aspect, as expected (see Fig. S1). Relevant direct effects of topography on bird abundance (non-mediated by local temperature and vegetation) are indeed hardly to occur, especially for elevation and aspect.

Statistical analysis

We used structural equation modelling to disentangle the direct effect of microclimate on bird abundance and the indirect effects of local temperature, mediated by vegetation characteristics. Structural equation modelling consists in a multivariate linear regression analysis, which allows assessing the effect of a predictor on a dependent variable through mediating variables. Therefore, following this approach, it is possible to model simultaneously both direct and indirect effect of a predictor on a dependent variable (Rosseel, 2012). In the model structure we adopted, variables describing vegetation characteristics represented the mediating variables in the indirect pathways relating local temperature and bird abundance (see the path diagrams in Fig. 2). Models were fitted using the sem function of package lavaan version 0.6-7 (Rosseel, 2012) in program R version 3.5.3 (R Core Team, 2019); we used the Satorra-Bentler maximum likelihood test statistics, which provides model fit measures that are robust to data non-normality (Satorra & Bentler, 1994). We considered standardized path coefficient as significant at p ≤ 0.1, to reduce the risk of type II error associated with the adopted statistical approach (Shipley, 2016). An indirect effect was considered as significant when all single path coefficients of the indirect pathway were significant at p ≤ 0.1. We compared direct and indirect effects of local temperature by comparing the standardized coefficients of direct and indirect pathways.

To describe vegetation characteristics in our models, we used two principal components describing overall 60% of the variability in the vegetation data collected in the field (PC1 = 40%, PC2 = 20%). Principal component analysis was carried out using the prcomp function in the base R package. The first component was positively associated with forested areas and negatively with open areas (grassland and rocks), while the second component described a gradient in woodland composition (it was negatively associated with larch, one of the two dominant arboreal species), as well as the degree of grassland cover in open areas (positive association with rocks, negative with grassland; see Table 1). As local climatic predictor, we used the mean temperature calculated for each sampling point across the entire study period, because it allows describing the general climatic conditions at a location (Virkkala et al., 2008; Stralberg et al., 2009). Mean temperatures during the study period varied across the study area between 6.6 °C and 14.8 °C, and obviously showed a clear trend along the elevation gradient (Fig. S1). As a measure of bird abundance, for each species we used the maximum number of individuals detected at each sampling point during the three sampling sessions. Unlike other studies dealing with mountain species’ ecology, but not aimed at quantitative disentangling direct and indirect climate effects (e.g., Chamberlain et al., 2013; Ceresa et al., 2020a), we considered the points recorded over the complete elevational gradient for all species. This implied that, for each species, we included also areas with clearly unsuitable vegetation characteristics. Given that vegetation characteristics are largely determined by local climate, the selection of sub-samples of the sampling points based on vegetation characteristics (as frequently done in other studies with different aims), would have biased the estimation of the indirect temperature effects on bird abundance (likely leading to strong underestimation of such effects). We only considered those species detected at least once at a minimum of 20 sampling points, i.e., 15 species. These species (listed in Table 2) show a variety of different habitat preferences and specialization, as well as different nesting habits and migration strategies, therefore our results are based on a highly representative species set.

We evaluated model goodness-of-fit by a chi-square fit statistic, which indicates a significant lack of fit when p < 0.05 (Hooper, Coughlan & Mullen, 2008). We also used other fit metrics commonly adopted in structural equation modelling: Comparative Fit Index (CFI), Tucker-Lewis Index (TLI), root mean square error of approximation (RMSEA) and standardized root mean square residual (SRMR). A good model fit is indicated by CFI and TLI > 0.95, RMSEA < 0.07 and SRMR < 0.05 (Hooper, Coughlan & Mullen, 2008). Using and reporting such a variety of fit metrics is recommended in structural equation modelling, because each one reflects a different aspect of model fit (Crowley & Fan, 1997; Hooper, Coughlan & Mullen, 2008). In order to test for spatial autocorrelation, we calculated Moran’s I for model residuals using ArcGis 10.8 (ESRI, Redlands, CA, USA). Moran’s I values > 0.3 are considered as relatively large (Lichstein et al., 2002).

Direct effects of local climate for birds in the Alps

All the seven species showing significant direct temperature effects (see Table 2) are woodland or shrubland birds and are not mountain specialists. The distribution of some of these species in Italy widely overlaps with the main mountain chains, but they cannot be considered as mountain specialists because they breed also in hilly and low elevation areas, and in nearby countries they are also widely distributed across lowlands (coal tit Periparus ater, crested tit Lophophanes cristatus, Eurasian bullfinch Pyrrhula pyrrhula, goldcrest Regulus regulus, mistle thrush Turdus viscivorus; species distributions are available in BirdLife International (2021)). Given the positive direct temperature effect on their abundance, at least part of the woodland within our study area is likely colder than the thermal optimum of these species; our inner-Alpine study area is indeed characterized by an especially severe climate, with long-lasting snow cover (Adler et al., 2015). Such temperature effect could be emphasized by the early breeding of these species: warmer sites likely provide better conditions early in the reproduction season, especially for what concerns the availability of key poikilothermic prey (invertebrates). Two other species showing a positive direct temperature effect (Eurasian blackcap Sylvia atricapilla and European robin Erithacus rubecula) are generalist species with very broad habitat niche, and are very common at low and middle elevations; our study area is probably located towards the coldest extreme of their thermal niche.

In the two open habitat species for which we obtain a sufficient sample for the analysis (water pipit and northern wheatear Oenanthe oenanthe), we found a significant negative temperature effect on abundance, but exclusively indirectly via vegetation. This suggests that, at least in inner Alpine areas like our study site, maintaining and correctly managing high-elevation pastures and grasslands would probably represent effective conservation measures for these species, as proposed also in other studies (e.g., Chamberlain et al., 2013; Brambilla et al., 2018, 2020). Differently from what we expected, our results indeed did not show a direct temperature effect for the water pipit, a mountain-specialist grassland bird that had been reported to be associated with cold climates in other studies focusing on species occurrence (e.g., Brambilla et al., 2017, 2019a), and to experience variation in breeding success according to nest-site orientation (Rauter, Reyer & Bollmann, 2002). The most likely explanation for such differences is that, in this cold, inner-Alpine area, the climatic conditions of alpine grasslands are close to the thermal optimum of this species. On the other side, the previous studies reporting an effect of local climate on this species’ occurrence were carried out at a much larger scale, including also low-elevation, warmer massifs (Brambilla et al., 2017), or took place in more southern and definitely warmer mountain chains (i.e., Apennines; Brambilla et al., 2019a, Brambilla et al., 2020). The large differences in temperature ranges between our and those previously investigated study areas probably explains the differences with our results, which revealed no direct effect of temperature on abundance (besides the different focus on occurrence vs. local abundance, and the different statistical approach). This highlights the importance of spatial scale when studying the influence of climatic factors on species distribution and abundance, as different or even contrasting results may arise from analyses carried out at different scales (e.g., Franklin et al., 2013; Brambilla et al., 2019a).

While the proportion of variance in bird abundance explained by our models was good for most species (see Table 2), the proportion of unexplained variance suggests the more or less important effect of additional factors not considered in this study, such as spatial variation in predation intensity and other biotic interactions, microhabitat characteristics, and possibly also terrain slope in some species (see, e.g., Thompson, 2007; Sherry et al., 2015; Freeman & Montgomery, 2016; Brambilla et al., 2019a). We could not consider also the possible effect of precipitations on bird abundance, because we did not collect precipitation data in the field. However, over such a restricted study area the spatial variation of precipitation is likely to be affected almost only by the reduction of temperatures with elevation (through orographic effect), which should be strictly related to our fine-scale temperature data. Therefore, precipitation data would likely add only limited information to our models.

Implications for research and conservation

Our results are largely consistent with those reported by Duclos, DeLuca & King (2019), who investigated a set of forest bird species breeding in a North American mountain area and found an indirect effect of climate via vegetation in all species, while directs effects occurred only in part of the species and strongly varied in their magnitude. This suggests that such climate-species relationship patterns may be widespread in mountain bird communities, at least in temperate regions and in areas where the number of high-elevation specialists is often very low, as in our study system. This should be taken into account when investigating and predicting future bird distribution and abundance. This could be particularly relevant for generalist mountain birds, as in many cases they may be largely tracking suitable habitat, rather than thermal niches, as our results suggest for several species inhabiting different habitats and different elevation belts. Consequently, accounting also for future changes in vegetation structure and composition would likely allow more reliable predictions than using climatic factors only; this is a complex task, as both human land use and climatic factors should be considered to forecast vegetation changes. In addition, vegetation changes may lag largely behind climate change (Iverson, Prasad & Matthews, 2008; Stralberg et al., 2015): in e.g., the northeastern United States tree species shift has been forecasted to occur in some centuries (Wang et al., 2016). Such a complex and context-dependent scenario (species-specific direct and indirect climate effects, influence of human land use) may help to explain why elevational range shifts of mountain bird distribution assessed until now are poorly consistent across different studies in both magnitude and direction (i.e., uphill/downhill; see Scridel et al. (2018)), in spite of the global scale temperature increase.

While indirect climate change effects via vegetation may be strongly delayed or influenced by human activities, direct effects are likely to impact bird population rapidly, by affecting, e.g., breeding phenology, prey availability, survival and breeding success (e.g., Rodriguez & Bustamante, 2003; Dunn, 2004; Wesołowski et al., 2016; Bison et al., 2020; Ceresa et al., 2020b; Strinella et al., 2020). According to our results, at least in the innermost part of the European Alps the ongoing temperature increase may favour some generalist woodland and shrubland bird species (as suggested, e.g., by Solonen, 2005; Scridel et al., 2017). This could occur due to wider woodland areas with adequate local climatic conditions for these species, but also because of a possibly longer breeding season (Bison et al., 2020). According to Ceresa et al. (2020a), in our study area European robins select warmer areas for breeding during the first part of the reproductive season, thus higher spring temperatures would likely allow start breeding earlier across wider areas, with a consequently larger time window available for subsequent reproduction attempts.

The indirect effect of local climate we found in the two open habitat species here considered highlights the importance of conservation and management of Alpine grasslands for those species. Uphill treeline shift and shrubland expansion have been recorded in several mountain ranges, including the European Alps, and besides temperature increase they are often strongly promoted by land abandonment (Gehrig-Fasel, Guisan & Zimmermann, 2007; Myers-Smith et al., 2011). Therefore, habitat loss for these species needs to be counteracted by maintaining cattle grazing in mountain pastures, managed in an extensive way to avoid overgrazing, which leads to grassland and soil degradation (Garcia-Pausas et al., 2017) and can be detrimental to open-habitat bird species (Brambilla et al., 2020). In addition, the construction of new touristic infrastructures such as ski-pistes should be limited, because they negatively affect high-elevation grassland bird communities (Rolando et al., 2007; Caprio et al., 2011) and are predicted to increasingly overlap with the distribution of high-elevation specialist birds as a consequence of climate change (Brambilla et al., 2016). Other bird species connected to high-elevation open areas could be affected by climatic factors also directly, e.g., the white-winged snowfinch Montifringilla nivalis (Brambilla et al., 2019b), and unfortunately we did not obtain a sufficient sample for this species and other locally breeding high-elevation specialists (alpine accentor Prunella collaris, alpine chough Pyrrhocorax graculus). However, also in case of direct climate effects the aforementioned conservation measures would be recommendable, in order to try buffering against climate change. Given the high terrain complexity of mountain ranges, across large areas of suitable habitat some sites may maintain adequate microclimatic conditions for cold-associated species, possibly acting as local ‘refugia’ despite temperature warming at larger scale (Morelli et al., 2016). Therefore, avoiding habitat loss and degradation as previously described could increase the future amount of refugia maintaining both vegetation and microclimatic suitable characteristics. Furthermore, correctly managed cattle grazing could provide/maintain micro-habitats suitable for foraging for species potentially vulnerable to earlier snowmelt (Brambilla et al., 2018).