Wildlife following people: A multidisciplinary assessment of the ancient colonization of the Mediterranean Basin by a long-lived raptor

Abstract

Modern humans widely shaped present ecosystems through intentional and unintentional geographical redistribution of wildlife, both in historical and pre-historical times. However, the patterns of ancient human-mediated indirect changes in wildlife range are largely unknown, and the mechanisms behind them remain obscure.We used a multidisciplinary approach to (a) reconstruct the process of colonization of the Mediterranean Basin by a long-lived bird of prey, the Bonelli's eagle (Aquila fasciata), and (b) test the hypothesis that this colonization was unintentionally favoured by anatomically modern humans through a release of competition by dominant species, primarily golden eagles (A. chrysaetos).The fossil record of Bonelli's eagles in the Mediterranean Basin was restricted to the last c. 50 ky. This timing matches the period of modern human presence in Europe. Distribution modelling showed that Bonelli's eagles find more suitable conditions in interglacial periods, while glacial maxima are largely unfavourable unless in coastal refugia. In agreement with this, all Bonelli's eagle's fossils were found in coastal areas, and demographic inference from genetic data revealed a drop in the effective population size by around the last glacial maximum.In today's communities, we found a strongly asymmetric competitive relationship between (subordinate) Bonelli's and (dominant) golden eagles, with the former occupying far more humanized areas than the latter both at the landscape scale and the local (i.e. nesting cliff) scale. Moreover, the nesting habitat overlap analysis indicated that, in the absence of the other species, a notably higher population of Bonelli's eagle, but not of golden eagle, could be expected.Our findings are consistent with the human-mediated competitor release hypothesis, by which anatomically modern humans could have unintentionally favoured the large-scale colonization by Bonelli's eagles of a previously competitively hostile Mediterranean Basin. Reconstructing the role of ancient humans in shaping present ecosystems may help to understand the historical, current and future population trajectories of competing species of conservation concern under the ongoing scenario of global environmental change. It also illustrates how human-mediated apparent competition may promote large-scale redistribution and colonization of wildlife, including long-lived species.

Read the free Plain Language Summary for this article on the Journal blog.

INTRODUCTION

The modern human fingerprint in ecosystems through species extinctions and geographical redistribution has been pervasive since our lineage spread out of Africa (Boivin et al., 2016). Early humans contributed to precipitate the extermination of Pleistocene megafauna in several continents (Araujo et al., 2015; Barnosky et al., 2004; Dembitzer et al., 2022; Sandom et al., 2014), and we are currently driving an alarming biodiversity global crisis (Barnosky et al., 2012; Dirzo et al., 2014). Since the rise of agriculture and domestication, especially in historical times, there are numerous examples of human-caused reshaping of the distribution range of both domestic (Diamond, 2002; Zeder, 2008) and wild species (Davies, 2009; Lockwood et al., 2013). However, little is known about the patterns and mechanisms behind human-mediated indirect changes in wildlife demography and distribution that occurred in pre-historical and pre-domestication times (Boivin et al., 2016).

The Mediterranean Basin has undergone multiple and continued cultural exchanges (Valdiosera et al., 2018). Anatomically modern humans have been present in northern Africa since at least 160 ka (Smith et al., 2007), and they spread later into the European Mediterranean coast, with records since 56.8–51.7 ka for western Europe (Slimak et al., 2022). During the Pleistocene–Holocene transition, humans induced severe changes in Mediterranean landscapes, involving not only the extinction of megaherbivores (Rodríguez et al., 2004), but also the use of fire that created open habitats (Haws, 2012; Steward, 1956). Human impact on the Mediterranean Basin increased with the establishment of farming and herding societies (Naveh, 1975), which spread in Europe along an east (10.5 ka) to west (7.3 ka) wave and reached northern Africa 7 ka (Simões et al., 2023; Zeder, 2008). As a result of this ecological and socio-cultural context, human-related intentional and unintentional changes in species range have been common in the Mediterranean region in both ancient and recent times (Clavero et al., 2016; Gaubert et al., 2009; Gippoliti & Amori, 2006).

However, many past wildlife–human relationships must remain undiscovered in the Mediterranean and elsewhere—though there are signs that suggest an anthropogenic fingerprint in the present distribution of some species. This may be the case of the Bonelli's eagle (Aquila fasciata), a long-lived, large bird of prey whose presence in the Mediterranean Basin is suspected to be relatively recent. Although the status of the Bonelli's eagle during the Pleistocene in Europe has never been analysed in detail, fossils of this species in the Mediterranean Basin seem to be restricted to the period after the arrival of anatomically modern humans. In contrast, other representatives of the genus Aquila, such as golden (A. chrysaetos) and Spanish imperial eagles (A. adalberti), were present in this area well before (Sánchez-Marco, 2004). Also, the occurrence of the Bonelli's eagle in the Mediterranean is noteworthy from the geographical point of view, as the Mediterranean Basin represents the northwestern extreme of its current distribution range, which mainly falls into tropical and subtropical areas of southern Asia (Ferguson-Lees & Christie, 2001; IUCN, 2016; Appendix S1: Figure S1). Not surprisingly, this species is highly sensitive to low temperatures (Gil-Sánchez et al., 2004; Moreno-Rueda et al., 2009; Muñoz et al., 2005; Ontiveros & Pleguezuelos, 2003). The main documented competitor of Bonelli's eagle in Mediterranean regions is the golden eagle, a larger species that influences Bonelli's eagle's nest-site selection (Gil-Sánchez et al., 1996) and negatively impacts several demographic parameters (Carrete et al., 2006; Gil-Sánchez et al., 2004), including territory displacement, nest usurpation and even mortality through direct attacks (Bautista et al., 2013; Bosch et al., 2007). Both species widely segregate throughout their global distribution area, but overlapping is widespread in the Mediterranean Basin (IUCN, 2016; Appendix S1: Figure S1). All this raises a number of questions: Were Bonelli's eagles actually absent from the Mediterranean Basin before the arrival of humans? If so, what caused the establishment and spread of the Bonelli's eagle through a competitively unfavourable environment? Could humans have altered this competitive scenario?

In this study, we gather insights from different disciplines to test the hypothesis that the current distribution of a competitively subordinate species, the Bonelli's eagle, in the Mediterranean Basin is a product of the unintentional human-mediated release of competition by dominant species, especially the golden eagle. We followed a two-pronged approach, whereby we first describe the colonization pattern of the Mediterranean Basin by the Bonelli's eagle and then propose an ecological mechanism that could explain it. On one hand, we reconstructed the potential past geographical distribution of Bonelli's eagle and its main competitor in the Mediterranean Basin through spatial modelling, and used genetic data to trace the demographic history of the Bonelli's eagle in the Mediterranean Basin. This allowed us to explore the influence of climatic changes related to glacial ages on distribution and demographic patterns. We also revised the fossil record to estimate the time by which Bonelli's eagles established in the Mediterranean Basin. On the other hand, we disentangled the asymmetric competitive interactions between Bonelli's and golden eagles in relation to human presence, as a mechanism potentially behind the human-mediated competitive release hypothesis.

METHODS Modelling the present and past environmental favourability of Bonelli's and golden eagles

We modelled the distribution of environmentally favourable areas for Bonelli's and golden eagles in the Mediterranean Basin in two periods: present days and last glacial maximum (LGM; 27–19 Ka; Armstrong et al., 2019; Clark et al., 2009). We obtained the present distribution of Bonelli's and golden eagles from polygon shapefiles available from the IUCN website (The IUCN Red List of Threatened Species, version 2017.3, http://www.iucnredlist.org), projected to a 1°× 1° resolution grid system. These cells were used as operational geographic units (OGUs). Considering that the ranges of both species are well known at the examined scale, we modelled the presence and absence data. The use of absence data in well-surveyed areas, like our case, helps refining the model predictions, reducing false-positive errors and improving overall model accuracy (Chamorro et al., 2021; Lobo et al., 2010; Muñoz et al., 2005). According to previously available knowledge about the factors affecting the distribution of these eagles (Brown et al., 2017; Muñoz et al., 2013), we preselected a set of 12 environmental variables related to topography and climate to account for their favourable areas (Muñoz et al., 2015; Appendix S1: Table S1). Present and past (LGM) climatic variables were obtained from WorldClim version 1.4 (Hijmans et al., 2005). Topography is an influential factor in the distribution of mountain species (i.e. not a mere surrogate of climate; Muñoz et al., 2013), and has remained almost unchanged over the study period. Altitude was obtained from the US Geological Survey (1996), and slope was derived from a digital elevation model using the altitude variable. We assessed whether the past values of the climatic variables that entered the models of the two studied species were within the range of values that these variables had for the present.

We built models revealing the current distribution of environmentally favourable areas for Bonelli's and golden eagles by employing the Favourability Function (Real et al., 2006). To obtain favourability values (F), we used this equation:[Image Omitted. See PDF]where P is probability of occurrence, calculated through a forward-stepwise logistic regression, and n₁ and n₀ are the number of OGUs with reported presences and absences of the species, respectively. Favourability values range from 0 (minimum favourability) to 1 (maximum favourability) and refer to the environmental conditions that favour the presence of a species in a given area, discounting the effect of prevalence (Acevedo & Real, 2012).

The correlations between individual variables and occurrence areas were tested separately, as well as their response curves. To reduce multicollinearity, we explored Spearman's correlations between predictors. Whenever two variables were highly correlated within a model (|r| > 0.75), the one with the lowest contribution was excluded and the model was re-run. We assessed the discrimination capacity of the models using the area under the Receiver Operating Characteristic (ROC) curve (AUC; Lobo et al., 2008). We quantified the classification ability of the models through the proportion of correctly classified (a) presences (i.e. sensitivity) and (b) absences (i.e. specificity; Fielding & Bell, 1997).

The so-obtained models were updated for the LGM period, using the predictions of the three available global circulation models (CCSM, MIROC and MPI), on the invariant topographic variables and the projected climatic predictors, re-estimating all the coefficients.

Genetic sampling and genotyping of Bonelli's eagle

We sampled 14 subpopulations of Bonelli's eagle in the western Mediterranean Basin: France (n = 1 population), Spain (n = 11), Portugal (n = 1) and Morocco (n = 1; see Appendix S2: Table S1 and Figure S1 and Data Sources). Feathers were obtained during the period 1994–2013, either from sampled individuals (see Mira et al., 2013; Resano-Mayor et al., 2014) or at roosting sites. We only analysed one feather (i.e. one individual) per breeding territory. Genomic DNA was extracted from blood and from the superior umbilicus of the feather shaft (Horváth et al., 2005) using the QIAGEN DNeasy Blood & Tissue Kit in accordance with the manufacturer's instructions. A set of 26 microsatellites were genotyped, comprising two from Martínez-Cruz et al. (2002), 15 from Mira et al. (2005) and nine new microsatellite markers developed in this work. To search for these new microsatellite markers in Bonelli's eagle genome, a pool of 14 samples from different individuals was sent to Genoscreen, France, for microsatellite development through 454 GS-FLX Titanium pyrosequencing of enriched DNA libraries (Malausa et al., 2011). Total DNA was enriched for eight repeat motifs (AG, AC, AAC, AAG, AGG, ACG, ACAT and ATCT). Briefly, GS-FLX libraries were constructed following the manufacturer's protocols (Roche Diagnostics) and sequenced on a GsFLX-PTP. The bioinformatics program QDD (Meglécz et al., 2010) was used to filter for redundancy, resulting in a final set of 3883 sequences from which 353 primer pairs were initially designed. From those, we tested 24 markers that produced specific and reliable amplifications, but only nine were polymorphic (Appendix S2: Table S2). These were the nine new markers that have been incorporated into this work (see Appendix S2: Section S1 and Data Sources for details on DNA amplification).

Then, genetic structure and diversity analyses (see Appendix S2: Section S1 and Data Sources) were conducted to confirm all the analysed samples as belonging to a panmictic population (see Section 3). Therefore, they were considered in the subsequent demographic analyses as a single population.

Inference of past demographic processes of Bonelli's eagle using current allele frequencies

We reconstructed the demographic history of the Bonelli's eagle in the western European Mediterranean Basin through Bayesian coalescent approaches. In particular, we used Msvar 1.3 (Beaumont, 1999; Storz & Beaumont, 2002) to infer past changes in the effective size of the population as would be expected because of bottlenecks or founding events. This Bayesian coalescent method is based on the observed allele distributions and allele frequencies, assuming that a stable population with a determinate size (N₁) decreased (or increased) linearly or exponentially some time ago (T) towards the current population size (N₀). Based on this assumption and using the MCMC approach, the method allows estimating N₀, N₁ and T. Mutation rates (θ = 2N₀μ) were assumed under the stepwise mutation model of evolution (SMM), and, as for the rest of parameters, we considered wide uninformative priors under a lognormal distribution. We used two estimates of Bonelli's eagle's mean generation time (g): 9 and 14.04 years. The first estimate corresponds to the lowest value used by Mira et al. (2013); range: 9–13 years). The second one was calculated using the equation of Lande et al. (2003): g = α + (s/(1 − s)), where α is the age of first reproduction in females (4 years for Bonelli's eagle; Ferguson-Lees & Christie, 2001) and s is the expected adult survival rate (taken from Hernández-Matías et al., 2013). Each analysis was repeated three times under the exponential growth model (Appendix S2: Table S5; Data Sources), which is appropriate for testing bottleneck and founder effects (Agudo et al., 2010).

The total number of iterations in each analysis was 5 × 10⁷, with a thinning interval of 10³ iterations. We discarded the first 10% of total iterations to avoid bias in parameter estimation due to starting conditions. The distribution of the remaining data was plotted against prior distributions to see the consistency of the results over the different runs and used to obtain the lower (10%), median (50%) and upper quantile (90%) of the posterior distributions.

Fossil record of Bonelli's eagle

We reviewed the fossil record of large raptors in the Mediterranean Basin, especially the Iberian Peninsula (Spain, Portugal and Gibraltar, UK). We focused on the Quaternary, between ca. 2.6 Ma and the Bronze Age (i.e. the onset of pre-history). First, through a bibliographic search, we identified 51 Quaternary localities with remains attributed to Hieraaetus fasciatus (=Aquila fasciata) and other large raptors, mainly A. chrysaetos, in the Mediterranean Region. There are no systematic palaeornithology analyses for the African side of the Mediterranean Basin for this period. However, we obtained some avian information published within the faunal lists from relevant African and Middle East fossil sites (see Appendix S3: Table S1 and Data Sources for details on the reviewed articles).

Second, we checked for possible misidentifications of the specimens extracted from the bibliographic search, especially regarding specimens assigned to Aquila sp. Following different atlas and osteological guides (Baumel, 1993; García-Matarranz, 2013; Schmid, 1972) and some synthesis works for the studied palaeontological records (Hernández-Carrasquilla, 1993; Mlíkovský, 2002; Steward & Hernández-Carrasquilla, 1997; Tyrberg, 1998), we performed a morphometric comparative analysis of diagnostic anatomical elements to distinguish Bonelli's eagles (Aquila fasciata = Hieraaetus fasciatus) from other large Mediterranean accipitrids that might cause confusion. This allowed us to assess whether species assignment made by the authors of the reviewed studies could contain errors due to confusion with similar species. The selected diagnostic anatomical elements for this analysis were: cranial traits, sternum, coxal bone and long bones (see Appendix S3: Section S1 and Figure S1 for details). The examined material came from the collections of the Natural History Museum and the University Complutense of Madrid, both located in Madrid (Spain).

Third, following new data and recent scientific advances, such as novel dating techniques and more complete stratigraphic studies, we re-assessed the dates of those localities that included Bonelli's eagle (Aquila fasciata = Hieraaetus fasciatus) remains.

Competition between Bonelli's and golden eagles: Abundance relationships at the landscape level

We explored presence and abundance relationships between Bonelli's and golden eagles at two spatial scales: (a) provinces of peninsular Spain with the presence of nesting Bonelli's eagles (n = 32 out of 47 peninsular provinces), and (b) 13 well-delimited and well-known breeding nuclei of Bonelli's and/or golden eagles in southern Spain (autonomous community of Andalusia; c. 88,600 km²).

First, we used Pearson's correlations to explore the relationship (a) between the number of pairs of each species (extracted from del Moral, 2006, 2009), (b) between the ratio Bonelli's:golden eagle pairs and mean human footprint (HF; a proxy for human pressure, available on https://wcshumanfootprint.org/; Venter et al., 2016) and (c) between the density of Bonelli's and golden eagles (pairs/km²) and mean HF for each Spanish province. For each province, we excluded those areas that were too cold for holding Bonelli's eagle pairs (i.e. <4°C mean minimum temperature in January; see Gil-Sánchez et al., 1996), as these are areas where interspecific competition for nesting sites is prevented. Calculations of province surface and HF were performed in ArcGIS (ESRI, 2016).

Second, we analysed the relationship between the densities of nesting pairs of Bonelli's and golden eagles in the 13 breeding nuclei of southern Spain (Appendix S4: Figure S1), as well as the factors affecting such densities (Appendix S4: Table S1). Thanks to a long-term monitoring programme (Bautista et al., 2006; Moleón, 2006), we obtained accurate data on the number of nesting pairs of both species. Each area shows different and exclusive environmental features, including cliff availability and human density (Appendix S4: Table S1). Using data for 2016, the number of cliffs present within each area was calculated as the number of UTM 1 × 1 km squares that had at least one suitable cliff for nesting, that is a cliff with more than 10 m height with cavities and ledges able to support nests (Ontiveros, 1999). Some squares had more than one cliff, but they were not taken into account since the minimum observed distance of interspecific tolerance in the study area was ca. 1 km, so cliffs within the same 1 × 1 km square would be redundant. We also excluded redundant cliffs located <1 km apart at different UTM squares. We used the statistical reports of the Andalusian government (https://www.juntadeandalucia.es/institutodeestadisticaycartografia/temas/est/tema_poblacion.htm) to calculate the number of humans living in each area in 2012. We used Pearson's correlations to describe the relationship between the nesting pair densities of both species. We then used a multi-model inference approach to study the influence of human density, cliff availability (i.e. a proxy for breeding pairs carrying capacity; Newton, 1979) and the density of nesting pairs of the sympatric species on the density of nesting pairs of the focal species (first, Bonelli's eagle; then, golden eagle) in the studied breeding nuclei. Before model implementation, we checked for multicollinearity between our predictors (see above). Based on these analyses, we removed the density of golden eagle pairs (in the case of Bonelli's eagle's models) and the density of Bonelli's eagle pairs (in the case of golden eagle's models), and retained the density of humans and the density of cliffs. We also evaluated a null model (i.e. without explanatory variables) for each response variable. We used generalized linear models (GLMs) with Gaussian distributions and identity link functions. Model selection was based on Akaike's information criterion corrected for small sample sizes (AIC_c), and considered those models with ΔAIC_c < 2 to have similar support (Burnham & Anderson, 2002). For each selected model, we calculated the deviance (D²) explained (Burnham & Anderson, 2002). These analyses were done in R Studio software v4.3.2 (R Core Team, 2023).

Competition between Bonelli's and golden eagles: Local-scale segregation and competitive exclusion

We studied the local-scale segregation between Bonelli's and golden eagles in the province of Granada (c. 12,500 km²; Appendix S4: Figure S1), where an intensive monitoring programme of both species has been ongoing since 1993 (e.g. see Ballesteros-Duperón et al., 2009; Gil-Sánchez & Moleón, 2018). We built a nesting habitat selection model separately for each species, by means of binomial GLMs (logit link functions) that compared occupied versus non-occupied cliffs. We used data for a pool of 99 non-occupied cliffs (see previous section for details on cliff identification), including some cliffs that held ancient, abandoned eagle nests. We used 94 occupied cliffs from 55 golden eagle nesting pairs, and 45 occupied cliffs from 45 Bonelli's eagle nesting pairs randomly selected from the total known Bonelli's eagle population in Granada (72 pairs; Gil-Sánchez & Moleón, 2018; authors' pers. observ.). We only used 45 occupied cliffs by Bonelli's eagles because only 45 non-occupied cliffs were available once we excluded those cliffs located in areas of <4°C mean temperature in January, which is the minimum temperature required for the breeding of this species in Granada (Gil-Sánchez et al., 1996). For both eagle species, we considered those cliffs in which breeding was confirmed at least twice during the monitoring period as occupied. For neighbouring occupied cliffs, we only considered cliffs separated >2.5 km, that is the observed minimum distance between intraspecific neighbouring pairs in the study area (authors' pers. observ.). We randomly selected 94 non-occupied cliffs for the golden eagle model. As environmental predictors (none were redundant, i.e. |R| > 0.75), we selected those with the highest biological significance to both species based on our local knowledge (Gil-Sánchez et al., 1996, 2004), which are basically related to human presence and cliff accessibility (Appendix S4: Table S2). AICc-model selection and D² estimation were conducted as described in the previous section.

To estimate the nesting habitat overlap, we applied a species' model to the other species. For each response variable (the presence of nesting Bonelli's eagles and the presence of nesting golden eagles), we focused on the selected model (ΔAIC_c < 2) with the highest D². For each cliff occupied by species a, we used the parameters obtained in the model for species b and the variable values quantified for the cliff occupied by species a, and vice versa. Through this methodology, we obtained a ‘probability of reciprocal occupancy’, that is the probability of occupancy of one cliff by species b but that is actually occupied by species a, and vice versa. Then, we explored the asymmetry in nesting habitat overlap by comparing the average value of the probability of reciprocal occupancy obtained for each species. To estimate the demographic impact of habitat exclusion in an interspecific scenario of dominance, we calculated the number of potential nesting pairs of species a that we could find in favourable cliffs (threshold p = 0.7) for species a but actually occupied by species b, assuming such cliffs would be occupied by species a in the absence of species b. We took into account intraspecific effects, by excluding the nesting cliffs located at <2.5 km. For the overlap analysis, we excluded those golden eagle cliffs located out of the favourable climatic range for Bonelli's eagle (see above). Finally, we estimated the global overlap by the Cole's index (Cole, 1949), as follows:[Image Omitted. See PDF]where overlap is the Cole's index (% of overlap), est. BE is the estimated number of Bonelli's eagle pairs in the absence of golden eagles; est. GE is the estimated number of golden eagle pairs in the absence of Bonelli's eagles; BE is the total number of actual Bonelli's eagle pairs; and GE is the total number of actual golden eagle pairs.

Ethics statement

The monitoring of nesting sites of Bonelli's and golden eagles, as well as the feather collection of Bonelli's eagles, were conducted to conform to the legal requirements of competent organisms.

RESULTS The biogeographic pattern: Present and past environmental favourability of Bonelli's and golden eagles

The most important variables associated with the presence of Bonelli's eagle were those related to slope, precipitation and temperature seasonality. For the golden eagle, both slope and altitude played an important role, together with annual precipitation and temperature seasonality (see Appendix S1: Table S2 and Data Sources). Overall, models showed that the present environmental conditions are much more favourable for Bonelli's eagles (Figure 1b) than during the LGM, with the exception of several coastal refugia (Figure 1a). In contrast, environmentally favourable areas for golden eagles in the Mediterranean Basin have not changed much between the LGM (Figure 1c) and present days (Figure 1d). The past values of the climatic variables that entered the models were within the range of present values for golden eagle; in the case of Bonelli's eagle, only 2.3% of past OGUs were not within the present range for the variable Temperature Seasonality. These cells are located in the north of our study area, in places where this species is not currently present, and is not predicted to have been present in the past. The discrimination capacity was high for the models of both species (AUC = 0.951 for Bonelli's eagle and 0.892 for golden eagle). The obtained models also showed high sensitivity (SE = 0.786 for Bonelli's eagle and 0.746 for golden eagle) and specificity (0.827 for Bonelli's eagle and 0.771 for golden eagle), thus identifying with precision most presences and absences of both species.

FIGURE 1. Evidence for reconstructing the past distribution and demography of Bonelli's eagle (Aquila fasciata) in the Mediterranean Basin. In the top-left column, the present and past (LGM) environmental favourability for Bonelli's (panels a and b) and golden (panels c and d) eagles is shown. In the bottom, we show the demographic history of the Bonelli's eagle in the western Mediterranean Basin, as inferred from Msvar 1.3. Lines represent the posterior distributions (densities) of the estimated parameters obtained in six independent runs for the ancient (N1) and present (N0) effective population sizes (e), and the time (in years ago) since the population bottleneck or founder effect (f). For these analyses, we assumed generation times of 9 and 14.04 years (dashed and solid lines, respectively; see the main text for details). The wide range of uninformative priors employed in the analyses are represented by dotted lines. In the top-right column, the fossil record of Bonelli's eagle, golden eagle (A. chrysaetos) and other Aquila species (A. adalberti, A. heliaca, A. pomarina, A. rapax, A. clanga and unidentified Aquila spp. other than Bonelli's eagle) in the Mediterranean Basin is shown for the last 100 ky period (current shorelines are represented): (g) localities revised in the present study (see Appendix S2: Table S1 for abbreviations); (h) spatial location of fossil records (approximate age, when known, is represented by a colour key according to (i); (i) temporal distribution of fossil records (black circumferences and ellipses) of Bonelli's eagle (top), golden eagle (middle) and other Aquila species (bottom). Time scale is shown according to chronology (1), Epochs (2) and Lithic intervals (3).

The molecular pattern: Demographic history of Bonelli's eagle from genetic evidence

We detected a total of 111 alleles, ranging from 3 to 19 at Hf-P2E11 and Hf-C6C4 loci, respectively. Micro-Checker and Genepop indicated that no marker showed consistent signatures of null alleles across populations. Therefore, all loci were considered in the subsequent analyses. The ragged values of ΔK and the decreasing trend of lnPr(X|K) that resulted from the Bayesian clustering analyses indicated that the number of genetic clusters best fitting the dataset was K = 1 (Appendix S2: Figure S2; Data Sources). Supporting this pattern, we found similar levels of genetic diversity across subpopulations and weak levels of differentiation among them (with pairwise F_ST ranging from 0.01 to 0.09; Appendix S2: Table S6). The subsequent AMOVA corroborated the lack of genetic structure within the analysed samples by revealing that only 6% of the genetic variation was found among subpopulations and 8% among individuals within subpopulations, while the remaining 86% corresponded to individual variation (F_ST = 0.06; F_IS = 0.08; F_IT = 0.14; all p < 0.001). Therefore, for the inference of past demographic processes, we considered the analysed samples as a single panmictic population representative of the western Mediterranean population.

The reconstruction of the demographic history of the Bonelli's eagle in the western European Mediterranean Basin through Bayesian coalescent approaches reflected strong signatures of a bottleneck or a recent founder effect, as supported by the segregated posterior distributions of log(N₀) and log(N₁) (Figure 1e). The densities of these estimates, and also of log(T), were very different from the distributions of the priors used (Figure 1f). The population bottleneck or founder effect was dated in tens of thousand years, with medians of ca. 16 and 34 ka for generation times of 9 and 14.04 years, respectively. The population experienced a drop of two orders of magnitude, from the ancient (i.e. pre-founder or pre-bottleneck) effective population N₁ of 15,000–20,000 individuals to the current N₀ of around a thousand (Table 1; Figure 1e). These results are in agreement with the current number of breeding pairs estimated for western Europe (n = 1100–1200; BirdLife International, 2017).

TABLE 1 Past demographic inference of the Bonelli's eagle population in the western Mediterranean Basin using Msvar 1.3. N₀ and N₁ are the mean current and ancient effective population sizes, respectively. Time (T) represents the mean date of the change in population effective size from N₁ to N₀. Estimates are presented for two different generation times (g; see main text for details) and averaged across three independent runs for each. For each estimate, the lower (10%), median (50%) and upper quantile (90%) of the posterior distributions are shown.

	g (years)	10%	50%	90%
N0 (individuals)	9	385	894	2085
	14.04	436	981	2219
N1 (individuals)	9	5401	16,615	50,490
	14.04	7073	19,019	50,297
T (years ago)	9	5757	16,014	41,713
	14.04	13,966	34,159	85,174

The palaeontological pattern: Fossil record and establishment of Bonelli's eagles in the Mediterranean Basin

We found fossil remains attributed to Bonelli's eagle in seven Mediterranean sites: Gorham's Cave (55–48.4 and 24.4–23.2 ka; Gibraltar, UK), Figueira Brava (58.8–28.9, 41.3–18.8 or 31.6–30.2 ka, depending on the dating technique; Portugal), Devil's Tower (30 ka; Gibraltar, UK), Tell Mureybet (11–10.5, 10.5–9.8 and 9.8–8.8 ka; Syria), Toscanos (2.7 ka; Spain), Altissent (undated; Spain) and Ghar Dalam (undated; Malta; see Appendix S3: Section S2 and Data Sources for details on dating, and Figure 1g for the location of fossil sites). All Bonelli's eagle's fossil sites are located in coastal areas (see Figure 1). Thus, the age and location of the earliest remains of Bonelli's eagles in Europe matches very well with the arrival of anatomically modern humans to this continent (Slimak et al., 2022). In contrast, fossils of golden eagle and other Aquila species (including unidentified Aquila spp. other than Bonelli's eagle) were more frequent and widespread, and included numerous cases older than 50 ka (Figure 1g–i). Following Sánchez-Marco (Sánchez-Marco, 2004), the oldest golden eagle fossil was from >300 ka, and the oldest Aquila clanga and Aquila sp. (other than Bonelli's eagle) fossils were from near 2 Ma. According to our anatomical comparisons and the particular diagnostic traits of A. fasciata (see Appendix S3: Section S1 for details), misidentification of specimens in the reviewed sites can be dismissed.

The mechanism behind the patterns: Evidence for competition between Bonelli's and golden eagles, and human-mediated competitive release

In relation to the landscape scale, the number of nesting pairs of both species in the Spanish provinces with the presence of nesting Bonelli's eagles showed a triangular relationship: we found some provinces with a high abundance of Bonelli's or golden eagle pairs and other provinces with a low abundance of both species, but no provinces where the abundance of both species was high (Figure 2a). Also, we found a positive relationship between HF and the ratio Bonelli's:golden eagle pairs (ln-transformed; R = 0.411; p = 0.019; Figure 2b), indicating that there were more Bonelli's eagle pairs compared to golden eagle pairs as the landscape became more anthropized. This was because the density of golden eagles was negatively related to HF (ln-transformed; R = −0.385; p = 0.030; p > 0.05 for Bonelli's eagle; Figure 2c). With regards to the analysis of breeding nuclei in Andalusia, overall we observed a negative but non-significant relationship between the nesting abundance of both eagles (R = −0.389; p = 0.189; Figure 2d). According to the generalized linear models, the nesting density of both Bonelli's and golden eagles was in general related to human and cliff densities (Appendix S4: Table S3). However, while the Bonelli's eagle density was positively related to human density (R = 0.690; p = 0.009) and, mostly, cliff density (R = 0.899; p < 0.0001; Figure 2e; Appendix S4: Tables S3 and S4), golden eagle density was inversely related to the density of people (R = −0.753; p = 0.003) and unrelated to the density of cliffs (R = −0.187; p = 0.540; Figure 2f; Appendix S4: Tables S3 and S4).

FIGURE 2. Abundance relationships at two spatial scales in peninsular Spain. In the left column, we show the relationships between the number of nesting pairs of Bonelli's eagles and golden eagles (a), between the mean human footprint (HF) and the ratio Bonelli's:golden eagle pairs (b), and the density of Bonelli's and golden eagles (pairs/km2) and mean HF (c) in the 32 Spanish provinces with the presence of nesting Bonelli's eagles. In the right column, we show the relationships between the density of nesting pairs of Bonelli's eagles and golden eagles (d) and between the density of nesting pairs of each species and the density of humans and cliffs (e and f) in 13 breeding nuclei of southern Spain (Andalusia). Lines are shown for significant relationships.

In relation to the local scale, golden eagles occupied cliffs located in rougher areas and farther from human presence than Bonelli's eagles (Table 2). This was reflected in the models of nesting-cliff selection for each species (Appendix S4: Tables S5 and S6). Regarding the analysis of nesting habitat overlap, the average reciprocal probability of occupancy was 0.95 for Bonelli's eagle (SD = 0.12, N = 31 golden eagle's nesting cliffs within the favourable climatic area for Bonelli's eagle; Figure 3a) and 0.39 for golden eagle (SD = 0.24, N = 72 Bonelli's eagle's nesting cliffs; Figure 3b). In the absence of the other species and accounting for intraspecific effects, 23 more pairs of Bonelli's eagle (Figure 3c) and only one more pair of golden eagle (Figure 3d) could be expected. This led to a moderate global overlap (Cole's index = 37.8%; (2*(23 + 1)/55 + 72)*100).

TABLE 2 Values (mean ± SE) of the variables related to the topography and degree of anthropization of cliffs used and not used for nesting by Bonelli's and golden eagles in southeastern Spain (province of Granada). The number of cliffs analysed (n) is also shown. AltiRange: Altitude range within a 500 m radius around the nesting cliff; DistUnpaved: Distance to the nearest unpaved road; DistPaved: Distance to the nearest paved road; DistHouse: Distance to the nearest isolated house; DistVillage: Distance to the nearest village.

Species	Nesting cliffs	n	AltiRange (m)	DistUnpaved (km)	DistPaved (km)	DistHouse (km)	DistVillage (km)
Bonelli's eagle	Yes	72	291.4 ± 110.7	0.36 ± 0.26	1.61 ± 1.34	0.96 ± 0.74	2.60 ± 1.51
	No	45	209.4 ± 89.9	0.18 ± 0.10	0.84 ± 1.00	0.41 ± 0.24	2.09 ± 1.71
Golden eagle	Yes	94	344.2 ± 136.2	0.63 ± 0.56	3.18 ± 2.14	1.88 ± 1.35	4.80 ± 2.25
	No	99	261.5 ± 107.3	0.28 ± 0.24	1.25 ± 1.24	0.70 ± 0.68	3.05 ± 2.41

FIGURE 3. Effects of the potential competition for nesting cliffs between Bonelli's eagles and the larger golden eagles in southeastern Spain (Granada province). According to the average reciprocal probability of occupancy, most of the cliffs currently occupied by golden eagles are suitable for Bonelli's eagles (a), while less than half of cliffs currently occupied by Bonelli's eagles are suitable for golden eagles (b). Thus, many more potential new pairs of Bonelli's eagle could be expected in the absence of the competitor species (c) than potential new pairs of golden eagle (d). Photographs by M. Otero.

DISCUSSION The past pattern of Bonelli's eagle presence and distribution in the Mediterranean Basin

We provided evidence supporting that Bonelli's eagle has a relatively recent history in Europe, with fossils restricted to the period of modern human presence in this continent (last c. 50 ky). Distribution modelling showed that Bonelli's eagles find more suitable conditions in interglacial periods, while glacial periods are largely unfavourable except in coastal refugia, thus supporting the existence of a strong temperature limitation for this species (Gil-Sánchez et al., 2004; Moreno-Rueda et al., 2009; Muñoz et al., 2005; Ontiveros & Pleguezuelos, 2003). This is concordant with genetic inferences, which revealed a drop in the effective population size by around the LGM. According to the fossil record, Bonelli's eagles started colonizing the European side of the Mediterranean around 50 ka, that is before the LGM. Strikingly, this period was much cooler than current and preceding interglacial conditions—the coldest window between the two last interglacial periods occurred around 50–13 ka, including the LGM around 27–19 ka (Armstrong et al., 2019; Clark et al., 2009; the penultimate interglacial period took place ca. 130–115 ka; Bakker et al., 2014). Then, Bonelli's eagles could have (a) become temporarily extinct during the LGM in this area, or (b) overcome the environmentally unfavourable LGM confined to the coast until warmer conditions returned. In both scenarios, there is the possibility that the immigration of additional Bonelli's eagles could have reinforced residual populations or established new ones. Though Bonelli's eagle fossils might be recovered or identified in future archaeo-palaeontological surveys, the large number of reviewed sites without Bonelli's eagle remains—but with fossils of other similar raptors—before human arrival supports that the pattern described here may be considered as plausible.

The mechanism: The modern human-mediated competitor release hypothesis

We built upon current ecological relationships among Bonelli's eagles, golden eagles and humans to infer a plausible mechanism explaining the ancient establishment by the Bonelli's eagle in an a priori hostile area that had continuously been occupied by larger eagles for many glacial cycles (Nebel et al., 2015). The results of this study are consistent with the hypothesis that anatomically modern humans facilitated the establishment of Bonelli's eagles in the climatically favourable areas of Mediterranean Europe by releasing them from the competitive exclusion of their main competitors, notably golden eagles (see Figure 4).

FIGURE 4. The modern human-mediated competitor release hypothesis to explain the colonization of the Mediterranean Basin by the Bonelli's eagle. (a) SPECMAP oceanic chronology according to the record of ocean calcite δ18O (from Imbrie et al., 1993). Terminations of major glaciations are indicated by Roman numerals. Thick arrow represents the period of modern human presence in the Mediterranean Basin (grey: presence in northern Africa; black: presence in Europe; red: agriculture and domestication period in the Basin). Thin arrows and Arabic numerals in circles indicate scenarios depicted in (b). (b) During maximum glacial periods, environmental favourability for Bonelli's eagle is low due to low temperatures and elevated competition with similar species such as golden eagle (A. chrysaetos; Scenario 1). The approximate location of ice sheets and the southern limit of the permafrost are coloured in blue and represented by a blue dotted line, respectively. Glacial peninsular refugia are indicated by red R symbols (from Provan & Bennett, 2008). During minimum glacial periods, the ice and permafrost limits move northwards (as indicated by blue arrows), and environmental favourability for Bonelli's eagle is higher (but still relatively low) than in glacial periods due to milder weather conditions (Scenario 2). Scenarios 1 and 2 have been recurrent during the Quaternary period. Interspecific competition is reduced as modern humans spread through the Mediterranean Basin because larger Aquila species were likely persecuted and avoided areas around human settlements (black dots; Scenario 3). This competitor release would facilitate the colonization of relatively humanized patches by the Bonelli's eagle, which is less sensitive than larger Aquila species to human presence and activity. Environmental favourability for Bonelli's eagle increases as the human population grows and after the establishment of agriculture and domestication (Scenario 4).

Our ecological assessment provides a possible explanation for the establishment of Bonelli's eagle populations in Europe, despite the generally unfavourable climatic conditions that prevailed during the last glacial period. In the present day, we found that golden eagles strongly exclude Bonelli's eagles from nesting sites, with the latter being restricted to the most anthropized areas, especially in the presence of a large availability of suitable nesting cliffs that may relax competition (Gil-Sánchez et al., 1996). In fact, golden eagles are able to actively displace (and even kill; Bosch et al., 2007) Bonelli's eagles and usurp their nesting cliffs and nests (Bautista et al., 2013), while the opposite has not been documented. The great difference in the reciprocal probability of occupancy shown by both species demonstrated a strong asymmetry in nesting habitat overlap, thus revealing an important limitation of Bonelli's eagles by golden eagles within the overlapped niche and a competitive release of Bonelli's eagles in those cliffs with the highest human presence, which were largely avoided by golden eagles. This suggests that golden eagles would dominate in human-free landscapes such as those characterizing the European side of the Mediterranean Basin before 50 ka. Then, a growing human population could have locally reduced the abundance of golden eagles, particularly after the LGM (human population size in Europe, after its minimum during the LGM, reached pre-LGM levels at 13 ka; Tallavaara et al., 2015) and the onset of sedentary human societies by 11–7 ka (Zeder, 2008), thus favouring the settlement of breeding Bonelli's eagles. The competitor release effect is frequently found in modern ecosystems (Ritchie & Johnson, 2009; Segre et al., 2016) and explains many secondary or invader-facilitated invasions (Grosholz, 2005; O'Loughlin & Green, 2017). Also, the positive, large-scale demographic effects on subordinate species of the experimental reduction in the population of competitively dominant species have already been demonstrated in raptors (Sergio & Hiraldo, 2008; Wiens et al., 2021). Our findings suggest that this may also have occurred in ancient ecosystems.

Today, golden eagles and many other large eagles are highly sensitive to human presence. This may be due to recurrent persecution by humans, which usually perceive them as predators of game and domestic animals (Watson, 2010). Thus, it is reasonable to presume that golden eagles were particularly sensitive to humans and/or subject to human persecution also in pre-historic times. For the early Palearctic people, golden eagles had a high value in symbolic behaviour and as ornaments (Laroulandie et al., 2020; Watson, 2010), which could have led to further killing. Moreover, golden eagles are especially vulnerable to trapping and poisoning because they are highly attracted to carrion (Sánchez-Zapata et al., 2010; Bonelli's eagles very rarely scavenge; Moleón et al., 2009). Thus, the expansion of modern humans through the Mediterranean Basin possibly generated an increased number of patches around human settlements where golden eagles and other large eagles became rarer. These larger eagle-free patches could have been subsequently colonized by Bonelli's eagle, a competitively inferior species but notably more tolerant of human presence compared to larger eagles (Gil-Sánchez et al., 1996, 2004; Moreno-Rueda et al., 2009; Figure 4). Indeed, Bonelli's eagles display a particularly cautious behaviour (e.g. reduced soaring flight activity compared to golden eagles, and directional flights at ground level to approach their nests in the presence of humans) that helps them live in close proximity to humans while going unnoticed (authors' pers. observ.).

Caveats and alternative explanations

We should recognize several potential caveats. First, regarding the contrasting results we found in the fossil record of Bonelli's eagles compared to larger eagles, it could be argued that fossils of smaller avian species are generally found less frequently than those of larger birds (Gardner et al., 2016). However, these differences become apparent only when comparing species that substantially differ in body size (i.e. twofold or more); moreover, the fossils of smaller species are better preserved than those of larger species (Gardner et al., 2016). Thus, it is unlikely that the relatively small differences in body size between the Bonelli's eagles and other Aquila species (see Appendix S3: Section S1) are sufficient to explain the paucity and, especially, the temporal distribution of the former species' fossil record. Second, applying current ecological patterns to past scenarios may be risky. In particular, we cannot rule out that golden eagles may have had a higher tolerance for human presence with fewer people and lower levels of persecution. However, the elimination of breeding pairs by humans and the consequent abandonment of nesting sites by golden eagles likely still occurred in the past. Moreover, size-dependent competitive interactions between Bonelli's and golden eagles must have remained unchanged.

Concerning alternative explanations, the relevant question is why Bonelli's eagles were absent in the Mediterranean Basin in preceding interglacial periods, which likely were climatically favourable for this species too (Bakker et al., 2014). Alternatives may include potential shifts in prey composition and availability around 50 ka. However, Bonelli's eagles are able to kill a wide array of medium-sized prey (Ferguson-Lees & Christie, 2001). Moreover, the key prey of current Bonelli's and other eagles in western Europe is the European rabbit (Oryctolagus cuniculus; Moleón et al., 2009), which emerged in the Iberian Peninsula in the mid-Pleistocene (López-Martínez, 2008). Indeed, this prey was so abundant during both glacial (though probably restricted to coastal areas; Branco et al., 2002) and interglacial periods that it led to the super-specialization of two predators, the Spanish imperial eagle and the Iberian lynx (Lynx pardinus), and their divergence from their Eurasian relatives since ca. 1 Ma (Ferrer & Negro, 2004). This suggests that Bonelli's eagles probably found a favourable prey base during the last glacial cycles. In addition, the prevailing temperature in Europe during the arrival of the earliest Bonelli's eagles was notably more unfavourable than those characterizing interglacial periods (see above). Thus, climatic conditions and prey availability are hardly able to compete with the human-mediated competitor release hypothesis. In other words, the most obvious environmental factor that distinguished the Mediterranean European landscape in the last 50 ky in relation to preceding interglacial and glacial periods was the presence of anatomically modern humans (Blondel et al., 2010).

Finally, we still must recognize that Bonelli's eagles could have been able to settle in Europe irrespective of modern humans, which implies that competition with larger eagles could have been insufficient to prevent the spread of immigrant Bonelli's eagles. In such a case, the simultaneous arrival of modern humans and Bonelli's eagles to this continent would be a mere coincidence, though the question of why the latter arrived around 50 ka, and not before, would remain unresolved. In any case, whether or not the arrival of Bonelli's eagles was linked to modern humans, our ecological assessment suggests that humans indirectly created a more favourable scenario for the establishment of Bonelli's eagle populations in Europe.

Final remarks

The Mediterranean Basin has a long history of human-favoured changes in the distribution and demography of species (Blondel et al., 2010), including vertebrates (Eddine et al., 2020; Gaubert et al., 2009; Gippoliti & Amori, 2006) and invertebrates (Clavero et al., 2016). According to our findings, the Bonelli's eagle should probably be added to this list. However, unlike other cases in the Mediterranean and elsewhere, in which humans directly transported (either intentionally or not) exotic species to new environments (Bax et al., 2003; Carrete & Tella, 2008) or provided resources (e.g. food and shelter) that indirectly facilitated the colonization by new species (Agudo et al., 2010; Livezey, 2009; Tomasevic & Marzluff, 2017), here we suggest a mechanism by which the non-native species took advantage of the gaps that competitively superior but highly human-sensitive species left after modern humans started to populate the landscape. This example of ancient human-mediated apparent competition illustrates how humans may promote large-scale redistribution and colonization of wildlife, which encourages wider consideration of indirect effects in biogeography and invasion ecology (Bhattarai et al., 2017; Pearson et al., 2018). Also, unlike other well-known ecological cascades that were triggered by the complete extinction of megafaunal species due to overexploitation by humans (Galetti et al., 2018), Bonelli's eagles benefited from the local exclusion of their main competitors (see, e.g. Segre et al., 2016 for a present-day example). All this expands the role of ancient humans in shaping present ecosystems (Boivin et al., 2016).

Overall, we provide an example of how cross-disciplinary integration may shed light on challenging palaeo-ecological and biogeographical questions (Barrientos et al., 2014; Boivin et al., 2016; Clavero et al., 2016; Livezey, 2009; Tella, 2011). Our study could inspire other researchers to explore similar ancient, large-scale cases of non-native species spreading in the Mediterranean and elsewhere, not only driven by humans, but also by any natural enemy able to trigger a competitor release effect. Paradoxically, the initial advantages of living near humans have transformed into an ecological trap for Bonelli's eagles in contemporary times. Presently, the populations of this species in the western Palearctic are significantly jeopardized by threats such as electrocution and habitat degradation, which stem from the extensive intensification of human activities (Real et al., 2001; Viada, 2021).

AUTHOR CONTRIBUTIONS

Marcos Moleón, José M. Gil-Sánchez and José A. Sánchez-Zapata conceived the idea. Marcos Moleón, Eva Graciá, Nuria García, José M. Gil-Sánchez, Raquel Godinho, Joan Real, Antonio Hernández-Matías, A. Román Muñoz and Eneko Arrondo collected the data. Marcos Moleón, Eva Graciá, José M. Gil-Sánchez, Raquel Godinho, A. Román Muñoz and Eneko Arrondo performed the analyses. Marcos Moleón led the writing of the manuscript, and all authors critically contributed to revise it.

ACKNOWLEDGEMENTS

J.F. Sánchez-Clemot, M. Otero, G. Valenzuela, F. Molino, J. Jaramillo and E. Ávila helped during fieldwork.

FUNDING INFORMATION

M.M. acknowledges financial support through the Severo Ochoa Programme for Centres of Excellence in R+D+I (SEV-2012-0262) and a research contract with Ramón y Cajal from the MINECO (RYC-2015-19231). R.G. was supported by a research contract by the Portuguese Foundation for Science and Technology (2021.00647.CEECIND). This study was partly funded by the Spanish Ministry of Economy and Competitiveness, the Spanish Ministry of Science and Innovation and EU ERDF funds through projects CGL2012-40013-C02-01/02, CGL2013-41565-P, PGC2018-093925-B-C33 and PID2019-105682RA-I00.

CONFLICT OF INTEREST STATEMENT

The authors declare no conflict of interest.

DATA AVAILABILITY STATEMENT

Data are available from the Dryad Digital Repository: https://doi.org/10.5061/dryad.qrfj6q5q4.