Ecological and geographical marginality in rear edge populations of Palaearctic forest birds

The centre–periphery hypothesis predicts that habitat suitability will decrease at the edge of a species’ range, a pattern often questioned by empirical data. Here we explore if habitat suitability decreases southwards and shapes the abundance distribution of rear edge populations of forest birds within the restricted geographical setting of the south‐western Palaearctic. We also test if birds endemic to the area fit more poorly to the latitudinal decrease in habitat suitability due to the putative effect of adaptations to regional conditions.

and species shifts poleward, it has been suggested that rear edge populations could be extirpated from their low latitude, droughtprone Mediterranean areas of the Palaearctic, a prediction that often overlooks some particularities affecting their resilience (Vilà-Cabrera et al., 2019). For instance, rear edge populations often show taxonomic differences, suggesting reproductive isolation and putative adaptations to regional conditions (García-Ramos & Kirkpatrick, 1997;Kawecki, 2008). In addition, they occur in refuge areas from which these species have shifted northwards after the Quaternary ice ages (Hampe & Petit, 2005;Hewitt, 2000), which suggests the ecological centrality of these southern refuges and the concomitant independence of latitudinal environmental gradients (Pironon et al., 2017). Thus, the patterns of habitat suitability at the southern edge of the Palaearctic for rear edge populations, and the potential vulnerability of these populations to increasing drought and temperature (Bucchignani et al., 2018;Giorgio & Lionello, 2008), remain unclear.
North-western Africa is occupied by Palaearctic forest birds that, similar to the woodlands they inhabit, reach here the southern edge border ( Figure 1A). Forest birds are particularly abundant in the woodlands of Central Europe from where they decrease outwards (Forsman & Mönkkönen, 2003;Tellería et al., 2003). This pattern suggests that habitat suitability could decrease along northern Africa where Palaearctic woodlands and their birds reach their range edge in the Saharan expanses ( Figure 1). However, as is common in rear edge populations (see above), many north-African forest birds are taxonomically differentiated from the large European populations, a pattern usually related to the barrier effect of the Strait of Gibraltar (e.g. Doña et al., 2016;Perktas et al., 2011;Potti et al., 2016). In addition, these southern woodlands have been considered refuge areas from which forest birds expanded northwards after glacial times (Pérez-Tris et al., 2004). Therefore, it can be alternatively argued that these populations will be adapted to regional conditions and that the southward decline in habitat suitability predicted for the Palaearctic will not affect their distribution. To test these hypotheses, we will study the distribution of 11 Palaearctic forest passerines (O. Passeriformes; Table 1) by testing the following predictions:

| Geographical versus ecological marginality
Habitat suitability in north-western Africa will decrease southwards as predicted by the centre-periphery hypothesis for the whole western Palaearctic. Interestingly, the idiosyncratic north-south distribution of woodland patches along the study area ( Figure 1) provides an exceptional opportunity to explore this prediction. We agree, however, that this geographical pattern may be subtle and difficult to detect because of the small latitudinal range covered by our study and the patchy distribution of woodlands resulting from environmental constraints and/or human impact (Cheddadi et al., 2015).

| Geographical distribution of abundance
The reported patterns of habitat suitability will shape the actual distribution of forest bird abundance at the edge of the Palaearctic. This effect can be masked, however, by the difficulties to reach isolated or distant habitat patches (Soberón, 2007). In this context, the idiosyncratic north-south distribution of woodlands within the F I G U R E 1 Location of the study area and distribution of mountain ranges (a, in grey) and woodlands (b, in green) in the western Palaearctic. Distribution of sampling points (c) and line transects (d) in which the occurrence and abundance of the studied forest birds has been evaluated. The situation of the Rif and Atlas Mountains and the Strait of Gibraltar is also indicated  Figure 1B) suggests that forest bird populations could be connected through stepping stone movements from some source areas (e.g. Hannah et al., 2014;Saura et al., 2014). As woodlands around the Strait of Gibraltar show a widespread cover and diversified forest bird assemblages , it can be conjectured that this area will operate as a source area of forest birds within the region. Thus, we will also explore whether bird abundance is negatively influenced by the distance to the Strait of Gibraltar.

| Taxonomic correlates
Populations endemic to North Africa fit more poorly to the proposed large-scale latitudinal decrease in habitat suitability than populations of less differentiated species. To test this effect, we will use the subspecies and/or species endemic to north-western Africa in Clements et al., (2019) as a surrogate of their actual degree of population differentiation (Phillimore & Owens, 2006).

| Study area
The study area occupies the north-western corner of Africa ( Figure 1A) and covers a gradient that spans from Mediterranean climate areas in the north to the Saharan expanses in the south. The region is crossed by the Rif and Atlas Mountains, a huge north-south mountain range (Toubkal mountain, 4,164 m a.s.l., is the highest elevation) that strongly affects climate (it captures the rains coming TA B L E 1 Relationship between habitat suitability and distance to the Strait of Gibraltar Results show the coefficients (b ± SE) and associated t-tests (P) of the best generalized least squares mixed models (GLS) selected according to AICc in which habitat suitability indices (HS) have been regressed against distance and distance 2 to the Strait of Gibraltar. The table shows the results of two analytical approaches in which the effect of the spatial correlation has not been considered (A) or included (B).
from the Ocean) and forest cover ( Figure 1A).  Figure 1B). The most humid and populated north-western lowlands are covered mainly by cereal fields and olive groves, while the rest of the country is covered by seasonal grasslands or bare areas (Tellería et al., 2014). The entire region is dominated by a Mediterranean climate characterized by a strong summer drought.

| Bird sampling
Bird sampling was designed to work with two independent datasets recorded during the second half of April and May of 2015, 2016 and 2017, a period that encompasses the bulk of the breeding activity of the birds within the study area (Thévenot et al., 2003). A first sampling was designed to record the occurrence of the study species within a network of itineraries along small roads or tracks covering the study area. In this way, we geo-referenced all birds detected by sight and sound in 8975 waypoints ( Figure 1C) from which we extracted the occurrences of forest passerines (Table 1). These occurrences were used to model species distribution in order to use the outputs as an index of habitat suitability (Weber et al., 2016). As we arbitrarily established a minimum of 40 occurrences per species, we discarded some rare birds from further analyses (e.g. Sitta europaea n = 24, Ficedula speculigera n = 3).
To assess the effect of habitat suitability and other drivers on bird abundance (e.g. distance to the north), we used the data recorded in 190 line transects scattered within the study area (Figure1D). It is important to note here that local bird abundance results from multi scale-dependent processes strongly influenced by fine-grained effects related to local habitat structure (Seoane et al., 2004). To cope with this effect, we also assessed the local structure of vegetation We selected two components related to increasing forest cover (PC1, eigenvalue: 2.04; explained variance: 50.94%; factor loadings, grass layer: −0.06, shrub cover: 0.88, tree cover: 0.77, tree and shrub species: 0.81) and grass cover (PC2, eigenvalue: 1.01; explained variance: 25.27%; factor loadings, grass layer: 0.99, shrub cover: −0.08, tree cover −0.04, tree and shrub species: −0.08). In this way, factor scores within the components were used as comprehensive indices of habitat structure. In a previous exploratory Spearman rank order correlation analysis (rs), we detected that forest complexity was positively and significantly (p < 0.001) related to the abundance of all species while grass cover only had a positive effect on Fringilla coelebs (rs: 0.196, p = 0.007) and Turdus merula (rs: 0.182, p = 0.012), two species that usually rely on grasses and meadows. Thus, for the sake of simplicity (we aimed to work with just one variable per driver type), we selected forest complexity as the only surrogate of habitat structure.

| Species distribution model
We http://www.earth env.org/). In all cases, the selected shapes were at 30-arc-second resolution (~1 km). These variables were weakly correlated, as shown by Pearson's correlation coefficients (for all pairwise correlations, |r| ≤ 0.621; see Table S1.1). All variables were normalized (mean centring and dividing by one standard deviation).
We used Maxent (Phillips et al., 2017) to model habitat suitability and potential geographical distribution of forest birds within the study area. This model produces spatially explicit predictions based on occurrence locations and environmental data. We used a grid resolution of 1 × 1 km to remove records within the same grid (i.e. only one occurrence record per grid square of 1 × 1 km). To select the area with which to calibrate the model, we used the minimum convex polygon (MCP) produced by the full set of sampling points ( Figure S1). . The best models were selected considering significance (partial ROC, with 500 iterations and 50% of data for bootstrapping), omission rates (E = 5%) and model complexity (models within two AICc units of the minimum value among the candidate models). Partial ROC and omission rates were evaluated on models resulting from training occurrences, whereas AICc values were calculated on models created with the full set of occurrences (Warren & Seifert, 2011). We then created the best fitted final models using the full set of occurrences and the selected parameterization (Table   S1-2). Finally, we carried out an additional external validation of the best fitted models with partial ROC and similar omission rates using an independent dataset (presence and absence data collected in our line transects; Figure 1C). We used the raw output of Maxent models in further analyses because this is the output with the most straightforward interpretation (Merow et al 2013). We used 'kuenm_cal', 'kuenm_ceval' and 'kuem_mod' functions from the KUENM package to run the models (Cobos et al., 2019).

| Geographical versus ecological marginality
We generated 1000 random selected sampling points within the range of each individual species defined within an MCP covering their occurrences ( Figure S1.1). In these random points, we sampled habitat suitability scores and the distance to the Strait of Gibraltar (it was strongly related to latitude, r:0.99, p < 0.001, n = 1000). We then used three sequential approaches to show the geographical patterning of habitat suitability. First, we carried out generalized least squares (GLS) mixed models in which habitat suitability indices were regressed against distance and square distance to the Strait of Gibraltar. In this way, by selecting the best models based on AICc, we explored if habitat suitability decreased southwards monotonically or, alternatively, showed a unimodal (bell-shaped) pattern with the best sites located within an inland area. Second, we assessed the spatial autocorrelation of habitat suitability by means of Moran's index. Given the large spatial scale of our approach, this index was used to assess the geographical clumping of habitat suitability within the study area and to highlight on the potential effect of spatial autocorrelation on further analyses (Diniz-Filho et al., 2003;Dormann et al., 2007). Finally, to remove the potential effects of the spatial autocorrelation on the reported patterns, we repeated the GLS analyses by considering five different alternatives spatial correlation structures (exponential, Gaussian, spherical, linear and rational quadratic; Dormann et al., 2007). We used AICc to select the best model. These analyses were performed with the 'nlme' package (Pinheiro et al., 2020)   The table also shows the results of multivariate hierarchical Bayesian species zero-inflated distribution models between abundance forest cover, distance to the Strait of Gibraltar and habitat suitability. Asterisks show that the 95% credible interval around the parameter mean did not include zero (i.e. a significant effect).

| Geographical distribution of abundance
We performed a multivariate approach for a comprehensive view on the interacting effects of forest cover, geographical location and habitat suitability on bird abundance. After testing for multicollinearity (the highest VIF score was 2.66 in multivariate lineal models for all the species), we controlled for the effect of spatial autocorrelation and, because of the scarcity of most forest birds within the study area, we coped with the effect of a zero-inflated structure in our data. To solve both problems, we performed a zero-inflated Poisson model with an intrinsic conditional autoregressive model (iCAR) process in which we accounted for the effect of the zeroinflated structure on abundance and controlled for the potential effects of spatial autocorrelation. These steps were achieved using a hierarchical Bayesian model (Latimer et al., 2006). This type of model integrates two processes: (a) the process for which the species is present or absent in a location and (b) the process determining the number of individuals observed at suitable locations (i.e. abundance process). The first (a) process was modelled using a binomial distribution, while the second (b) was modelled 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 (Plumptre et al., 2016). For the parameter inference in a Bayesian framework, we used non-informative priors with large variance: Normal (mean = 0, variance = 10e6), except for the variance parameter of the spatial random effects, for which we used a weak informative prior: Uniform (min = 0, max = 10). We ran two parallel Markov chain Monte Carlo for each parameter. These analyses were performed with the 'hSDM' package (Vieilledent et al., 2014) in R program (R Development Core Team, 2019).

| Habitat suitability models
Species distribution models of forest birds performed well in north-western Africa according to the internal and external validations (Table S1.2). Forest/tree density displayed the most relative contribution and climatic variables reported the positive effect of precipitation on habitat suitability (Table S1.3). As a rule, the most suitable sites fitted the existing north-south distribution of woodlands (Figures 1 and 2). Within this spatial patterning, some birds displayed a wider distribution (Fringilla coelebs, Turdus merula…) while others crowded in mountain ranges (e.g. Certhia brachydactyla, Erithacus rubecula, Periparus ater, Regulus ignicapilla; Figure 2). Just one species (Periparus ater) adjusted to a unimodal distribution of habitat suitability (Table 1A) with the most suitable sites in inland mountains (Figures 2 and 3). Despite the differences, 10 of the 11 species displayed a significant monotonic north-south decrease in habitat suitability (Table 1A).

| Geographical versus ecological marginality
All species displayed spatial autocorrelation in habitat suitability (I's Moran; Table 1) as suggested by its clumped distribution within the area (Figure 3). After accounting for the effect of spatial autocorrelation, the results supported a decrease southwards of habitat suitability in all species except Periparus ater and Regulus ignicapilla (Table 1B). Thus, overall, most species displayed a significant loss of habitat suitability southwards while only two did not report any significant trend after controlling for the effect of spatial autocorrelation.

| Geographical distribution of abundance
Forest birds were recorded in only a few of the 190 randomly distributed line transects (Table 2, Figure 4), a result that agrees with the low habitat suitability reported by Maxent (Figure 3). After accounting for the effect of spatial autocorrelation and zero inflated data, multivariate approaches displayed significant patterns in all species except Troglodytes troglodytes (Table 2). Habitat suitability was positively related to abundance of all cases, with 7 of 11 species displaying a statistically significant relationship (  Figure 4). Forest cover was the least influential variable and displayed different trends according to the habitat preference of birds (we will not discuss these results as forest cover has been used just to control its local effects on bird abundance). It can be thus concluded that abundance increased with habitat suitability but displayed different patterns with distance to the Strait of Gibraltar, with most birds decreasing and just one species increasing abundances southwards.

| Taxonomic correlates.
Results in this paper do not support the prediction that birds endemic to North Africa fit more poorly to the predicted large-scale, southward decrease in habitat suitability. Six of eight endemic populations reported a significant decrease in habitat suitability after accounting for the effect of spatial autocorrelation (but see Periparus ater; Table 1). Abundance of most species tracked habitat suitability and three endemic populations displayed a negative effect of distance on abundance. Only Periparus ater (endemic) displayed the opposite trend (Table 2). It can be thus concluded that these results do not support any idiosyncratic effect of taxonomic differentiation on the distribution of habitat suitability and abundance in rear edge populations of forest passerines within the study area.

| DISCUSS ION
Understanding the factors that shape species distribution is a key issue in ecology, biogeography and conservation, but the actual input of regional studies to this field is usually flawed by several methodological shortcomings. For instance, drivers involved in species distribution may differ between spatial and environmental settings or between different model species making the results difficult to extrapolate; different methodological approaches may provide different outlooks on species distribution; and deficient data can blur the patterns on which to the predictions are tested (reviewed in Santini et al., 2019). These methodological problems have probably affected the conclusions of the most comprehensive reviews by blurring the actual relationships between habitat suitability, abundance and geographical distribution (Dallas et al., 2017;Pironon et al., 2017;Sagarin & Gaines, 2002).
Our study explores the regional trend of habitat suitability within North-western Africa to test if it decreases southwards as predicted by the Palaearctic distribution of forest birds. Our results suggest that the distribution of forest birds is affected by two idiosyncratic traits to be considered before extending our conclusions to other environmental and geographical settings. First, the latitudinal arrangement of woodlands provides a suitable scenery in which to detect the relationships between geographical and ecological marginality at the southern edge of the Palaearctic that may not occur in other areas. Second, the sharp environmental patchiness related to the presence of refuge areas suggests a meta-population dynamic in which the connectivity can strongly affect the actual distribution of populations (Hannah et al., 2014). From a methodological perspective, it is also important to point out that we used standardized field data recorded for this purpose and that we have tried an integrative approach (habitat suitability, connectivity, abundance) to clarify the effect of multiple drivers on the actual distribution of birds (Pironon et al., 2017;Sexton et al., 2009).

| Geographical versus ecological marginality
Species distribution models supported the view that the best places for forest birds were related to moist, tree-covered areas (Table S1), a pattern well-known within the xeric environmental setting of the south-western Palaearctic Tellería & Santos, 1993). As a result, habitat suitability displayed similar patterns in all  Figure 3). This species has been reported as a conifer specialist in Morocco so that such habitat selection could be related to its crowding in the extensive cedar woodlands of the Atlas Mountains (Thévenot et al., 2003). Thus, beyond one species linked to very particular habitat preferences, habitat suitability of most species decreased southwards supporting the relationship between ecological and geographical marginality predicted by the centre-periphery hypothesis (Pironon et al., 2017).

| Geographical distribution of abundance
Abundance did not display similar patterns to habitat suitability (Figures 3, 4) because under a given threshold of habitat suitability species are unable to persist without a continuous input of individuals (Pulliam, 1988). This explains why forest birds were absent in most of the study area and displayed low regional abundances (Table 2). However, over a given threshold, it can be conjectured that abundance will be positively affected by habitat suitability (Weber et al., 2016) and negatively affected by the distance to some putative source areas (Saura et al., 2014). Multivariate analyses designed to prevent the statistical effect of zero inflated data supported these conjectures. First, abundance was positively related to habitat suitability in most species (  Figure 4). This suggests that, in addition to the effect of large-scale patterns in habitat suitability, some species increased abundance in inland areas of north-western Africa. These numerical arrangements, which have usually been related to the local availability and quality of preferred habitats (Pearson & Dawson, 2003; F I G U R E 3 Distribution of habitat suitability in one thousand random distributed points within the range of the species along increasing distances to the Strait of Gibraltar. The highest scores show the most suitable areas for each individual species. Red lines show the patterns resulting from a polynomial function Seoane et al., 2004), suggest the effect of a regional meta-population dynamics on the abundance distribution of birds.

| Taxonomic correlates
Birds did not display different distribution patterns according to their taxonomic differentiation in north-western Africa (Tables 1 and 2). We acknowledge that this classification can offer a simplified view of the actual differentiation of species and that further studies could alter some aspects of the taxonomical proposal in Clements et al., (2019). For example, while the specific differentiation of Cyanistes (Illera et al., 2011) and the subspecific division of Fringilla (Svensson, 2015) has been recently tested, it seems that the African population of Sylvia atricapilla may differ from European population (Delmore et al., 2020). However, despite the uncertainties, our results do not support that endemic birds fail to adjust the predictions of the centre-periphery hypothesis. This is the case, for instance, of the African blue tit (Cyanistes teneriffae) and the Atlas chaffinch (Fringilla coelebs africana), two North African endemic birds that reported a significant monotonic loss southwards of habitat suitability (Table 1, Figure 3). Something similar can be said of abundance since all populations displayed similar patterns (Table 2). Thus, these results suggest a homogenous reaction of forest birds to the environment independently of their degree of population differentiation, a pattern that could be explained by the resilience to changes of niche in evolutionary processes (Peterson, 2011 Barbet-Massin et al., 2012;Cuervo & Møller, 2013). If habitat suitability is a main driver of distribution, it is commonly agreed that those variables with the greatest influence will represent dimensions along which the populations can be disturbed (Lee-Yaw et al., 2018).
This is the case of tree cover and precipitations, two main correlates of the distribution of forest birds in the south-western Palaearctic.
These two variables display synergetic interactions since forest cover may buffer climate effects allowing organisms to remain stable under increasingly adverse climate conditions (Suggitt et al., 2018).
Thus, the fate of forest birds will balance between woodland recovery resulting from forestry or rural abandonment (Pereira & Navarro, 2015) and woodland depletion resulting from human land use and climate-mediated processes (e.g. physiological stress, wildfires; Allen et al., 2010). Under the effect of increasing drought and temperature in this xeric region, the northwards and upwards retreat of less drought-tolerant woodlands has been predicted (Ruiz-Labourdette et al., 2012). A similar trend has been predicted for Palaearctic birds arrived to winter in this marginal area (Tellería et al., 2016(Tellería et al., , 2020. In this scenario, the relationships between ecological and geographical marginality displayed in this study support a similar retreat of forest birds in which those populations inhabiting the southernmost woodlands will be the most vulnerable. This suggests that forest remnants located in the north and in the west-facing slopes of the Atlas Mountains are a conservation priority for the maintenance of forest bird species in the region.

ACK N OWLED G EM ENTS
This paper is a contribution to the project CGL2017-85637-P (Life at the border: population differentiation of forest birds south of the Palaearctic) granted by the Spanish Ministry of Science, Innovation and Universities. REHL has been granted by Consejería de Educación de la Junta de Castilla y León and Fondo Social Europeo (EDU/556/2019). The field work consisted of bird watching on public lands where permits were not required.

DATA AVA I L A B I L I T Y S TAT E M E N T
All data supporting the results are available at https://doi.