Landscape effects on the conservation of rear-edge populations of forest passerines (Passeriformes): current patterns and prospects

Summary

Landscape changes affect species abundance and drive biodiversity loss. Here we explored if habitat amount and patch aggregation shape the abundance of forest passerines within the south-western Palaearctic (Morocco). As forests in this region are affected by increasing drought and temperature, we also forecasted their trends according to current predictions of climate change and explored how landscape changes could affect bird distribution. We recorded geo-referenced occurrences of seven forest passerines that were modelled with a set of environmental variables with Maxent to predict their distribution. The occurrence probabilities provided by the models were used as surrogates for the current distribution of habitat amount and patch aggregation within the country. In addition, 190 500-m line transects scattered within the country were used to estimate local bird abundance. Results showed that bird abundance recorded in line transects was positively correlated with habitat amount and patch aggregation of landscape around transects. This supports the idea that changes in these landscape metrics affect the abundance of the study species. Climate-change projections suggest that habitat amount and patch aggregation will decline in southern sectors but will be maintained or will increase at higher elevations. Given their relationship to abundance, these landscape changes suggest that forest birds will have to shift to the northernmost and elevated sectors. These results showed that landscape management can play an important role in the conservation of rear-edge populations of forest birds and suggest that any increase in forest amount and connectivity will improve bird resilience under a global change scenario.

Introduction

Species distributions are influenced by biotic and abiotic factors, which may operate at different scales and vary between different geographical settings (Guisan and Thuiller 2005, Sirami et al. 2017). Therefore, in the context of global change, unravelling the main drivers of species distribution will be useful for land managers and conservation practitioners in developing guidelines aimed at conserving biodiversity (Pullin and Steward 2006, Brambilla et al. 2020).

The Mediterranean region is a biodiversity hotspot affected by more frequent droughts and rising temperatures (Lionello and Scarascia 2018) in which rear-edge populations of northern organisms occur despite the effects of increasingly hostile environmental conditions (Hampe and Petit 2005, Cuttelod et al. 2008). This suggests that the current distribution of these southern populations may change or vanish in the most vulnerable sectors under climate change (Hampe and Jump 2011, Merker and Chandler 2020). This pattern is expected for forest birds, a set of species dependent on woodlands in the core of the western Palaearctic that, according to predictions, will be increasingly scarce in the southernmost ranges (Schippers et al. 2011, Barbet-Massin et al. 2012). However, the actual environmental constraints affecting the distribution of rear-edge populations are far from fully understood (Vilà-Cabrera et al. 2019). Forest birds seem to be distributed according to the centre–periphery hypothesis (reviewed in Pironon et al. 2017), which suggests a progressive deterioration of habitat suitability at the periphery of their range (Tellería et al. 2003, 2021). In this context, because of the idiosyncratic north–south distribution of woodland patches within this area (Figure 1A), it has been suggested that the rich bird assemblages of the northernmost woodlands around the Strait of Gibraltar (Correia et al. 2015) could operate as a regional source of individuals moving southwards (Tellería et al. 2021). This suggests meta-population dynamics may be occurring in this system and the resilience of these populations to climate change could be linked to their ability to track the spatial patterning of suitable habitats (Hannah et al. 2014, Mair et al. 2014), a process positively related to the size and spatial adjacency of available habitat patches (Taylor et al. 1993, Hanski 1998).

In this paper, we explore the effect of the spatial configuration of habitat patches on the abundance distribution of seven forest passerines during the breeding period in the south-western Palaearctic (Morocco) and the way these patterns could vary as a result of climate change. The selected species occur in tree-covered habitats where they mainly feed on foliage (Common Firecrest Regulus ignicapilla, Western Bonelli’s Warbler Phylloscopus bonelli, African Blue Tit Cyanistes teneriffae, and Great Tit Parus major), trunk (Short-toed Treecreeper Certhia brachydactyla), and ground invertebrates (European Robin Erithacus rubecula and Common Chaffinch Fringilla coelebs) during the breeding period (e.g. Talavera and Tellería 2022 for a quantitative approach to their feeding substrata in Mediterranean woodlands). The African Common Chaffinch has the widest habitat niche as it occurs in a broad variety of woodland landscapes, while the rest tend to occur in denser, tree-covered forests (Hanzelka et al. 2019).

We addressed our objectives through two complementary approaches.

Landscape effects on bird abundance

First, we assessed the effects of the current landscape structure on bird abundance in the south-western Palaearctic. More explicitly, we tested the effect of habitat amount (HA) and the spatial adjacency of habitat patches on the observed abundance of individual species. Our approach was as follows: 1) we used geo-referenced occurrences of each study passerine recorded within the study area to generate species distribution models (SDMs) (Elith and Leathwick 2009), and used the resulting occurrence probabilities as surrogates of habitat suitability on which to assess HA and aggregation index (AI), an index of spatial adjacency (McGarigal et al. 2002); 2) we used a set of line transects to test, once the effect of the fine-grained forest structure of sample sites was accounted for in the model, the positive effect of these two landscape metrics on bird abundance.

Prospects for change

If our results suggested that HA and AI had a positive, significant effect on bird abundance, we would forecast changes in these landscape metrics under a climate-change scenario. More explicitly, we explored: 1) how HA and AI may vary in the future, and 2) how these predicted changes may vary with elevation and distance to source areas in the north, the two main gradients related to forest bird distribution within the region (Tellería et al. 2021). Thereby, we aimed to detect the most vulnerable areas for which we propose management guidelines to conserve rear-edge populations of forest birds.

Methods

Study area

The study area spans the Mediterranean climate areas in the north to the dry Saharan expanses in the south. The Atlas Mountains (Toubkal Mountain, the highest elevation: 4,164 m a.s.l.) cross the region from north to south, creating a rain shadow and strongly affecting climate and forest cover (Figure 1A and B). The mountain range divides the region into a moist sector in the west suitable for woodlands, and a dry sector in the east where tree cover is almost non-existent. The western side is covered by humid woodlands in the north (Quercus suber, Q. canariensis, and Q. pyrenaica), cedar forests (Cedrus atlantica) and fir (Abies maroccana) in mountainous areas, drought-tolerant woodlands (Q. ilex, Juniperus thurifera, and Tetraclinis articulata) in lowland areas, and argan (Argania spinosa) and acacia (Acacia spp.) open woodlands in the driest sectors of the south (Quézel 1983). Woodlands have been strongly affected by human use, particularly in lowland agricultural areas in the west, resulting in a fragmented distribution (Cheddadi et al. 2015).

Figure 1. (A) Tree cover in the western Palaearctic and the location of the study area. (B) Spatial pattern of tree cover and elevation over 1,500 m asl within the study area. (C) Projected changes in future habitat suitability of seven forest bird species. The range shifts for all species are shown as: gained (pink), stable (grey), and lost (purple) surface areas. Suitability maps are average calculations of the two periods (2050 and 2070) under two emissions scenarios (RCP 4.5 and RCP 8.5).

Bird sampling and local habitat structure

We conducted fieldwork during the second half of April and May of 2015, 2016, and 2017. These months encompass the breeding season of passerine birds (Thévenot et al. 2003). We designed this fieldwork to produce two independent data sets that have already been described in Tellería et al. (2021); in https://datadryad.org/stash/dataset/doi:10.5061/dryad.34tmpg4kh. We used the first data to obtain geo-referenced occurrences of birds that were detected by sight and sound along small roads or tracks covering Morocco. In this way, we recorded birds at around 9.000 waypoints from which we extracted the occurrences of the seven species (Figure S1). We used these occurrences to model species distribution to develop surrogates of habitat suitability (see below). The second data set resulted from bird counts along 500-m line transects scattered within the study area. Within each line transect, we counted all birds within a 500-m belt during the early morning. In addition, we set two 25-m radius circles distributed at 150 m and 350 m from the starting point of each transect. In each circle we visually assessed grass layer, shrub (vegetation between 0.5 m and 2 m height) and tree (vegetation >2 m height) covers, and the number of tree and shrub (>0.5 m) species. We used these variables to perform a principal component (PC) analysis to determine the latent variables useful in describing local habitat structure. Finally, we selected the first component related to increasing forest complexity (FC) vs. grass layer (PC1, eigenvalue: 2.80; explained variance: 46.63%; factor loadings, grass layer: -0.37, shrub cover: 0.72, tree cover: 0.74, shrub species: 0.69, tree species: 0.78) that was used as an index of FC within each sample site. In this way, line transects provided us with information on bird abundance and forest structure.

Species distribution modelling

We used bird occurrences to explore the spatial pattern of habitat suitability by modelling the distribution of individual species using Maxent version 3.4.4 (Phillips et al. 2017). To do this, we selected 30 arc-second resolution (~1 × 1 km grid) shapes of environmental variables that were useful to describe the most productive, tree-covered areas selected by forest birds in the xeric environmental setting of the south-western Palaearctic (Tellería et al. 2021). This scale is suitable for working with these birds as, according to the observed densities, 1 km² hosts a mean of 63.2 (min–max: 0–880) and 4.8 (0–40) individuals of the densest (Common Chaffinch) and scarcest (Western Bonelli’s Warbler) species (Tellería et al. 2021). We used climatic traits related to the current strength of the Mediterranean drought (i.e. Bio10: mean temperature of warmest quarter, Bio12: annual precipitation, and Bio17: precipitation of driest quarter; CHELSA http://chelsa-climate.org/) (Table S1).

We forecasted the effects of climate change by using the average projection of these climate variables for two periods (2050 and 2070), two different emission pathways, i.e. representative concentration pathway (RCP) 4.5 and 8.5, and two global circulation models (CCSM4 and HadGEM2-ES (Gent et al. 2011, Martin et al. 2011).

To include forest cover in the bird models, we modelled the current and predicted distribution of tree-covered areas using Maxent (see above). In this case, we used the global tree density map data set (Crowther et al. 2015) as a source of tree occurrences. We selected climate variables that have been described as constraining woodland development in the Mediterranean region (Bio 4: temperature seasonality, Bio10: mean temperature of warmest quarter, Bio12: annual precipitation, and Bio17: precipitation of driest quarter from CHELSA) (Table S5). In addition, we included in the model some soil variables that affect tree distribution such as bulk density, clay, pH, etc. (Table S3) from the ISRIC – World Soil Information database (https://www.isric.org/explore/soilgrids). We modelled the predicted changes in forest distribution with the climate variables resulting from the aforementioned emission pathways and global circulation models.

To select an area in which to calibrate bird models, we used the minimum convex polygon (MCP) resulting from all the occurrences recorded per species (see Figure S1). For each species, we used the kuenm R package (Cobos et al. 2019), in R version 3.6.1 (R Development Core Team 2019), to select the best model provided by Maxent of a series of candidates arranged by different combinations of parameter settings (Merow et al. 2014). Herein, we created 119 candidate models by combining the whole set of independent variables, 17 regularisation multiplier values (0.1–1.0 at intervals of 0.1, 2–6 at intervals of 1, and 8 and 10), and all the seven possible combinations of three feature classes (linear, quadratic, and product). Best models were selected by considering significance, i.e. partial receiver operating characteristic curve (ROC), omission rates (E = 5%), and model complexity, i.e. models within two Akaike’s information criterion corrected for small sample sizes (AICc) units of the minimum value among the candidate models. We then created the best-fitted final model using the full set of occurrences and the selected parameterisation (see Table S2). Finally, we carried out an additional external validation of the best-fitted models with partial ROC and similar omission rates using an independent data set. For the analyses, we used the cumulative output because this is the output recommended when the interpretations are related to range boundaries of species (Merow et al. 2013).

We estimated the future habitat suitability by projecting the fitted models on to future projections of environmental variables for two time periods (2050 and 2070) under different emission pathways (RCP 4.5 and 8.5). We quantified the degree of change in suitable habitats for each future scenario relative to the current situation, classifying habitat as gained, unchanged, or lost. To identify suitable habitats from unsuitable ones, we first converted model output into binary maps (unsuitable: 0 and suitable: 1) based on the equal training sensitivity and specificity cumulative threshold. This threshold maximises the agreement between observed distribution of species and predicted distribution calculated by the model (Liu et al. 2005).

Landscape metrics

We assessed the effect of the spatial configuration of habitat landscapes for each individual species by considering the 1 × 1 km grids surrounding the 500-m line transects within a radius of 5,000 m from the centre of each line transect. This radius size was selected to be within the usual range for landscape studies (e.g. Saura 2021). Within this circle, we calculated HA in km² that defined the extent of the suitable habitat for each study species. We also obtained the AI of suitable habitat patches for each species under current and future climate conditions. AI is a metric of spatial adjacency and expresses the frequency at which different pairs of habitat patches appear side-by-side in a landscape. This landscape metric ranges from 0, when the patches are fully disaggregated, to 100 when patches are maximally aggregated into a single, compact patch (McGarigal et al. 2002). Both metrics express different landscape features but often show correlated patterns since aggregation decreases as HA diminishes (Wang et al. 2014). All metrics were calculated using the landscapemetrics package in R (Hesselbarth et al. 2019).

Current tracking by birds of landscape patterns

We tested the effect of HA and AI on bird abundance provided by line transects after accounting for the local effect of FC. In these analyses, we controlled for the effect of spatial autocorrelation, and because of the scarcity of forest birds within the study region, we accounted for the potential effect of a zero-inflated structure in our data. We solved both problems with a Zero-inflated Poisson model with an Intrinsic Conditional Autoregressive (iCAR) model. We conducted these analyses using a hierarchical Bayesian model (Latimer et al. 2006). This model integrates two processes: 1) the process for which the species is present or absent in a location (i.e. suitability process), and 2) the process determining the number of individuals observed at suitable locations (i.e. abundance process) to determine whether the species is present in a location and to quantify the number of individuals. We modelled the first process using a binomial distribution, while we modelled the second one with a Poisson distribution. Both processes included an iCAR model which assessed the spatial configuration of the eight nearest neighbouring cells to measure the spatial autocorrelation. The iCAR component assumes that the probability of the species’ presence and abundance at one site depends on the probability of the species’ presence and abundance at the eight nearest neighbouring 1 × 1 km cells (Plumptre et al. 2016). These analyses were performed with the hSDM package (Vieilledent et al. 2014) in R.

Prospects for change

To explore the landscape shifts related to climate change, we generated 500 random sampling circles of a 5,000 m radius within the range of each individual species. This range was defined within a multispecialty community provider (MCP) covering their occurrences (Figure S2). In these circles, we calculated the HA and AI of suitable habitat under current and projected conditions (see above). We evaluated the changes with a Generalised Least Squares (GLS) model, with HA and AI being the dependent variables, period (current and future) a fixed factor, and elevation and distance to source forests in the north the covariates. In these analyses, we were particularly interested in the period–elevation and period–distance interactions, as they highlight the potential changes in landscape of two main geographical drivers of bird distribution in the study area (Tellería et al. 2021). To remove the potential effects of the spatial autocorrelation on the GLS, we considered five alternative spatial correlation structures (exponential, Gaussian, spherical, linear, and rational quadratic) and used AICc to select the best model (Dormann et al. 2007). These analyses were performed with the nlme package (Pinheiro et al. 2020) in R.

Results

Spatial patterning of habitat suitability

Species distribution models reported by Maxent suggested that forest cover and precipitation had the largest effect on the distribution of suitable habitats for forest passerines (Table 1). This means that the most suitable areas were in the north and along the moist western slopes of the Atlas Mountains (Figure 1C). When mapping the projected distribution of habitat suitability according to climate-change projections, the results suggest that forest passerines will maintain most of their current ranges but will retreat in some lowlands and expand in some highlands (Figure 1C). This is particularly so for the Common Chaffinch and the African Blue Tit, two abundant and widespread birds within the study area (Table 2).

Table 1. Estimates of the relative contribution of environmental variable models predicting habitat suitability for forest bird species. The values represent the percentage contribution importance of each variable in the model. Percentage contribution indicates the change in regularised gain by adding the corresponding variable. Values are averages and SDs over 10 replicate runs. Bio10, mean temperature of warmest quarter; Bio12, annual precipitation; Bio17, precipitation of driest quarter; Forest: forest cover. Symbols in parentheses show the trend of the response curves for the variables: + increase; − decrease; Ω hump-shaped; = no trend.

Table 2. Distribution of mean abundance of species based on results from 190 line transects (n = number of line transects with the particular species detected) and the results of multivariate hierarchical Bayesian species zero-inflated distribution models between abundance and local tree cover within the 500-m sampling plot and the potential combinations of habitat amount (HA) and aggregation index (AI) of suitable habitat within a radius of 5 km around the sampling point. Asterisks indicate that the 95% confidence interval around the parameter mean did not include zero (i.e. a significant effect).

Landscape effects on bird abundance

The two landscape metrics (AI and HA) were positively correlated in all species suggesting high aggregation of habitat patches in landscapes with a high percentage of suitable forest (Table 2). Abundance of all birds except for the African Blue Tit showed a significant positive relationship with FC at the study sites (Table 2), but once this effect was accounted for, all species had a significant positive relationship with AI and five out of the seven species were significantly positively correlated with HA (Table 2).

Prospects for change

The scores of both landscape metrics (HA and AI) for currently existing distributions decreased southwards in most species and showed fewer discernible increases along elevation (Figure 2). Climate-change projections suggested that HA and AI will show different trends along elevation and distance to the north (see period–elevation and period–distance interactions in Table 3). According to climate projections, HA showed a parallel, non-significant decrease southwards while AI experienced a significant decline in half of the investigated species (Table 3, Figure 2). In the case of elevation, projections suggested that both landscape metrics will experience similar increases in highlands, with significant elevation–period interactions in two species in HA and six species in AI (Table 3). In all species, these projections predicted the depletion of both parameters in lowlands, while they were predicted to maintain or increase their values in highlands (Figure 2).

Figure 2. Relationship between habitat amount (HA) and the aggregation index (AI) with elevation and the distance to the Strait of Gibraltar. Results for two scenarios for the time period are also shown, with blue for current and pink for future time periods.

Table 3. Geographical configuration of suitable habitats. Results are for the coefficients (b ± SE) of the best Generalised Least Squares (GLS) mixed models selected according to Akaike’s information criterion corrected for small sample sizes (AICc), in which the habitat amount (HA) and aggregation index (AI) have been regressed against elevation, distance to the Strait of Gibraltar, and period (current and future conditions). The data have been recorded in 500 random selected circles distributed within the range of species (see Methods). The interactions of period, elevation, and distance have also been calculated to detect potential changes. Significance ^P <0.10, *P <0.05; **P <0.01; ***P <0.001.

Discussion

Habitat and climate change are the two main drivers of the worldwide loss of biodiversity (Sala et al. 2009). However, even though the damaging effects on biodiversity produced by landscape shifts can be exacerbated in some scenarios due to climate change, both processes have been primarily addressed independently (Travis 2003, Opdam and Wascher 2004, Mantyka-Pringle et al. 2012, Sirami et al. 2017, Guo et al. 2018, Pettorelli et al. 2021). This appears to be the case for forest organisms at the south-western edge of the Palaearctic, where woodlands resulting from post-glacial climate changes (Hewitt 1999) and deforestation (Cheddadi et al. 2015) are now coping with the effects of ongoing climate warming (Carnicer et al. 2011, Ruiz-Labourdette et al. 2012). In this context, our results provide some evidence of the synergetic interactions of these two drivers that could exacerbate the difficulties of forest birds to remain in the southern edge of their range.

Effects of landscape structure on bird abundance

We found a positive relationship between bird abundance and HA and AI after accounting for the effect of the fine-grain habitat structure of the sampling sites (FC; Table 2). This relationship is usually explained by the fact that large patches of suitable habitat harbour more individuals than small ones, and because the high connectivity of patches improves the recovery of populations after local extinctions (Hanski 2015). We acknowledge that there is a debate on the relative contribution of these landscape metrics since it has been suggested that in fragmented landscapes HA is the actual driver of species distribution (e.g. Fahrig 2013, 2015, Saura 2021, Basile et al. 2021). However, our study has not been designed to disentangle this hypothesis, and we merely show the strong collinearity of both parameters and their effect on bird abundance (Table 2). Thus, as these relationships between HA and connectivity may differ among landscapes (e.g. Wang et al. 2014), our results showed the idiosyncratic landscape patterns of a region where woodlands preferentially occur in northern and mountainous areas and show an increasingly fragmented pattern moving outwards from those areas (Figures 1B and 2). These patterns of landscape depletion emphasise the view provided by the observed decrease southwards of habitat suitability of most forest birds (Tellería et al. 2021).

Prospects for change

Climate-change projections suggest that the two landscape metrics (HA and AI) that affect forest bird abundance will experience different trends along elevation and distance to source areas in the north (Figure 2). The loss of HA and connectivity in southern areas forecasts a northward retreat of the range edge in several species, which is a usual prediction on the effect of climate warming in temperate regions (Lenoir and Svenning 2015). In this context, elevation will play a different role because, despite the decrease in HA and connectivity in lowlands, projections show that some species will maintain or even increase these landscape metrics in elevated areas (Figure 2). Habitat shifts upwards in elevation related to mountain ranges have been detected in several species within the south-western Palearctic (Wilson et al. 2005). However, beyond the usual withdrawal and reduction of HA in mountain tops (Sekercioglu et al. 2008), our results suggest that the hypsographic configuration of north-western Africa (with large expanses above 1500 m) (Figure 1B) will allow areas to expand suitable habitats (Figures 1C and 2). A similar projection has been reported for wintering ranges of some passerines in the south-western Palaearctic, which will expand current ranges upwards as winter temperatures increase (Tellería et al. 2016). Thus, it suggests that shifts of suitable habitat for birds along the western mountain slopes could improve the resilience of bird populations to climate change, a process that explains the role of these mountains as refuge areas during the major climatic fluctuations of the past (Molina-Venegas et al. 2017).

Conclusions

Our study highlights the importance of landscape management in the conservation of rear-edge populations of forest passerines under a scenario of global change. We acknowledge, however, that the actual consequences of the projected environmental changes are contingent and difficult to predict as they rely on synergetic interactions between human decisions affecting woodland structure, the shortcomings of climate predictions, and the idiosyncratic reaction of species to changes at range edges (Angert et al. 2020, MacLean and Beissinger 2017, Northrup et al. 2019). Despite this, our results suggest some basic guidelines aimed to protect forest birds and other forest organisms (Larsen et al. 2012) at the southern edge of their range. First, in a region strongly affected by human land use, the prevention of forest loss and fragmentation (with the concomitant decrease of habitat amount and connectivity) seems the most obvious approach to protect forest organisms. Second, the uneven geographical patterning of landscape structure advises a preventive/proactive management of southern and lowland woodland landscapes to increase cover and patch connectivity. As in other countries of the south-western Palaearctic, this recovery of woodlands could be assisted in the future by the human abandonment of less productive, marginal areas (Pereira and Navarro 2015). Finally, the large mountain ranges involved in these processes can increase forest cover and connectivity in the most elevated areas. This means that any approach to improve tree expansion in highlands could improve the regional resilience of forest organisms to the pervasive effects of climate warming in this peripheral area of the Palaearctic.